Introduction

Microbial interactions are highly relevant in oenology, as they shape the biodiversity of the fermentation process and, consequently, the organoleptic characteristics of wine (De Gioia et al., 2022). Within this microbial diversity, certain spoilage yeasts can cause sensory deviations and negatively impact wine quality (Mehlomakulu et al., 2014; Rojo et al., 2015; Kuchen et al., 2019). Traditionally, sulphur dioxide (SO₂) has been used in winemaking to control these undesirable populations; however, its use is increasingly discouraged due to its adverse effects on human health (Vally et al., 2009; Relaño de la Guía et al., 2025).

Biocontrol, understood as the ability of one organism to prevail, dominate, or biosuppress another, has long been proposed as an alternative strategy (Baker & Cook, 1974) and remains highly relevant today, particularly as a means of reducing the use of SO₂ in winemaking (Berbegal et al., 2017; Canonico et al., 2023a; Canonico et al., 2023b). It is an anthropic approach based on natural microbial interaction mechanisms (Boynton, 2019) and is designed according to the technological objectives of the process. For example, in the context of wine fermentation, Wickerhamomyces anomalus has been successfully applied against Zygosaccharomyces rouxii during the pre-fermentative stage (Kuchen et al., 2023a; Kuchen et al., 2024). Other species, such as Metschnikowia pulcherrima and Torulaspora delbrueckii, have also been proposed for use in the same phase (Canonico et al., 2023a; Canonico et al., 2023b). Nevertheless, biocontrol yeasts should not alter the conditions required for the establishment and fermentation kinetics of the fermentative yeast (Albertin et al., 2017). Achieving this dynamic balance between both species is crucial to combine biocontrol efficacy with the technological stability of the wine.

Among the factors that promote a favourable environment for biocontrol, the priority effect stands out, because it describes how the order and timing of species establishment shape subsequent interactions (Stroud et al., 2024). In this sense, the sequential co-inoculation of yeasts, such as Wickerhamomyces anomalus and Saccharomyces cerevisiae, represents a controlled application of the priority effect. By adjusting the inoculation timing and population ratios, it is possible to guide competitive interactions and define the outcome of the process, promoting biocontrol, fermentation stability, or the formation of desirable metabolites (Gallo et al., 2024). This strategy has been extensively studied in oenology and remains an active field of research (Toro & Vázquez, 2002; Su et al., 2024). However, sequential co-inoculation may present challenges related to interspecies competition, which could slow down fermentation, if not properly managed (Albertin et al., 2017). These interactions can also be modulated by physicochemical factors and interventions in the fermentation process, such as temperature and the presence of SO₂, which directly influence growth rate and the maximum population density achieved (Savadó et al., 2011; Baiano et al., 2012; Balsa-Canto et al., 2020). Such factors affect resource competition and the production of primary and secondary metabolites with inhibitory potential (Pommier et al., 2005; Mazzuco et al., 2019; Kuchen et al., 2023b; Georgescu et al., 2024), as well as the stability of these compounds in the medium and the sensitivity of the sensitive yeast (Kuchen et al., 2021).

The modelling of microbial interactions represents the first step in applying process engineering to biotechnology and the wine industry. Although grape must constitutes an ecologically complex system (Lax & Gore, 2023), and more general frameworks are sometimes required to describe microbial communities (Ruiz et al., 2023), this study aims to quantitatively elucidate the interaction mechanisms that emerge during sequential co-inoculation within the interaction dynamics. Mechanistic primary models can be used to determine the magnitude of interesting response variables, such as the maximum specific growth rate or the time required to reach a given level – in this case, one or more interacting populations (Schoener, 1973; Gilpin & Ayala, 1973). However, the parameters growth rate, carrying capacity and interspecific competition can depend on environmental conditions, and secondary models describe this dependence (McMeekin et al., 2002), thereby enabling predictions under untested conditions (Hellweger et al., 2016; Balsa-Canto et al., 2020). Evaluating yeast growth kinetics while simultaneously considering multiple physicochemical factors can be both laborious and costly. Therefore, the use of experimental designs that minimise the experimental workload while maximising the information yield is particularly relevant (Pedrozo et al., 2024).

Few studies have explored sequential co-inoculation in winemaking through ecological interaction models, and those available focuses exclusively on Saccharomyces species (Balsa-Canto et al., 2020). To our knowledge, mathematical modelling approaches have not yet been applied to investigate the timing of co-inoculation and other physicochemical factors between non-conventional yeasts and S. cerevisiae with the aim of extending biocontrol effects during the early stages of winemaking without affecting the fermentation process itself.

The objective of this study was to evaluate the relevant physicochemical factors and the interaction kinetics of sequential mixed co-inoculation between the biocontrol yeast W. anomalus and S. cerevisiae from an ecological modelling perspective. Specifically, two main goals were pursued: to maximise the presence of the biocontrol yeast in the medium and to ensure that the fermentative kinetics of S. cerevisiae would remain unaffected.

Materials and methods

1. Experimental design

A Box–Behnken design (Fisher, 1935) was employed with three replicates at the central point, totalling 15 experiments. This design allowed the simultaneous evaluation of three factors: temperature (15, 17.5, and 20 °C), SO₂ concentration (0, 0.1, and 0.2 ppm of molecular SO₂), and sequential co-inoculation time (0, 24, and 48 h). The response variables analysed were the growth kinetics (for 120 h) of the two yeast populations involved: the biocontrol strain Wickerhamomyces anomalus BWa156 and the fermentative strain Saccharomyces cerevisiae V22 previously identified at the molecular level (Kuchen et al., 2019; Petrignani et al., 2024).

2. Experimental procedure

Pre-inoculum: the yeast strains were activated on YPD-agar. Pre-inoculum of both strains were prepared in 250 mL Erlenmeyer flasks containing 200 mL of concentrated/diluted must adjusted to 21 °Brix and supplemented with 1 % yeast extract. Cultures were incubated for 24 h at 110 rpm and 25 °C.

Fermentation trials: the 15 fermentations were conducted in 250 mL Erlenmeyer flasks containing 200 mL of concentrated/diluted must adjusted to 21 °Brix and supplemented with 0.1 % yeast extract. The initial inoculum for both strains was set at 1.5 × 10⁶ cells/mL. Sampling was performed after 24, 48, 72, 96 and 120 h respectively, in addition to the initial inoculation of both yeast strains; i.e., there was a total of six sampling time points per experiment. Given there were 15 experiments and two yeast species, a total of 180 data points were collected. Yeast counts were conducted on WLN agar plates using appropriate dilutions to ensure colony counts of between 50 and 250 CFU per plate.

3. Mathematical modelling

Two models based on the Lotka–Volterra conceptual framework were proposed to predict the population density of both species in mixed cultures. This model describes intra- and inter-specific competition (Murray, 2001). Intraspecific competition is defined as competition among individuals of the same population for the exploitation of resources or space (Schoener, 1973), while interspecific competition refers to facilitative or inhibitory interactions between populations (Toju et al., 2018). However, Gilpin and Ayala (1973) highlighted the limitations of the Lotka–Volterra framework when dealing with nonlinear patterns of intraspecific regulation.

All these models can be embedded in the following mathematical formulation:

$\frac{dxi}{dt}=G\left(xi,xj\right)=x_i{\mu }_i*\left(1-{\left(\frac{x_i}{K_i}\right)}^{{\theta }_i}-{\alpha }_{i,j}\left(\frac{x_j}{K_i}\right)\right) i,j= Wa,Sc i\ne j \left(\mathrm{E}\mathrm{q}.\mathrm{ }1\right)$

where i is either S. cerevisiae (Sc) or W. anomalus (Wa) and j either Wa or Sc with i ≠ j; $x_i$ corresponds to the cell density for both species; μ_i corresponds to the specific growth rate; K_i is the carrying capacity for species i; and θ_i ≥ 1 controls the degree of non-linearity in intraspecific growth regulation. Coefficients α_i,j measure the competitive strength of species j on i. If α_i,j > 0 species or strains are said to be competing.

4. Secondary models

Based on previous data on the dependence of μ on temperature (Salvadó et al., 2011), two secondary models with a maximum of two parameters each were proposed: the Ratkowsky or square root model (SQRT) (Ratkowsky et al., 1983) (Equation 3) and the Arrhenius equation (Equation 4):

${\mu }_i\left(T\right)={\mu }_{i,0}^2{\left(T-{\mu }_{i,1}\right)}^2 \left(\mathrm{E}\mathrm{q}.\mathrm{ }2\right)$

${\mu }_i\left(T\right)= {\mu }_{i,0}*e^{\left(-\frac{{\mu }_{i,1}}{T}\right)} \left(\mathrm{E}\mathrm{q}.\mathrm{ }3\right)$

where, in Ratkowsky, µ_i,j,0 is the sensitivity of the response to temperature, and µ_i,j,1 is the minimum growing temperature. In Arrhenius, µ_i,j,0 is the maximum growing temperature while, µ_i,j,1 is the sensitivity. Note that the Arrhenius equations can also explain the case in which the biological parameter is constant with the temperature (with µ_i,j,1 = 0).

The other two factors, sequential co-inoculation time and SO₂ concentration, were included as constants. The first, timeCo, defines the presence or absence of the yeast S. cerevisiae. The second, β₁ and β₂, quantifies the sensitivity of each strain to molecular SO₂, acting as a secondary mortality rate.

The resulting system of equations is:

$\frac{{dx}_{sc}}{dt}=timeCo*x_{sc}*{\mu }_{sc}\left(T\right)*\left(1-{\left(\frac{x_{sc}}{K_{sc}}\right)}^{\theta sc}-{\alpha }_2*\frac{x_{wa}}{K_{sc}}\right)-{\beta }_1*{SO}_2 \left(\mathrm{E}\mathrm{q}.\mathrm{ }4\right)$

$\frac{{dx}_{wa}}{dt}=x_{wa}*{\mu }_{wa}\left(T\right)*\left(1-{\left(\frac{x_{wa}}{K_{wa}}\right)}^{\theta wa}-{\alpha }_1*\frac{x_{sc}}{K_{wa}}\right)-{\beta }_2*{SO}_2 \left(\mathrm{E}\mathrm{q}.\mathrm{ }5\right)$

5. Theoretical methods: the modelling process

Modelling was approached from a system identification perspective, which included the following steps: formulation of candidate models, parameter estimation from multiple experiments, and model reduction and selection.

6. Formulation of candidate models

The Gilpin–Ayala model was employed as the initial candidate owing to its higher level of complexity and considering the biology of the competition system (primary model) and different phenomenological submodels to account for temperature. The model solution depends on the given initial conditions, and the values of a series of unknown parameters.

7. Parameter estimation

The aim was to calculate the unknown parameters, constants related to growth, and kinetic parameters that minimise the distance between the data and the model predictions, in other words, the error.

One way to minimise the distance is through an objective function that includes the measure of distance using least squares:

$J\left(\theta \right)=\frac{1}{2}\sum _{i=1}^{n_{exp}}\sum _{j=1}^{n_{obs\left(i\right)}}{r_{ij}\left(\theta \right)}^T{Q}_{expmax,ij} r_{ij}\left(\theta \right)\mathrm{ }\left(\mathrm{E}\mathrm{q}.\mathrm{ }6\right)$

where n_exp is the number of experiments and n_obs is the number of observables. For each experiment and residuals are given by:

$r_{ij}\left(\theta \right)=y_{ij}-{\tilde{y}}_{ij}\left(\theta \right)\mathrm{ }\left(\mathrm{E}\mathrm{q}.\mathrm{ }7\right)$

here, $y_{ij}$∈ ${\mathbb{R}}^{n_{t,ij}}$ are the experimental data samples for observable j in experiment i, and ${\tilde{y}}_{ij}\left(\theta \right)$∈ ${\mathbb{R}}^{n_{t,ij}}$ are the corresponding model predictions. The parameter vector is denoted by $\theta$ ∈ ${\mathbb{R}}^{n_{\theta }}$.

Instead of using the identity matrix as weight, we adopted a normalisation strategy based on the maximum experimental value per observable and experiment:

$Q_{expmax, ij}= \frac{1}{{\left(y_{ij}^{max}\right)}^2} I , y_{ij}^{max}={max}_k{y}_{ij}\left(t_k\right) \left(\mathrm{E}\mathrm{q}.\mathrm{ }8\right)$

This choice scales the residuals by the maximum amplitude observed $y_{ij}^{max}$ in each time series $t_k$, ensuring a comparable contribution across observables with different units or magnitudes.

8. Model selection and reduction

The models were iteratively compared based on their ability to fit the experimental data. Since models with more parameters often have better fits but tend to overfit, the number of parameters was also considered. The Akaike Information Criterion (AIC) was used for this purpose, aiming to balance parsimony and the relative loss of information among candidate models by penalising the number of parameters (Burnham & Anderson, 2004). For least squares is defined as follows:

${AIC}_M=nd * \log \left(\frac{RSS}{nd}\right)+\mathrm{ }2\mathrm{ }np+\mathrm{ }1\mathrm{ }\left(\mathrm{E}\mathrm{q}.\mathrm{ }9\right)$

where, np is the number of parameters, nd is the number of data points, and RSS is the sum of squared residuals. The process started with the most complex candidate models, and after fitting the data, less influential parameters were iteratively eliminated and evaluated through the AIC strategy.

The minimum AIC_M value AIC_min was used to rescale AIC. The rescaled value ΔM = AIC_M – AIC_min was used to measure the model’s merit: models with ΔM ≤ 2 provide substantial support, models with 4 ≤ ΔM ≤ 7 provide considerably less support, and models with ΔM > 10 have no support (Burnham & Anderson, 2004).

9. Numeric tools

To automate the modelling process, we used the AMIGO2 toolbox (Balsa-Canto et al., 2020). AMIGO2 is a MATLAB-based tool focused on the identification and optimisation of parametric models, including sensitivity and identifiability analysis. It provides a set of numerical methods for both simulation and optimisation. From the available options, we selected CVODES (Hindmarsh et al., 2005) to solve the model equations and enhanced Scatter Search (eSS) (Egea et al., 2009) to find optimal parameter values within a reasonable time frame.

10. Initial inoculum simulations

The final model was used to develop simulations for the different co-inoculation times: 0, 24, and 48 h. The simulations (400) were performed using the Runge–Kutta fourth-order method, based on varying initial inoculums from 5 × 10⁶ to 1 × 10⁸ cells/mL of each yeast isolate.

11. Simulation validation

Model predictions were validated under experimental midpoint conditions (0.1 ppm molecular SO₂ and 17.5 °C). Validation was performed using three assays in which the initial inoculum concentrations were set to 5 × 10⁶ cells/mL of W. anomalus and 5 × 10⁵ cells/mL of S. cerevisiae, applied at sequential co-inoculation times of 0, 24, and 48 h.

Results and discussion

Mixed sequential inoculation is sought for various purposes: to enhance the organoleptic profile (Albertin et al., 2017), reduce ethanol concentration (Maturano et al., 2019), or biocontrol spoilage microorganisms (Canonico et al., 2023a; Canonico et al., 2023b). According to our literature review, this work represents an application of classical population dynamics modelling to the oenological context, specifically addressing co-inoculation processes for spoilage yeast biocontrol and their possibly influencing factors.

1. Obtained kinetics

Figure 1 presents the 15 kinetics of mixed sequential co-inoculation of Wickerhamomyces anomalus and Saccharomyces cerevisiae. These kinetics were conducted based on a Box–Behnken design, varying the initial conditions of co-inoculation times (0, 24, and 48 h), molecular SO₂ concentrations (0, 0.1, and 0.2 ppm), and temperatures (15, 17.5, and 20 °C).

The kinetics were developed up to day 6, as W. anomalus is expected to be present during this pre-fermentation period, contributing to biocontrol and even to the organoleptic profile, before significantly decreasing in population (Kuchen et al., 2024). Moreover, as a facultative aerobic yeast (Mehlomakulu et al., 2021), its persistence is strongly influenced by oxygen concentrations in the medium, which, under conditions of strong and predominant kinetics of S. cerevisiae, would force it to shift its metabolism to a less efficient state for this microorganism (Qin et al., 2024).

2. Model fitting

A total of n = 180 data points were used to fit the models. The Gilpin–Ayala model with submodels, described in the system of equations (3–6), was evaluated.

2.1. Parameter screening

The influence of temperature on the specific growth rate parameter was evaluated using submodels applied in previous studies (Salvadó et al., 2011; Balsa-Canto et al., 2020). Each best fit was obtained using a combination of Arrhenius or Ratkowsky submodels (Equations 3 and 4) (Ratkowsky et al., 1983), applied to describe the maximum specific growth rate of each species (Table 1). This evaluation was performed based on the goodness of fit for each combination of submodels with the Gilpin–Ayala model modified by us (Equations 5 and 6).

Table 1. Submodels used and goodness of fit.

Parameter	Submodels combination
µsc	Arrhenius	Ratkowsky	Arrhenius	Ratkowsky
µwa	Arrhenius	Arrhenius	Ratkowsky	Ratkowsky
Best fit	0.77	1.11	1.48	1.91

The model incorporating both Arrhenius-type submodels showed a better fit than the alternatives, consistent with its widespread use to describe the effect of temperature on biological parameters (Balsa-Canto et al., 2020). However, as ${\mu }_{i,1}$ = 0 the resulting exponential terms were equal to one, reducing the response “temperature” to a constant value and leading to the submodel removal. Biologically, this outcome is attributed to the narrow experimental temperature range (15–20 °C), initially considered representative of white wine fermentation (Petrignani et al., 2024). However, this choice limited the interpretation of the thermal effect on yeast interactions, constituting a relevant constraint of the study, especially considering red wine fermentation (Kuchen et al., 2021). From a mathematical standpoint, the Ratkowsky model describes the basal response of an organism during the phase preceding the increase in its growth rate, whereas the Arrhenius equation reflects its behaviour as it approaches the maximum growth rate. This suggests that the evaluated conditions are close to the optimal range for both yeasts, which can be considered low-temperature mesophiles, consistent with the findings of Petrignani et al. (2024) for the V22 strain.

Table 2. Parameter estimation, uncertainty, and screening summary.

Gilpin–Ayala models						Lotka–Volterra models
Parameter	Optimal value	Uncertainty	Parameter	Optimal value	Uncertainty	Parameter	Optimal value	Uncertainty	Parameter	Optimal value	Uncertainty
µsc,0	0.049	129 %	µsc	0.051	129 %	µsc	0.054	35.10 %	µsc	0.0453	21.40 %
µsc,1	0	Inf %	µwa	0.033	178 %	µwa	0.034	37.70 %	µwa	0.0302	26.40 %
µwa,0	0.028	111 %	α1	1.92	160 %	α1	1.55	125 %	α1	2.78	101 %
µwa,1	0	Inf %	α2	0.012	359 %	α2	0.0019	215 %	Ksc	7.64	2.01 %
α1	1.92	160 %	Ksc	7.66	6.69 %	Ksc	7.66	2.12 %	Kwa	7.81	4.77 %
α2	0.011	359 %	Kwa	7.68	13.10 %	Kwa	7.61	5.50 %
Ksc	7.66	6.69 %	θsc	1	176 %
Kwa	7.68	13.10 %	θwa	1	219 %
θsc	1.12	176 %	β1	0	Inf %
θwa	1.25	219 %	β2	0	Inf %
β1	1.03	120 %
β2	2.01	211 %

The reduced model, without temperature influence (submodels), was fitted using 10 parameters (Table 2). Intraspecific competition $\theta$, occurs within the same population, and can exhibit nonlinearity due to population density and resource consumption effects (Schoener, 1973). However, the estimated parameters converged to 1, indicating linear inhibition consistent with the original Lotka–Volterra formulation; therefore, the $\theta$ terms were removed. It may be possible to explain this outcome by the nutrient abundance in the must during the six-day fermentation period of this study, which likely prevented the occurrence of resource competition and, consequently, the detection of intraspecific inhibition effects (Fronhofer et al., 2024).

On the other hand, the parameters associated with SO₂, β₁, and β₂ converged to zero, indicating that it had no significant effect on the growth of either yeast. This could be attributed to the low concentration tested (0.2 ppm) and the known tolerance of W. anomalus (up to 0.4 ppm; Kuchen et al., 2019) and S. cerevisiae V22 to this compound (0.51 ppm; Petrignani et al., 2024). Although these limits were qualitatively known, it was assumed that a lower concentration would reveal changes in the kinetics, which was not observed. Therefore, this experimental choice represents another relevant limitation of the study.

Then, the Lotka–Volterra model was fitted considering only the co-inoculation time as a factor (R² = 0.774). Both interaction parameters α were positive, indicating competing organisms (Table 2). The value of α₂, representing the effect of W. anomalus on S. cerevisiae, was close to zero, indicating a minimal influence of the biocontrol yeast on S. cerevisiae. By contrast, α₁, on the order of unity, revealed that S. cerevisiae exerted strong competition over W. anomalus, indicating a predominantly unidirectional interaction.

These results show that the temporal sequence of inoculation significantly influenced the overall dynamics of W. anomalus. Such behaviour could be interpreted within the framework of the priority effect, defined as the phenomenon by which colonisation success and population dynamics, in this case W. anomalus, depend on the order of arrival to a local site or niche (Stroud et al., 2024). The inhibition of W. anomalus by S. cerevisiae could be attributed, at least in part, to competition for oxygen as a limiting resource. Since W. anomalus is a facultative aerobic yeast (Mehlomakulu et al., 2021), reduced oxygen availability would affect its metabolism and decrease its population. This effect could be further intensified by oxygen consumption by S. cerevisiae, limiting its availability to W. anomalus and reflecting potential competition for this resource (Hagman & Piškur, 2015).

For simplicity, robustness, and applicability, the interaction parameter of W. anomalus on S. cerevisiae α₂ was removed, yielding a fit of 0.781, which is slightly lower than the classical Lotka–Volterra model. The fitted parameters are shown in Table 2, and the model provides a good description of the data. Nevertheless, the comparison using the Akaike Information Criterion (AIC) (Table 3), performed to avoid overparameterisation, indicated low relevance of the physicochemical factors, except for co-inoculation time, and no improvement after removing α₂; therefore, it was retained in the model.

Table 3. Akaike information criterion for model evaluation.

Model	Gilpin–Ayala (timeCo, temperature, SO2)	Gilpin–Ayala (timeCo, SO2)	Lotka–Volterra	Lotka–Volterra unidirectional
	12	10	6	5
AICc	–230.3678	–242.4008	–251.9148	–252.4124
ΔAICc	12.033	10.0116	0.4976	0
	No support	No support	No considerably less support	Reduced model

2.2. Debugged model

$\frac{{dx}_{sc}}{dt}=timeCo * x_{sc} * {\mu }_{sc} * \left(1-\left(\frac{x_{sc}}{K_{sc}}\right)-{\alpha }_1*\frac{x_{wa}}{K_{sc}}\right) \left(\mathrm{E}\mathrm{q}.\mathrm{ }10\right)$

$\frac{{dx}_{wa}}{dt}=x_{wa} * {\mu }_{wa} * \left(1-\left(\frac{x_{wa}}{K_{wa}}\right)-{\alpha }_2*\frac{x_{sc}}{K_{wa}}\right) \left(\mathrm{E}\mathrm{q}.\mathrm{ }11\right)$

2.3. Simulated kinetics

2.4. Correlation analysis

There is a correlation between the interaction parameter of Saccharomyces cerevisiae on W. anomalus α₁ and the carrying capacity of the biocontroller yeast Kwa (Figure 3). This fact, could be biologically explained though the competition between both species likely increases cell mortality, which may reduce the maximum growth rate of W. anomalus. However, the release of amino acids associated with this cell death could contribute to an increase in carrying capacity by enhancing organic nitrogen availability (Jones, 1999; Toro & Vazquez, 2002). Furthermore, the commonly assumed positive relationship between equilibrium population density and competitive ability is often invalid (Fronhofer et al., 2024).

3. Time integral evaluation – Wickerhamomyces anomalus presence

The kinetic profiles showed better establishment and persistence of W. anomalus when co-inoculation was performed at 48 h, followed by 24 h. On average, co-inoculation at 48 h resulted in the highest growth and persistence of the biocontrol yeast; taking this condition as 100 %, co-inoculation at 24 h reduced its population approximately by 35.3 % and simultaneous inoculation (0 h) by 75 %.

Maintaining a high W. anomalus population over time is crucial for effective biocontrol against spoilage yeasts (Canonico et al., 2023a; Canonico et al., 2023b). However, this behaviour is likely not a constant behaviour, as the growth kinetics of W. anomalus tend to decline around the fourth day, probably due to the decrease in dissolved oxygen (~80 %) (Pérez-Magariño et al., 2023), which limits its aerobic metabolism (Mehlomakulu et al., 2021).

The impact of W. anomalus co-inoculation on S. cerevisiae was minimal, consistent with Gallo et al. (2024), who also found no relevant effects of Hanseniaspora vinae on the fermentative yeast. Nevertheless, Onetto et al. (2024) reported a significant interference of Hanseniaspora uvarum with S. cerevisiae, suggesting that the magnitude of the interaction depends on the specific non-Saccharomyces species involved and on the medium conditions.

4. Simulation and validation of different inoculum scenarios

Figure 4 illustrates the population dynamics across different initial inoculum concentrations and co-inoculation times: 0 h (Figure 4A), 24 h (Figure 4B), and 48 h (Figure 4C). This analysis is conceptually similar to the nullcline analysis (Zeeman, 1995). It evaluates the temporal evolution of both populations across 400 simulations, enabling visualisation of their overall trends.

The experimental validations (Figure 5), performed with high initial concentrations of W. anomalus and low concentrations of S. cerevisiae, confirmed the model simulations: in all three co-inoculation times, the fermenting yeast reached a carrying capacity close to 7.64 log₁₀ cells/mL (4.7 × 10⁷ cells/mL), converging to the same final sink state. This indicates that the order of arrival, or the population size of the colonising organism, did not negatively affect the competitive outcome, deviating from an inhibitory priority effect and instead suggesting a neutral or even facilitative effect (Stroud et al., 2024). Unlike natural systems described by Toju et al. (2018), where early-arriving species modified the niche and inhibited subsequent ones, here the anthropic introduction of S. cerevisiae imposed a predictable competitive direction.

Conclusion

The modified Lotka–Volterra model was successfully fitted and validated to describe the dynamics of sequential mixed co-inoculation between the non-conventional yeast and S. cerevisiae, considering different inoculation times and initial population sizes.

Within the evaluated experimental range, temperature and SO₂ concentration showed no significant effects on population dynamics, whereas co-inoculation timing was identified as the only determining factor.

Delaying the addition of S. cerevisiae transiently favoured the establishment of W. anomalus; however, the system rapidly shifted toward the competitive dominance of the fermentative yeast.

The order of arrival did not affect the final competitive outcome or sink state, contrary to expectations for an inhibitory priority effect. Even at low initial proportions, S. cerevisiae consistently prevailed in the experimental validations, reaching the same equilibrium point predicted by its carrying capacity (4.7 × 10⁷ cells/mL).

The greater persistence of W. anomalus and its lack of interference with S. cerevisiae under the conditions outlined in this study highlight its potential as a biocontrol agent, representing a promising strategy to reduce SO₂ use and promote more sustainable winemaking practices.

Acknowledgements

The authors wish to thank the national scientific community, particularly the National Scientific and Technical Research Council (CONICET), the Industrial Technological Institute (INTI), and the National University of San Juan (UNSJ), institutions to which we are proud to belong. We also thank ENAV S.A. for the grape must samples.

Declaration of conflicting interests

The authors declare that they have no known financial conflicts of interest or personal relationships that could have influenced the work presented in this article.