# Modeling reveals the role of aging and glucose uptake impairment in L1A1 Listeria monocytogenes biofilm life cycle

Eva Balsa-Canto, Carlos Vilas, Alejandro López-Núñez, Maruxa

Mosquera-Fernández, Romain Briandet, Marta L. Cabo, Carlos Vazquez

► To cite this version:

Eva Balsa-Canto, Carlos Vilas, Alejandro López-Núñez, Maruxa Mosquera-Fernández, Romain Briandet, et al.. Modeling reveals the role of aging and glucose uptake impairment in L1A1 Listeria monocytogenes biofilm life cycle. Frontiers in Microbiology, 2017, 8, 10.3389/fmicb.2017.02118. [https://doi.org/10.3389/fmicb.2017.02118](https://doi.org/10.3389/fmicb.2017.02118)

HAL Id: hal-02619989

https://hal.inrae.fr/hal-02619989

Submitted on 25 May 2020

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
Biofilms are microbial communities that appear in biotic or abiotic surfaces (Costerton et al., 1987). *L. monocytogenes* can form biofilms on common food contact surfaces, such as plastic, polypropylene, rubber, stainless steel, and glass (Silva et al., 2008). In biofilms, *L. monocytogenes* finds protection from a variety of environmental factors, such as temperature, sugar, salt or pH (Moretro and Langsrud, 2004) and tolerates better biocides, hampering the process of decontaminating surfaces (Carpentier and Cerf, 2011). In fact, many *L. monocytogenes* strains have been isolated from food processing plants (Rodríguez-López et al., 2015) despite the programs to sanitize industrial facilities.

Resistance to stress, including resistance to biocides, is largely associated with biofilm structure (Costerton et al., 1987; Donlan and Costerton, 2002). Most biofilms exhibit complex structures in that patches of cell aggregates are scattered throughout an exopolysaccharide matrix, creating channels.

Plate counts have been widely used to analyze biofilms. However, plate counts may be misleading in the study of the progress of biofilms (Daims and Wagner, 2007), as they include only viable culturable cells and do not inform about the structure.

Imaging techniques allow for a more comprehensive study of biofilms. They allow mapping viable and damaged or dead cells (Tawakoli et al., 2013) or the distribution of extracellular polymeric substance but also reconstructing three-dimensional structures. Quantitative image analysis provides further insights by computing parameters which characterize structures of biofilms (Yang et al., 2000).

Recent works suggest several alternative work-flows and software tools for the systematic analysis of microscopy images. IMARIS (commercial software) enables the reconstruction of 3D structures. COMSTAT (Heydorn et al., 2000), ISA 3D (Beyenal et al., 2004), or PHLIP (Mueller et al., 2006) allow quantifying confocal laser scanning microscopy (CLSM) images. BIOFILMDIVER (Mosquera-Fernández et al., 2014) permits the quantification of 2D images for epifluorescence and CLSM. Vyas et al. (2016) used machine learning algorithms to analyse scanning electron microscopy images.

Modeling is a complementary tool for studying biofilm dynamics. The emphasis is paid into qualitative validations of developed models, which are able to recover experimentally observed structures. In that pursue, the new generation of biofilm models offer detailed descriptions of the formation of heterogeneous structures with clusters and mushrooms (see the reviews by Picioreanu et al., 2004; Wanner et al., 2006; Horn and Lackner, 2014). Continuous and hybrid models have been proposed that incorporate various mechanisms to describe biomass growth, spreading and detachment as well as nutrients transport and conversion.

Biofilms of *L. monocytogenes* show a variety of structures: mono-layers of adhered cells, flat unstructured multi-layers, honeycomb structures or clusters (Chae and Schraft, 2000; Djordjevic et al., 2002; Marsh et al., 2003; Rieu et al., 2008; Bridier et al., 2010; Pilchová et al., 2014). Recently, Guilbaud et al. (2015) reconstructed CLSM images to observe *L. monocytogenes* intra-species diversity in forming biofilms. The work considers the biofilm structures formed by 96 isolates and concludes that most strains form complex honeycomb-like structures at 48 h.

Previous works (Mosquera-Fernández et al., 2014, 2016) used quantitative image analysis throughout time to study the life cycle of biofilms formed by three *L. monocytogenes* strains. The analysis showed the presence of, at least, three phases: separate clusters which evolve to interconnected clusters, honeycomb-like or flat structures and a final detachment. The rates at which these phases occur vary significantly among strains.

In this study we focus our attention in the biofilms formed by L1A1 strain (serotype 1/2a-3a, lineage II). L1A1 biofilms exhibited a clearly distinctive dynamics: a fast transition toward a flat biofilm, which increases in thickness for a couple of days to massively detach after 4 days (Mosquera-Fernández et al., 2016). The mechanisms that explain this characteristic dynamics are not yet fully known.

The aim of this work is to provide some novel insights into the mechanisms that drive the dynamics of the biofilms formed by L1A1 *L. monocytogenes* strain. For this purpose, we combined quantitative image analysis, cell counts, nutrient uptake tests and modeling through a data-based model identification loop. We selected a deterministic continuous modeling framework. Deterministic reaction-diffusion models (RDM) offer the advantage of the reproducibility (Wanner et al., 2006). Besides, this type of models can be solved with advanced numerical techniques to guarantee the computational efficiency required for model identification through optimization based techniques (Balsa-Canto et al., 2010; Vilas et al., 2017). We formulated various candidate reaction-diffusion models to describe the system. The candidate models incorporate alternative mechanisms for growth, nutrient consumption, and detachment. Each model was reconciled with the measured data through optimization based data fitting. The best model was selected attending to a best compromise between the number of unknown parameters and its capability to quantitatively reproduce the measurements.

# 2. MATERIALS AND METHODS

## 2.1. Experimental Methods

### Bacterial Culture Conditions

Bacteria tested was *L. monocytogenes* L1A1 (serotype 1/2a/3a, lineage II) isolated from thermal gloves used in the fishing industry by the Microbiology and Marine Technology Products research group at IIM-CSIC.

Bacterial stock cultures were kept at -80°C in tryptone soy broth, TSB (BD Difco, USA), containing 50% glycerol in the ratio 1:1(v/v). For each experiment, L1A1 cells were grown in consecutive subcultures in TSB medium for 8 and 16 h at 37°C (optimal growth temperature). The subcultures were adjusted to an OD700 nm of 0.1(±0.001) which corresponds to approximately $10^8$ CFU/ml.

### Biofilm Formation

A volume of 50 µl of the overnight subculture was added to 200 µl of fresh TSB in each well of polystyrene of a 96-well microtiter plate (Greiner Bio-one, France). Microplates were incubated at 25°C under static conditions.
After 1 h of adhesion, the wells were rinsed to eliminate any non-adherent bacteria before being refilled with 200 µl of fresh medium. Biofilm harvesting was performed eight times per strain at 1, 4, 24, and every 24 h up to 120 h.

## Image Acquisition by Confocal Laser Scanning Microscopy (CLSM)

Images acquisition was performed using the high throughput method described by Bridier et al. (2010). The method relies on the use of microtiter plates with a µ-clear base allowing for a high-resolution imaging of biofilms. For visualization, we used the FilmeTracer LIVE/DEAD Biofilm Viability kit (Invitrogen, USA). The kit contains two different fluorescent nucleic acid markers, SYTO and propidium iodide (PI). SYTO fluorochrome penetrates in all cells while PI only penetrates in damaged cells. The fluorochromes generate a bi-color labeling in such a way that viable cells appear in green and non-viable cells in red.

A Leica SP2 AOBS confocal laser scanner microscope (Leica Microsystems, France) was used at the INRA MIMA2 Imaging platform. Scans were obtained at 400 Hz using 63 × 0.8 NA oil immersion objective with a 488 nm argon laser set at 25% of intensity.

One horizontal cross-section of biofilm, corresponding to a real biofilm area of 238 × 238 µm², defines one image or slice. Horizontal cross-sections were acquired consecutively along z-axis using a scanning step size of 1 µm. The collection of horizontal cross-sections corresponds to a stack. We collected seven stacks at diverse randomly selected locations in the biofilms. Two biofilm replicas per sampling time were considered (a total of 7 × 2 × 7 stacks). For each image acquisition, we recorded dual (green and red) emissions.

## Image Visualization

We used IMARIS software (www.bitplane.com/imaris/imaris) for the visual inspection of the images at different times and heights. IMARIS offers several view modes. The slice 2D mode allows the visualization of one horizontal cross-section at a particular height, i.e., at a given slice of the stack. The blend 3D view mode provides an aerial view of the observed region, thus helping to identify biofilm structures. The section view mode lets the user inspect several slices together. Moreover, the gallery view mode displays all slices as a continuous series of images allowing the user to analyze how the structure evolves along the z-axis. Green and red channels—GCh, green channel, corresponding to live cells and RCh, red channel, to non-viable cells—can be visualized separately allowing to easily detect the appearance of non-viable cells during the biofilm life cycle.

## Number of Adherent Cells

The number of adherent cells was determined according to Herrera et al. (2007). Samples were collected at 1 and 96 h from squares removed from the microtiter cavities and immersed in 10 ml of phosphate-buffered saline for 10 s to release non-adherent cells. Adherent cells were collected with peptone water-moistened swabs. After the squares had been rubbed twice with the swabs, they were transferred to 10 ml peptone water and subjected to 1 min of vortexing. The number of adherent cells was

determined by plating the appropriated serial dilutions on tryptic soy agar (Cultimed, Spain) after incubation at 37°C for 24 h.

## Analytical methods

At each sampling interval, the bulk was poured to an eppendorf and centrifuged at 9,000 g during 10 min. Supernatant was used to determine the carbon and nitrogen sources, sugar (Bernfeld, 1951) and protein (Lowry et al., 1951) contents for nutrient analysis. Sediment was washed twice with distilled water and dried to constant weight at 106°C for biomass estimation.

## 2.2. Theoretical Methods

### Image Quantitative Analysis

We used IMARIS software to compute the maximum thickness (MxT) of formed biofilms. MxT values were extracted for the eight replicas at each sampling time (up to 120 h) when cells were organized forming a clearly differentiable biofilm. The mean values obtained out of the eight replicas were used as the basis for model identification.

BIOFILMDIVER (Mosquera-Fernández et al., 2014) was used to obtain the area covered by cells (CA) for each slice. CA corresponds to the ratio between the number of colored pixels and the total number of pixels. In this way, $CA = 0\%$ is equivalent to the absence of cells while $CA \sim 100\%$ corresponds to a densely occupied slice. The CA value was computed for the biofilm including both channels (GCh and RCh) but also for both channels separately to analyse when and where non-viable cells appear in the biofilm.

### Modeling Approach

Model building was formulated as an iterative approach. To elucidate the mechanisms explaining the life cycle of the biofilms formed by L1A1 *L. monocytogenes* strain we proposed a set of candidate models M1–M4. Each model incorporated different mechanisms which in turn called for various unknown parameters. Unknown parameters were estimated using data fitting techniques within the AMIGO2 toolbox (Balsa-Canto et al., 2016).

For the sake of computational efficiency and reproducibility, all candidate models correspond to deterministic reaction-diffusion models in one dimension. The models consist of a set of (non-linear) partial differential equations (PDEs) which describe the spatio-temporal dynamics of biomass and nutrients.

### Model Simulation

The non-linear nature of the model candidates made the analytical approach impossible. Therefore, numerical techniques were employed to simulate the models. The domain of interest was discretized into a number *N* of smaller sub-domains. The use of a low order polynomial interpolation in each sub-domain allows to approximate the PDEs by a set of algebraic or ordinary differential equations (AEs or ODEs).

We considered two numerical approaches. The first combines the finite differences scheme in space - with centered differences for the nutrients and a backwards-forward space for the biomass- and the Crank-Nicolson approach in time (Balsa-Canto et al., 2017). A Newton-Raphson algorithm was used to
solve the resulting non-linear equations. The second combines the finite differences scheme as implemented in MATMOL (Vande Wouwer et al., 2014) (www.matmol.org) with ode15 s (implicit ODE solver included in MATLAB). Both approaches were handled in AMIGO2 (Balsa-Canto et al., 2016) using the black-box model option. Results obtained with both implementations coincide thus confirming the reliability of the numerical methods.

## Model Identification

Model identification consisted of two steps. First, we estimated unknown model parameters for each model candidate; second, we selected the final model attending to the goodness of fit.

Parameter estimation was formulated as a non-linear optimization problem to find the unknown model parameters which minimize the distance between the model predicted values and the available data. We solved the model under the specified experimental conditions for given parameter values and results were used to obtain maximum thickness (MxT) and average nutrient (AvgN) in the bulk liquid.

We computed MxT taking into account the location of the moving boundary, the interface between the biofilm and the bulk. We computed the average nutrient (AvgN) in bulk as mean[C_N(idx)] being idx the indexes that correspond to the bulk liquid.

AMIGO2 toolbox (Balsa-Canto et al., 2016) was used to solve the parameter estimation problem. In particular, we selected the least squares approach. The problem is formulated as: Find θ to minimize:

$$J_{ml} = \sum_{i}^{n_s} (mMxT_i - MxT_i(\theta))^2 + \sum_{i}^{n_s} (mAvgN_i - AvgN_i(\theta))^2 \quad (1)$$

subject to the dynamic constraints (the model) and bounds on the parameters θ_L ≤ θ ≤ θ_U.

In Equation (1), n_s regards the number of sampling times, mMxT and mAvgN the measured MxT and AvgN; θ corresponds to the vector of unknown model parameters.

The global optimizer eSS (Egea et al., 2010) was used to solve the parameter estimation problem due to its well-recognized efficiency and robustness.

Input files for simulation and parameter estimation are publicly available at AMIGO2 web page: https://sites.google.com/site/amigo2 toolbox/examples. Parameter ranges are included in Supplemental Data.

## Parametric Sensitivities

Local parametric sensitivities for a given observable o and at a sampling time t_s are defined as follows:

$$S_p^o(t_s) = \frac{\partial y^o}{\partial \theta_p}(t_s); \quad p = 1 \dots n_\theta \quad (2)$$

where y^o refers to either MxT or AvgN.

The corresponding relative sensitivities, S_p^o = Δθ_p / Δy^o, can be used to assess the individual local parameter influence or importance, that is to establish a ranking of parameters.

As an overall measure of the parametric sensitivity we will use:

$$\delta_p^o = \sum_{s=1}^{n_s} |s_p^o(t_s)| \quad (3)$$

# 3. RESULTS AND DISCUSSION

## Quantitative Image Analysis, Cell Counts, Nutrient Consumption Data and Modeling Can Be Combined to Investigate the Dynamics of Biofilm Structures

The scheme in Figure 1 is intended to reconcile model candidates with experimental data through parameter estimation and model selection. The approach combines experimental data from different sources: CLSM quantitative image analysis considering two colored channels, green for viable cells and red, for damaged or dead cells; cell counts; nutrients measurements through time and modeling to test different mechanisms that drive the life cycle of bacterial biofilms.

Recent works (Mosquera-Fernández et al., 2014, 2016) suggest that biofilm maximum thickness and covered area are the most informative parameters to quantify structure dynamics. In this work we selected IMARIS to reconstruct CLSM images, PHLIP to compute thickness throughout time and BIOFILDIVER to compute covered area as a function of time and z-axis. We also measured nutrients consumption in the selected sampling times as well as the number of adhered cells.

Alternative biological hypotheses generate alternative candidate models. The most appropriate model can be then selected, attending to its capability to explain the time-resolved data by model data fitting using the AMIGO2 toolbox (Balsa-Canto et al., 2016).

## L1A1 L. monocytogenes Forms Flat, Thick Structures and Experiments a Massive Detachment after 96 h

CLSM three-dimensional reconstructions (Figures 2a–f) show dense and homogeneous biofilms with scattered damaged or dead cells as stained by propidium iodide. Under the experimental setup used in this study, this strain forms rather flat, unstructured biofilms. After initial attachment, a thin biofilm is already present at 24 h. The flat structure is rather stable throughout time with a sustained thickness increase until 96 h. After 96 h the presence of damaged or dead cells is high.

The measured maximum thickness (mMxT) of the biofilms was computed for each of the experimental replicates over time. The distribution of maximum thickness values is shown in the Figure 2g in the form of box plots. The boxes indicate the degree of variability and skewness in the data. Results show how the median of the maximum thickness increases with time until 96 h. The median mMxT value at this time corresponds to 50 µm.

Rieu et al. (2008) showed the differences between the MxT value achieved in static and flow conditions for the AR0009 strain. A maximum value of 10 µm was achieved in static conditions after 48 h. Bridier et al. (2010) considered the maximum thickness at 24 h of the biofilms formed by 60
<figure><img src="image_2.png" /></figure>

**FIGURE 1** | Model identification scheme based on CLSM and nutrients consumption measurements. Biofilms were grown under static conditions. CLSM was used to gather image stacks in several sampling times. IMARIS allowed reconstructing 3D-structures and quantifying maximum biofilm thickness throughout time. BIOFILMDIVER enabled computing biofilm covered area as a function of time and z-axis. Nutrients consumed by cells were measured at each sampling time. We defined candidate models, estimated unknown parameters and selected the most appropriate model using data fitting in the AMIGO2 toolbox.

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

**FIGURE 2** | Dynamics of L1A1 *L. monocytogenes* biofilms during life cycle. (a–f) Present the three-dimensional reconstruction of the CLSM images captured at different times of the biofilms life cycle. (g) Presents the measured maximum thickness (mMxT) vs. time. Box plots represent the variability of maximum thickness over the different replicates. On each box, the central mark indicates the median, and the bottom and top edges of the box indicate the 25 and 75th percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers. The figure shows the median of MxT value increases with time up to 96 h. Maximum median value is of around 50 µm. After that MxT values decrease revealing a massive detachment episode between 96 and 120 h.
pathogens including 10 *L. monocytogenes* strains. Their results reveal values ranging from 16 to 35 µm for *L. monocytogenes* strains. Guilbaud et al. (2015) computed the mean thickness achieved by 96 different strains under static flow conditions. Their results reflect that most of the strains form biofilms with a mean thickness value in the range 15-25 µm at 48 h. The fact that L1A1 biofilms were around 36 µm at 48 h indicates that this strain forms quite thick biofilms as compared to other strains.

At 120 h the MxT decreases importantly indicating a massive detachment in that period of 24 h. The maximum relative variability among experimental replicates for a specific sampling time (t₅) can be computed as $\frac{[\max(MxT(t₅))]-\min(MxT(t₅))}{\min(MxT(t₅))}$. The maximum relative variability (19%) is found at 96 h. This may indicate that some replicates may experience a lighter detachment before 96 h. Massive detachment occurs for all replicates after 96 h.

## Covered Area Data Reveals that Live and Damaged or Dead Cells Are Arranged in a Multi-layer Structure

The covered area was calculated for the different slices, channels and sampling times. **Figures 3A,B** present the mean values

obtained over experimental replicates as functions of time and thickness (z-axis). The figures reflect that at each sampling time the biofilms are structured in layers with different CA values. **Figure 3A** presents the total covered area. The Figure shows how the maximum CA values are found in intermediate layers while CA is lower in the deepest layers and even lower toward the surface of the biofilm. Covered areas of <1% on the surface of the biofilm may reflect certain roughness of around a couple of µm or may correspond to weakly adhered cells. Lower covered areas in the deepest layers (CA ≤ 20%) might indicate the presence of voids or molecules that are not stained (including exopolysaccharides, proteins and eDNA). **Figure 3B** presents the results obtained for the red channel. The Figure shows the presence of damaged or dead cells from early times and in deep layers of the biofilm. Only at 120 h damaged or dead cells appear also in the surface of the biofilm.

**Figures 3C,D** present the evolution of the mean CA over all slices as a function of time for the different replicates. Box plots provide information on the variability among the different replicates. **Figure 3C** shows that mean CA is little at the earliest times. These low values indicate that at early times, cells tend to adhere on top of small clusters (of around 5 µm). Between 4 and 24 h there is a massive growth and deposition of cells which

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

**FIGURE 3** | Covered area vs biofilm thickness and time. **(A,B)** Present the mean CA values obtained over experimental replicates as functions of time and thickness (z-axis). The Figures show that biofilms are structured in layers of different CA values. Maximum CA is around 40%. At 120 h damaged or dead cells account for almost half of the area covered by the biofilm. **(C,D)** Present the mean CA values over all slices as a function of time for the different replicates. Box plots provide information on the variability among the different replicates. CA values are kept almost constant from 24 to 96 h. At 120 h CA increases substantially.
contribute to a substantial increase in biofilm dispersal. The CA reaches values in the range 15 – 20% which are kept until 96 h. The fact that covered area is maintained to a somehow constant value from 24 to 96 h may suggest that biofilm is only growing vertically. CA increases quite drastically from 96 to 120 h when a massive detachment occurs. This may suggest that some detached cells may be reallocated in the biofilm.

When analyzing the results obtained for the red channel (Figure 3D) we discover that a small amount of damaged or dead cells already appear at around 24 h in the deepest layers of the biofilm. The area occupied by damaged or dead cells is <6% for most of the time. Only at 120 h they appear distributed throughout the biofilm with CA values ranging between the 8% and the 12%. Remarkably at this time, damaged or dead cells account for almost half of the area covered by the biofilm, suggesting an episode of massive death which correlates with the large detachment described before.

Figures 2c,d show good reproducibility of results in most of the sampling times. A larger variability is observed in the red-channel values at 120 h. This may indicate that for the different replicates cells start dying massively at different times after 96 h thus contributing to a larger variability in the presence of damaged or dead cells.

Remarkably CA values are rather low throughout time whereas 3D-reconstructions (Figures 2b–f) show a large cellular distribution in the surface from 24 h. Damaged or dead cells also appear as early as 24 h. These results altogether may indicate the presence of small voids and macromolecules, including eDNA, supporting the hypothesis on the role of eDNA in L1A1 biofilms.

## Cells Prefer Glucose as Carbon Source But Stop Its Consumption Early in the Biofilm Life Cycle

Figure 4 presents the concentration of nutrients as measured in the bulk liquid. Data show that L1A1 consumes glucose as the primary carbon source.

Cells consume most of the glucose at early times. After 24 h glucose uptake stops even when the bulk still contains 0.4 mg/ml glucose.

## Glucose Impaired Uptake and Biofilm Aging Are Critical to Describe the Life Cycle of the Biofilms Formed by L1A1

Biofilms were grown in static flow conditions therefore our models are restricted to the hydrostatic case. Besides, since L1A1 biofilms are rather flat, we considered a one-dimensional domain, $\Omega = (0, L)$, with $L = 80$ µm, in which we distinguish two time-dependent regions. $\Omega_1(t)$, represents the bulk, and $\Omega_2(t)$, represents the part occupied by the biomass (See Figure 1 for the graphical representation). Both regions can be characterized by the value of the normalized biomass concentration, $C_B(t, x) = C_B/B_{max}$ and the normalized nutrients concentration, $C_N(t, x) = C_N/N_0$. $B_{max}$ corresponds to the maximum biomass concentration and $N_0$ corresponds to the initial nutrients concentration, both extracted from the experimental data.

We defined four candidate models (M1-M4), from the simplest reaction-diffusion model (RDM) with constant biomass diffusion and a linear biomass detachment, to the most sophisticated model, including glucose impaired uptake and a non-linear detachment due to biofilm aging. All candidate models can be embedded in the following mathematical formulation:

$$ \frac{\partial C_N}{\partial t} = \nabla \cdot (d_N(C_N)\nabla C_N) - K_1 G_N(C_N C_B) \quad (4) $$

$$ \frac{\partial C_B}{\partial t} = \nabla \cdot (d_B(C_B)\nabla C_B) + K_3 G_N(C_N C_B) - K_4 G_D(C_N, C_B) \quad (5) $$

Equation 4 describes the dynamics of nutrients. The first term in the right-hand side accounts for the diffusion of the nutrients and

<figure><img src="image_5.png" /><figcaption>FIGURE 4 | Nutrients in the bulk throughout time. (A,B) Present the absolute concentration of nutrients as measured in the bulk liquid. Boxplots reflect the variability among experimental replicates. (C) Shows the mean relative values using the concentrations at 1 h as the reference. Data show that L1A1 prefers glucose to protein as a source of nutrients. After 24 h cells stop consuming glucose.</figcaption></figure>
the second corresponds to nutrients uptake; $d_N$ is the diffusion coefficient which is different in the bulk and in the biofilm; $K_1 = \frac{\mu_B}{Y_{BN}} \frac{B_{max}}{N_0}$ is proportional to the ratio between the maximum specific growth rate of the microorganisms ($\mu_B$) and the nutrient growth yield $K_4 = Y_{BN}m_s$ and $K_4 = \frac{\mu_B m_s B_{max}}{Y_{bn}}$ which is proportional to the maintenance coefficient $m_s$.

Equation 5 describes the dynamics of the biofilm. In the right-hand side, the first term accounts for the diffusion of the biomass, the second for the production of biomass and the last one corresponds to the decay or detachment; $d_B$ is the diffusion coefficient and $K_3 = \mu_B/B_{max}$. We explored several possibilities for the definition of $G_N$ and $G_D$.

Equations 4 and 5 are similar to the ones used in the study of mass transfer and conversion in biofilms (Picioreanu et al., 2000; Eberl et al., 2001).

In our case bacteria adhere to the surface at $x = 0$, with a given thickness of 4.5 µm and there is neither flux of bacteria nor nutrients. Mathematically, boundary and initial conditions for biomass read as follows:

$$ \frac{\partial C_B}{\partial x}(t, 0) = 0, \quad t \in [0, T], \qquad (6) $$

$$ \frac{\partial C_B}{\partial x}(t, L) = 0, \quad t \in [0, T], \qquad (7) $$

$$ C_B(0, x) = \begin{cases} B_0, & \text{if } 0 \le x \le 4.5, \\ 0, & \text{if } 4.5 < x \le L. \end{cases} \qquad (8) $$

Moreover, a fixed nutrients concentration is fed at $t = 0$:

$$ \frac{\partial C_N}{\partial x}(t, 0) = 0, \quad t \in [0, T], \qquad (9) $$

$$ \frac{\partial C_N}{\partial x}(t, L) = 0, \quad t \in [0, T], \qquad (10) $$

$$ C_N(0, x) = 1, \quad x \in [0, L], \qquad (11) $$

All candidate models present the following features:

* There is a sharp front of biomass at the bulk/solid transition. This front can be used to determine the biofilm thickness.

* Biomass density can not exceed a maximum bound $B_{max}$ which is a parameter of the model restricted by the measured cell counts.

* Biomass production is due to nutrient consumption.

* Nutrients diffuse in the bulk and in the biofilm with different diffusion constants:

$$ d_N(C_N) = d_N, \quad x \in \Omega_1 \qquad (12) $$

$$ d_N(C_N) = d_{eff} d_N, \quad x \in \Omega_2, \qquad (13) $$

where $d_N$ corresponds to the constant glucose diffusivity in the bulk and $d_{eff}$ is the effective glucose diffusivity within the biofilm ($d_{eff} = 0.24$; Stewart, 1998).

We estimated the unknown parameters for the different candidate models using the time resolved data of biofilm maximum thickness and average nutrient concentration in the

bulk. Cell counts and cell dry weight were used to define initial conditions and bounds on the parameters. Details on the parameters to be estimated, bounds and problem formulation are included in the section Theoretical Methods (Table 1).

Figure 5 presents an overview of the candidate models and the corresponding best fit in terms of least squares error.

Model 1 (M1) assumes that biomass diffusivity is constant. Since there is an excess of nutrients in the medium, we assume that $G_N$ follows a mass action description $G_N = C_N C_B$ instead of the usual Monod formulation. Note that we also tested the Monod formulation without improvements in the sense of goodness of fit (data not shown). Detachment is assumed to be linear and dependent on $C_B$ as in, for example, Eberl et al. (2001).

Results in Figures 5A,B reveal that the standard reaction-diffusion models typically used to describe biofilms dynamics can not predict L1A1 biofilm life cycle. M1 predicts that cells consume all nutrients. Once the nutrients are exhausted the cells enter a phase of maintenance in which the thickness does not increase. Finally, the decay term dominates the growth term, and biomass depletion starts once cells maintenance is no longer possible. The maximum biofilm thickness is lower than the maximum found by quantitative image analysis, and decay starts earlier than expected.

<table><caption>TABLE 1 | Presents the parameters and the corresponding bounds considered for parameter estimation.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th><th>Model</th></tr></thead><tbody><tr><td>N<sub>0</sub> (mg/ml)</td><td>Initial glucose concentration</td><td>2.74 (Measured)</td><td>M1-M4</td></tr><tr><td>B<sub>max</sub> (mg/ml)</td><td>Maximum biomass</td><td>≥15.55 (Measured)</td><td>M1-M4</td></tr><tr><td>d<sub>N</sub> (m<sup>2</sup>/s)</td><td>Glucose diffusivity (bulk)</td><td>1 × 10<sup>-14</sup> – 6.6 × 10<sup>-10</sup></td><td>M1-M4</td></tr><tr><td>d<sub>eff</sub> (d)</td><td>Effective glucose diffusivity (biofilm)</td><td>0.24 (Stewart, 1998)</td><td>M1-M4</td></tr><tr><td>mu<sub>m</sub> (1/s)</td><td>Maximum growth rate</td><td>1 × 10<sup>-5</sup> – 5 × 10<sup>-2</sup></td><td>M1-M4</td></tr><tr><td>Y<sub>BN</sub> (d)</td><td>Biomass yield</td><td>1 × 10<sup>-3</sup> – 1</td><td>M1-M4</td></tr><tr><td>d<sub>B</sub> (m<sup>2</sup>/s)</td><td>Biomass diffusivity</td><td>1 × 10<sup>-16</sup> – 1 × 10<sup>-15</sup></td><td>M1,M3,M4</td></tr><tr><td>m<sub>s</sub> (1/s)</td><td>Maintenance coefficient</td><td>3 × 10<sup>-6</sup> – 5 × 10<sup>-5</sup></td><td>M1-M4</td></tr><tr><td>ε</td><td>Biomass diffusivity related constant</td><td>1 × 10<sup>-9</sup> – 1 × 10<sup>2</sup></td><td>M2</td></tr><tr><td>a (d)</td><td>Biomass diffusivity related constant</td><td>0 – 4 (Eberl et al., 2001)</td><td>M2</td></tr><tr><td>b (d)</td><td>Biomass diffusivity related constant</td><td>0 – 4 (Eberl et al., 2001)</td><td>M2</td></tr><tr><td>N<sub>min</sub> (d)</td><td>Threshold for glucose impaired uptake</td><td>0.128 – 0.153</td><td>M3-M4</td></tr><tr><td>k<sub>d</sub></td><td>Rate of activation of detachment</td><td>200 – 400</td><td>M4</td></tr><tr><td>D<sub>min</sub> (d)</td><td>% damaged or dead cells before detachment</td><td>0.05 – 0.1</td><td>M4</td></tr></tbody></table>

Parameter bounds considered for model identification. (d) Corresponds to dimensionless parameters.
<figure><img src="image_6.png" /><figcaption>**FIGURE 5** | Candidate models analysis. (A,B) Present the best fit to the data obtained for the different candidate models. (C,D) Show the spatio-temporal dynamics of the biomass and nutrients concentrations as predicted by the most successful model M4. (E) Shows the absolute value of the relative parametric sensitivities as computed for M4. (F) Presents the effect of modifying the parameter values on MxT.</figcaption></figure>

In order to better reflect maximum thickness and the time when it is achieved, we proposed the model 2 (M2). M2 includes the hypothesis that biomass spreading is significant once a certain density value is approached. To model this scenario we selected the expression proposed by Eberl et al. (2001):

$$d_B(C_B) = C_B^b \left( \frac{\epsilon}{B_{max} - C_B} \right)^a \quad (14)$$

where a and b are additional unknown parameters to be computed through data fitting.

The introduction of the non-linear diffusion of biomass induces a certain delay in nutrient consumption. However, the consumption becomes faster than in M1 after 20 h, thus causing thicker biofilms at early times. Both models converge to the same solution after around 60 h. The fact that the addition of two extra parameters (a and b) does not contribute to improving the quality of the fit implies that cell motility may be the driving force for cell spread at early times and not the cellular density.

Models M1 and M2 predict nutrient depletion in contrast to what it is experimentally observed. In Model 3 we hypothesize that L1A1 strain may present two glucose uptake systems with
distinct affinities. Parker and Hutkins (1997) described two glucose transport systems in *L. monocytogenes* Scott A. A low-affinity proton motive force (PMF) driven system and a high-affinity phosphotransferase mediated transport system (PTS). Both systems are active at high glucose levels (early times), whereas only the PTS system is active at low glucose levels (from 24 h). This fact may be related with the expression of the virulence regulator PrfA. On the one hand, PrfA is essential for biofilm formation in *L. monocytogenes* (Lemon et al., 2010; Zhou et al., 2011) and, on the other hand, Marr et al. (2006) showed that *L. monocytogenes* strains expressing high levels of the virulence regulator PrfA show a significantly reduced expression of PTS components of the glucose uptake. Moreover, carbon catabolite repression of at least two genes of the PfrA virulence regulon previously demonstrated in *L. monocytogenes* could explain the low level of biofilm formation observed at early times of culture (Milenbachs et al., 1997).

Model 3 (M3) is equal to M1 but includes the hypothesis of two glucose transport mechanisms, in such a way that up to a certain level of glucose ($N_{min}$), cells consume glucose following the mass action kinetics (as in M1) and once $N_{min}$ is achieved, glucose uptake is impaired.

M3 results in a better quality of fit (Figures 5A,B) since it can recover the glucose impaired uptake. However, decay is no longer recovered, and the predicted maximum thickness is lower than the one measured.

Model 4 is intended to describe detachment. The biological, chemical, and physical factors that drive detachment are complex and not fully understood yet. Chambless and Stewart (2007) reviewed the phenomena hypothesized as factors influencing detachment. The list includes fluid shear, the degradation of the EPS, absence of sufficient nutrients or quorum sensing, to name a few. From the point of view of modeling, fluid shear and substrate limited detachment have received notable attention (Xavier et al., 2005; Chambless et al., 2006). However, these mechanisms are not suitable to describe the system under consideration in the present work. Biofilms are grown in static conditions, so that fluid shear is assumed to be negligible and nutrients are not consumed.

We hypothesize that the large detachment observed after 96 h is related to cell death and the degradation of the extracellular DNA (eDNA), i.e., to biofilm aging. The presence of eDNA in biofilms appears to be associated with both lysis of cells and active secretion (Whitchurch et al., 2002). Following cell death, a sub-population of the dead bacteria lyse and releases DNA (Bayles, 2007). Harmsen et al. (2010) showed that *L. monocytogenes* eDNA might be the only central component of the biofilm matrix. Their work suggests that eDNA has a pronounced effect on the initial attachment of cells in static assays. Also the use of DNase I treatment resulted in an extensive removal of cell material at late stages of biofilm formation essays.

In Model 4 (M4) we introduced the hypothesis of the role of aging into the detachment term. As a measure of aging, we used the measured covered area of damaged or dead cells [$CBD(t)$] in such a way that the decay starts once a given value of damaged or

dead cells is present in the biofilm. Mathematically, this reads as follows:

$$G_D(C_B) = \frac{C_B}{1 + \exp(-k_d(CBD(t) - D_{min}))} \quad (15)$$

where $G_D(C_B)$ is a s-shaped function centered in $D_{min}$ and with slope $k_d$.

M4 is closer to the data (Figures 5A,B). Remarkably the least squares error is half the values achieved by M1-M3. The maximum biofilm thickness is in good agreement with the experimental data throughout time. The model predicts that the highest thickness is achieved at around 102 h. After that moment cells detach massively.

Figures 5C,D present the spatio-temporal dynamics of the biomass and nutrients in the domain as obtained for model M4. Results reveal that the biomass concentration is higher in layers which are deep in the biofilm, in particular between 10 and 40 h. The fast consumption of nutrients in the first hours (1–20 h) contributes to achieving the higher biomass concentration values at early times.

Parameter values corresponding to the optimal solution are reasonable. The diffusion coefficient for dilute glucose in water at 25°C is $6.6 \times 10^{-10}$; the glucose diffusion in the model corresponds to $d_N = 1.07 \times 10^{-11}$, which is reasonable taking into consideration that the bulk liquid contains other nutrients than glucose (e.g., proteins). The diffusivity in the biofilm is 24% that value (Stewart, 1998). To compensate for a rather low diffusivity of nutrients, the nutrient growth yield factor is rather large $Y_{BN} = 0.98$ which would indicate that the cells use very effectively the nutrients to grow and to produce EPS which would contribute to achieve high biomass concentrations at early times. The growth rate $\mu_B = 8.5 \times 10^{-3}$ is within the range of values published for *L. monocytogenes* species under different stress conditions (Augustin and Carlier, 2000). Note that this value does not correspond to the cellular growth but to the biofilm growth. The maintenance coefficient $Ms = 3.66 \times 10^{-5}$ would indicate that cells may survive up to around 7 h without access to nutrients. $N_{min} = 1.39 \times 10^{-1}$ corresponds to an average nutrient concentration of 0.38 which is a bit over the average nutrient value at the end of the process. $D_{min} = 5.21 \times 10^{-2}$ would indicate that a covered area of $CA = 5.21\%$ of damaged or dead cells would determine the beginning of the detachment which is rapidly linear as shown by the high value of $k_d = 311.5$.

To assess the influence of the different mechanisms into the biofilm dynamics, we computed the absolute relative parametric sensitivities at the optimum (see section Theoretical Methods). Results (Figure 5E) show that nutrients dynamics is critically affected by the glucose impairment parameter ($N_{min}$). The biofilms growth and the yield are also influencing the nutrients consumption. The biofilm thickness is affected by all parameters. The most relevant mechanisms are the biomass diffusion and growth as well as aging. Figure 5F shows how the dynamics of the maximum thickness varies when modifying the model parameters using a normal distribution with a 5% standard deviation. The model is quite robust to parameters'
modifications, the maximum thickness is consistently achieved at around 100 h presenting a low dispersion (around 10%). Furthermore, $D_{min}$ and $k_d$ have a considerable effect on the dynamics of detachment, both regarding intensity and velocity.

## 4. CONCLUSIONS

*Listeria monocytogenes* is a food-borne pathogen that can cause systemic infections in humans. In some cases food contamination has been linked to the capability of these bacteria to form biofilms, microbial structures with enhanced resistance to biocides. Therefore, it is of the highest interest to gain new insights into the life cycle of biofilms to design new elimination strategies.

These insights can be gained from a multidisciplinary approximation in which measurements and hypotheses are reconciled through a model identification procedure.

In this work, we proposed such a procedure combining quantitative CLSM image analysis, cell counts and nutrient consumption measurements with computationally efficient modeling and parameter estimation techniques.

The procedure allowed us to compare several candidate models, i.e., hypotheses, and to find the model, i.e., the mechanisms, that better describe *L1A1 L. monocytogenes* biofilms life cycle.

Remarkably glucose impaired uptake and biofilm aging are crucial to explain biofilms maximum thickness and the observed massive detachment in static cultures.

## REFERENCES

Augustin, J.-C., and Carlier, V. (2000). Mathematical modelling of the growth rate and lag time for *Listeria monocytogenes*. *Int. J. Food Microbiol.* **56**, 29–51.
doi: 10.1016/S0168-1605(00)00223-3

Balsa-Canto, E., Alonso, A., and Banga, J. (2010). An iterative identification procedure for dynamic modeling of biochemical networks. *BMC Syst. Biol.* **4**:11.
doi: 10.1186/1752-0509-4-11

Balsa-Canto, E., Henriques, D., Gabor, A., and Banga, J. (2016). AMIGO2, a toolbox for dynamic modeling, optimization and control in systems biology. *Bioinformatics* **32**, 3357–3359. doi: 10.1093/bioinformatics/btw411

Balsa-Canto, E., López-Núñez, A., and Vázquez, C. (2017). Numerical methods for a non-linear reaction-diffusion system modelling a batch culture of biofilm. *Appl. Math. Mod.* **41**, 164–179. doi: 10.1016/j.apm.2016.08.020

Bayles, K.-W. (2007). The biological role of death and lysis in biofilm development. *Nat. Rev. Micro.* **5**, 721–726. doi: 10.1038/nrmicro1743

Bernfeld, P. (1951). Enzymes of starch degradation and synthesis. *Adv. Enzymol. Relat. Subj. Biochem.* **12**, 379–428. doi: 10.1002/9780470122570.ch7

Beyenal, H., Lewandowski, Z., and Harkin, G. (2004). Quantifying biofilm structure: facts and fiction. *Biofouling* **20**, 1–23.
doi: 10.1080/0892701042000191628

Bridier, A., Dubois-Brissonnet, F., Boubetra, A., Thomas, V., and Briandet, R. (2010). The biofilm architecture of sixty opportunistic pathogens deciphered using a high throughput CLSM method. *J. Microbiol. Methods* **82**, 64–70.
doi: 10.1016/j.jimet.2010.04.006

Carpentier, B. and Cerf, P. (2011). Persistence of *Listeria monocytogenes* in food industry equipment and premises. *Int. J. Food Microbiol.* **145**, 1–8.
doi: 10.1016/j.ijfoodmicro.2011.01.005

Chae, M. S., and Schraft, H. (2000). Comparative evaluation of adhesion and biofilm formation of different *Listeria monocytogenes* strains. *Int. J. Food Microbiol.* **62**, 103–111. doi: 10.1016/S0168-1605(00)00406-2

## AUTHOR CONTRIBUTIONS

EB formulated the models; CVi and AL implemented the numerical simulations; CVi and CVá designed and supervised the numerical simulations; CVi and EB performed the optimizations; MM, MC, and RB designed and performed the experiments; EB designed the work and drafted the document; all authors read and approved the manuscript.

## FUNDING

EB and MC acknowledge funds received from the Axencia Galega de Investigación (GAIN). AL and CVá acknowledge the grants MINECO MTM2016-76497-R and Xunta de Galicia GCR2014/044, including FEDER funds. AL acknowledges the financial support from the Spanish MECD-FPU program (ref. FPU13/02191). MM acknowledges the financial support from the JAE-CSIC program co-founded by FEDER.

## ACKNOWLEDGMENTS

Julien Deschamps and Alexis Canette (INRA) are acknowledged for their help in biofilm cultivation and imaging.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2017.02118/full#supplementary-material

Chambless, J. D., Hunt, S. M., and Stewart, P. S. (2006). A three-dimensional computer model of four hypothetical mechanisms protecting biofilms from antimicrobials. *App. Environ. Microbiol.* **72**, 2005–2013.
doi: 10.1128/AEM.72.3.2005-2013.2006

Chambless, J. D., and Stewart, P. S. (2007). A three-dimensional computer model analysis of three hypothetical biofilm detachment mechanisms. *Biotech. Bioeng.* **97**, 1573–1584. doi: 10.1002/bit.21363

Cossart, P., and Lebreton, A. (2014). A trip in the “new microbiology” with the bacterial pathogen *Listeria monocytogenes*. *FEBS Lett.* **588**, 2437–2445.
doi: 10.1016/j.febslet.2014.05.051

Costerton, J. W., Cheng, K. J., Geesey, G. G., Ladd, T. I., Nickel, J. C., Dasgupta, M., et al. (1987). Bacterial biofilms in nature and disease. *Annu. Rev. Microbiol.* **41**, 435–464. doi: 10.1146/annurev.mi.41.100187.002251

Daimms, H. and Wagner, M. (2007). Quantification of uncultured microorganisms by fluorescence microscopy and digital image analysis. *Appl. Microbiol. Biotechnol.* **75**, 237–248. doi: 10.1007/s00253-007-0886-z

Djordjevic, D., Wiedmann, M., and McLandsborough, L. (2002). Microtiter plate assay for assessment of *Listeria monocytogenes* biofilm formation. *Appl. Environ. Microb.* **68**, 2950–2958. doi: 10.1128/AEM.68.6.2950-2958.2002

Donlan, R. M., and Costerton, J. W. (2002). Biofilms: survival mechanisms of clinically relevant microorganisms. *Clin. Microbiol. Rev.* **15**, 167–193.
doi: 10.1128/CMR.15.2.167-193.2002

Eberl, H. J., Parker, D. F., and van Loosdrecht, M. C. M. (2001). A new deterministic spatio temporal continuum model for biofilm development. *J. Theor. Med.* **3**, 161–175. doi: 10.1080/10273660108833072

Egea, J. A., Martí, R., and Banga, J. R. (2010). An evolutionary method for complex-process optimization. *Comp. Oper. Res.* **37**, 315–324.
doi: 10.1016/j.cor.2009.05.003

Guilbaud, M., Piveteau, P., Desvaux, M., Brisse, S., and Briandet, R. (2015). Exploring the diversity of *Listeria monocytogenes* biofilm architecture by high-throughput confocal laser scanning microscopy and the predominance
of the honeycomb-like morphotype. *App. Env. Microbiol.* 81, 1804–1810.
doi: 10.1128/AEM.03173-14

Harmsen, M., Lappann, M., Knochel, S., and Molin, S. (2010). Role of extracellular
dna during biofilm formation by *Listeria monocytogenes*. *Appl. Environ.
Microbiol.* 76, 2271–2279. doi: 10.1128/AEM.02361-09

Herrera, J. J. R., Cabo, M. L., Gonzalez, A., Pazos, I., and Pastoriza, L. (2007).
Adhesion and detachment kinetics of several strains of *staphylococcus aureus
subsp. aureus* under three different experimental conditions. *Food Microbiol.*
24, 585–591. doi: 10.1016/j.fm.2007.01.001

Heydorn, A., Nielsen, A. T., Hentzer, M., Sternberg, C., Givskov, M., Ersboll, B. K., et al. (2000). Quantification of biofilm structures by the novel computer program COMSTAT. *Microbiology* 146, 2395–2407.
doi: 10.1099/00221287-146-10-2395

Horn, H., and Lackner, S. (2014). Modeling of biofilm systems: a review. *Adv. Biochem. Eng. Biotechnol.* 146, 53–76. doi: 10.1007/10_2014_275

Kathariou, S. (2002). *Listeria monocytogenes* virulence and pathogenicity, a food safety perspective. *J. Food Protect.* 65, 1811–1829.
doi: 10.4315/0362-028X-65.11.1811

Lemon, K. P., Freitag, N. E., and Kolter, R. (2010). The virulence regulator PrfA promotes biofilm formation by *Listeria monocytogenes*. *J. Bacteriol.* 192, 3969–3976. doi: 10.1128/JB.00179-10

Lowry, O., Rosebrough, N., Farr, A., and Randall, R. (1951). Protein measurement
with the folin phenol reagent. *J. Biol. Chem.* 193, 265–275.

Marr, A. K., Joseph, B., Mertins, S., Ecke, R., Müller-Altrock, S., and Goebel, W. (2006). Overexpression of PrfA leads to growth inhibition of *Listeria monocytogenes* in glucose-containing culture media by interfering with glucose uptake. *J. Bacteriol.* 188, 3887–3901. doi: 10.1128/JB.01978-05

Marsh, E. J., Luo, H. L., and Wang, H. (2003). A three-tiered approach to differentiate *Listeria monocytogenes* biofilm-forming abilities. *FEMS Microbiol. Lett.* 228, 203–210. doi: 10.1016/S0378-1097(03)00752-3

Milenbachs, A. A., Brown, D. P., Moors, M., and Youngman, P. (1997). Carbon-source regulation of virulence gene expression in *Listeria monocytogenes*. *Mol. Microbiol.* 23, 1075–1085. doi: 10.1046/j.1365-2958.1997.2711634.x

Møro, T., and Langsrud, S. (2004). *Listeria monocytogenes*: biofilm formation and persistence in food-processing environments. *Biofilms* 1, 107–121.
doi: 10.1017/S1479050504001322

Mosquera-Fernández, M., Rodríguez-López, P., Cabo, M. L., and Balsa-Canto, E. (2014). Numerical spatio-temporal characterization of *Listeria monocytogenes* biofilms. *Int. J. Food Microbiol.* 182, 26–36.
doi: 10.1016/j.ijfoodmicro.2014.05.005

Mosquera-Fernández, M., Sanchez-Vizuete, P., Briandet, R., Cabo, M. L., and Balsa-Canto, E. (2016). Quantitative image analysis to characterize the dynamics of Listeria monocytogenes biofilms. *Int. J. Food Microbiol.* 236, 130–137. doi: 10.1016/j.ijfoodmicro.2016.07.015

Mueller, L. N., de Brouwer, J. F. C., Almeida, J. S., Stal, L. J., and Xavier, J. B. (2006). Analysis of a marine phototrophic biofilm by confocal laser scanning microscopy using the new image quantification software PHLIP. *BMC Ecol.* 6: 10.1186/1472-6785-6-1

Parker, C., and Hutkins, R. W. (1997). *Listeria monocytogenes* Scott A transports glucose by high-affinity and low-affinity glucose transport systems. *App. Environ. Microbiol.* 63, 543–546.

Picioreanu, C., van Loosdrecht, M. C. M., and Heijnen, J. J. (2000). A theoretical study on the effect of surface roughness on mass transport and transformation in biofilms. *Biotechnol. Bioeng.* 68, 355–369. doi: 10.1002/ (SICI)1097-0290(20000520)68:4<355::AID-BIT1>3.0.CO;2-A

Picioreanu, C., Xavier, J. B., and van Loosdrecht, M. C. M. (2004). Advances in mathematical modeling of biofilm structure. *Biofilms* 1, 1–12.
doi: 10.1017/S1479050505001572

Pilchová, T., Hernould, M., Prévost, H., Demnerová, K., Pazarová, J., and Tresse, O. (2014). Influence of food processing environments on structure initiation of static biofilm of *Listeria monocytogenes*. *Food Control.* 35, 366–372.
doi: 10.1016/j.foodcont.2013.07.021

Rieu, A., Briandet, R., Habimana, O., Garmyn, D., Guzzo, J., and Piveteau, P. (2008). *Listeria monocytogenes* EGD-e biofilms: no mushrooms but a network of knitted chains. *App. Environ. Microbiol.* 74, 4491–4497.
doi: 10.1128/AEM.00255-08

Rodríguez-López, P., Saá-Ibusquiza, P., Mosquera-Fernández, M., and López-Cabo, M. (2015). *Listeria monocytogenes*-carrying consortia in food industry. composition, subtyping and numerical characterisation of mono-species biofilm dynamics on stainless steel. *Int. J. Food Microbiol.* 206, 84–95.
doi: 10.1016/j.ijfoodmicro.2015.05.003

Silva, S., Teixeira, P., Oliveira, R., and Azeredo, J. (2008). Adhesion to and viability of *Listeria monocytogenes* on food contact surfaces. *J. Food Prot.* 71, 1379–1385.
doi: 10.4315/0362-028X-71.7.1379

Stewart, P. S. (1998). A review of experimental measurements of effective diffusive permeabilities and effective diffusion coefficients in biofilms. *Biotech. Bioeng.* 59, 261–272.
doi: 10.1002/(SICI)1097-0290(19980805)59:3<261::AID-BIT1>3.0.CO;2-9

Swaminathan, B., and Gerner-Smidt, P. (2007). The epidemiology of human listeriosis. *Microbes Infect.* 9, 1236–1243. doi: 10.1016/j.micinf.2007.05.011

Tawakoli, P. N., Al-Ahmad, A., Hoth-Hannig, W., Hannig, M., and Hannig, C. (2013). Comparison of different live/dead stainings for detection and quantification of adherent microorganisms in the initial oral biofilm. *Clin. Oral Investig.* 17, 841–850. doi: 10.1007/s00784-012-0792-3

Vande Wouwer, A., Saucez, P., and Vilas, C. (2014). *Simulation of ODE/PDE Models with MATLAB, OCTAVE and SCILAB: Scientific and Engineering Applications*. Dordrecht: Springer.

Vilas, C., Arias-Méndez, A., Garcia, M. R., Alonso, A. A., and Balsa-Canto, E. (2017). Towards predictive food process models: a protocol for parameter estimation. *Crit. Rev. Food Sci. Nut.* doi: 10.1080/10408398.2016.1186591.
[Epub ahead of print].

Vyas, N., Sammons, R. L., Addison, O., Dehghani, H., and Walmsley, A. D. (2016). A quantitative method to measure biofilm removal efficiency from complex biomaterial surfaces using SEM and image analysis. *Sci. Rep.* 6:32694.
doi: 10.1038/srep32694

Wanner, O., Eberl, H. J., Morgenroth, E., Noguera, D. R., Picioreanu, C., Rittmann, B. E., et al. (2006). *Mathematical Modelling of Biofilms*. Technical report, IWA Task Group on Biofilm Modelling.

Whitchurch, C. B., Tolker-Nielsen, T., Ragas, P. C., and Mattick, J. S. (2002).
Extracellular DNA required for bacterial biofilm formation. *Science* 295:1487.
doi: 10.1126/science.295.5559.1487

Wilks, S. A., Michels, H. T., and Keevil, C. W. (2006). Survival of *Listeria monocytogenes* scott a on metal surfaces: Implications for cross-contamination. *Int. J. Food Microbiol.* 111, 93–98. doi: 10.1016/j.ijfoodmicro.2006.04.037

Xavier, J. B., Picioreanu, C., and van Loosdrecht, M. C. M. (2005). A general description of detachment for multidimensional modelling of biofilms. *Biotechnol. Bioeng.* 91, 651–669. doi: 10.1002/bit.20544

Yang, X. M., Beyenal, H., Harkin, G., and Lewandowski, Z. (2000). Quantifying biofilm structure using image analysis. *J. Micro. Mthds.* 39, 109–119.
doi: 10.1016/S0167-7012(99)00097-4

Zhou, Q. C., Feng, F. F., Wang, L., Feng, X. Q., Yin, X. J., and Luo, Q. (2011). Virulence regulator PrfA is essential for biofilm formation in *Listeria monocytogenes* but not in *listeria innocua*. *Curr. Microbiol.* 63, 186–192.
doi: 10.1007/s00284-011-9964-7

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Balsa-Canto, Vilas, López-Núñez, Mosquera-Fernández, Briandet, Cabo and Vázquez. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.