Identifiability, Sensitivity, and Genetic Algorithms in Bacterial Biofilm Selection Models

TL;DR

This paper develops a mechanistic model of bacterial biofilm formation, using genetic algorithms and sensitivity analysis to optimize data collection and understand long-term evolutionary responses to environmental stressors.

Contribution

It introduces a structured biofilm model incorporating evolutionary dynamics and applies genetic algorithms for optimal data collection strategies.

Findings

Simplified biofilm dynamics by removing detachment

Optimized data collection schedule using genetic algorithms

Identified key parameters influencing evolutionary speed

Abstract

Bacteria often develop distinct phenotypes to adapt to environmental stress. In particular, they can produce biofilms, dense communities of bacteria that live in a complex extracellular matrix. Bacterial biofilms provide a safe haven from environmental conditions by distributing metabolic workload and allowing them to perform complex multicellular processes. While previous studies have investigated how bacterial biofilms are regulated under laboratory conditions, they have not considered (1) the data requirements necessary to estimate model parameters and (2) how bacteria respond to recurring stressors in their natural habitats. To address (1), we adapted a mechanistic population model to explore the dynamics of biofilm formation in the presence of predator stress, using synthetic data. We used a Maximum Likelihood Estimation framework to measure crucial parameters underpinning the biofilm formation dynamics. We used genetic algorithms to propose an optimal data collection schedule that minimised parameter identifiability confidence interval widths. Our sensitivity analysis revealed that we could simplify the binding dynamics and eliminate biofilm detachment. To address (2), we proposed a structured version of our model to capture the long-term behaviour and evolutionary selection. In our extended model, the subpopulations feature different maximal rates of biofilm formation. We compared the selection under different predator types and amounts and identified key parameters that affected the speed of selection via sensitivity analysis.

Tables

Table 1. Table 1: Example parameter estimation values at different noise levels For given values of standard deviation measured in the Log-normal output distribution. Section   3.1 contains the details on the generation of synthetic data and a discussion of error magnitudes. Values are estimated using MLE, as described in Section   3.3 .

$\sigma$	${\chi }_{PB}^{\max }$	${\chi }_{PB}^{\min }$	$a$	${\chi }_{BP}$
Nominal	3.20	0.0320	0.0800	0.120
0.0025	3.20	0.0294	0.0803	0.117
0.0050	3.59	0.2630	0.0539	0.298
0.0100	3.19	0.0212	0.0811	0.108
0.0500	3.13	-0.0345	0.0869	0.055

Table 2. Table 2: Parameter confidence intervals for uniform and optimal sampling strategies and levels of data. The 95% confidence interval widths are given for each parameter estimation. The columns titled U-N apply a uniform sampling schedule consisting of N N sample points, including the initial condition. The Op-N columns have the genetic algorithm-optimised sampling schedules for N N sample points, including the initial condition. Where no value could be estimated on the interval, the value NaN is given. The search interval was extended (see Figure   2 B, where the grey region indicates the extended range) to quantify a value in the following best-case scenario where possible (e.g., 51 sampling times for χ P ​ B min \chi_{PB}^{\min} ). The values are also visualised in Figure   4 .

Parameter	U-11	U-21	U-51	Op-11	Op-21	Op-51
${\chi }_{PB}^{\max }$	0.159	0.073	0.028	0.0177	0.0123	0.0107
${\chi }_{PB}^{\min }$	NaN	NaN	NaN	NaN	NaN	NaN
$a$	0.480	0.188	0.071	0.0901	0.0708	0.0473
${\chi }_{BP}$	NaN	0.417	0.216	0.500	0.373	0.216

Table 3. Table 3: Various sensitivity indices quantifying impacts on B ∗ B^{*} for select phenotype transition parameters. The table shows several metrics to quantify the dependence of the final biofilm size, B ∗ B^{*} , on the model’s parameters.

Index type	${\chi }_{PB}^{\max }$	${\chi }_{PB}^{\min }$	$a$	${\chi }_{BP}$
Intercept	0.218	0.426	0.263	0.460
Gradient	0.218	0.010	0.172	-0.023
Pearson	0.749	0.034	0.592	-0.078
Partial	0.935	0.118	0.901	-0.266
Spearman CC	0.014	0.005	-0.002	-0.001
PRCC	0.014	0.010	0.001	0.000
Sobol tot	0.616	0.001	0.400	0.007
Sobol 1st	0.590	0.001	0.374	0.006

Table 4. Table 4: Population model parameters comparsion to original values Table of each of the parameters present in Equation ( 1 )-( 5 ), and Figure   1 A. The values are a nondimensionalization of those presented in [ 55 ] ; see Section   6.1 for details on this process. A brief description of each nominal non-dimensional value is given. The parameter values, or our dimensionalised value equivalents, from [ 55 ] , and their units are given in the final two columns. The first set of parameters is fixed and well-known within our model throughout. The second set ( a a , χ P ​ B max \chi_{PB}^{\max} , χ P ​ B min \chi_{PB}^{\min} , χ B ​ P \chi_{BP} ) are the unknown parameters, which we consider the variations of. Their nominal values used to generate synthetic data are given. The final set of five values then provides the initial condition values used for simulations, unless otherwise specified.

Variable	Name	ND Value	Value	Units
	Fixed parameters
$e_b$	Bacterial growth efficiency	0.2
$r_P$	Planktonic growth rate	25.2	0.21	h-1
$H_P$	Planktonic Half-Saturation	3.0	1	$\mu {\text{g ml}}^{-1}$
$r_B$	Biofilm growth rate	0.84	0.007	h-1
$H_B$	Biofilm Half-Saturation	8.0	1.0	$\mu {\text{g ml}}^{-1}$
$g_S$	Planktonic predator growth rate	5.76	0.12	$\mu {\text{g ml}}^{-1}$
$e_S$	Planktonic predator growth efficiency	0.5
$H_S$	Planktonic predator half-saturation	3.0	1.0	$\mu {\text{g ml}}^{-1}$
$g_A$	Biofilm predator growth rate	4.32	0.09	h-1
$e_A$	Biofilm predator growth efficiency	0.33
$H_A$	Biofilm predator half-saturation	0.8	0.1	$\mu {\text{g cm}}^{-2}$
	Re-estimated parameters
a	Attachment strength parameter	0.08	0.01	$\mu {\text{g cm}}^{-2}$
${\chi }_{PB}^{\max }$	Biofilm attachment rate maximum	3.2	0.05	cm h-1
${\chi }_{PB}^{\min }$	Biofilm attachment rate minimum	0.032	0.0005	cm h-1
${\chi }_{BP}$	Biofilm detachment rate	0.12	0.005	h-1
	Initial conditions
$C\left(0\right)$	Initial media concentration	1.0	0.33	$\mu$gC ml-1
$P\left(0\right)$	Initial plankton concentration	1.0	0.33	$\mu$gC ml-1
$B\left(0\right)$	Initial biofilm concentration	0.01	0.00125	$\mu$gC cm-2
$S\left(0\right)$	Initial ciliate concentration	0.03	0.01	$\mu$gC ml-1
$A\left(0\right)$	Initial amoeba concentration	0.008	0.001	$\mu$gC cm-2

Equations

$$\frac{dC}{dt}$$

$$\frac{dP}{dt}$$

$$\frac{dB}{dt}$$

$$\frac{dS}{dt}$$

$$\frac{dA}{dt}$$

$$\frac{dC}{dt}$$

$$\frac{d{P}_i}{dt}$$

$$\frac{d{B}_i}{dt}$$

$$\frac{dS}{dt}$$

$$\frac{dA}{dt}$$

$$\mathcal{P}$$

$$\mathcal{B}$$

$$P_i=\frac{A}{{\sigma }_{\chi }}\exp \left(-\frac{{\left(i-i_c\right)}^2}{2{\sigma }_{\chi }^2}\right),$$

$$Y_i(t)=X_i(t){\mathcal{E}}_i,\text{where}{\mathcal{E}}_i\sim \text{Log-Normal}(\mu ,\sigma ).$$

$${\mu }_l=\log \left(\frac{{\mu }^2}{{\mu }^2+{\sigma }^2}\right),$$

$${\sigma }_l^2=\log \left(1+\frac{{\sigma }^2}{{\mu }^2}\right).$$

$$\ell =\sum _{i,j,k}\log [p(Y_{i,j}(t_k)/X_i(t_k;\mathbf{\theta }))].$$

$${\mathbf{\theta }}_{opt}=\underset{\mathbf{\theta }}{\mathrm{argmax}}\sum _{i,j,k}\log [p(Y_{i,j}(t_k)/X_i(t_k;\mathbf{\theta }))].$$

$${\bar{\chi }}_{\max }=\frac{1}{\mathcal{B}}\sum _j{\chi }_j{B}_j.$$

Full text

Keywords: Parameter Identifiability, Optimal Experimental Design

1 Introduction

Many bacteria are cosmopolitan, but can survive under extreme conditions such as high or low temperature, humidity, salinity, and pH [65, 43, 31]. Bacteria have evolved to employ several key strategies to maintain the flexibility needed to address such various environmental challenges [9, 42, 60]. One such strategy is to adopt distinct ecotypes, observable differences between members of the same species in response to their environment. One type of adaptive response is the ability to form biofilms in response to both abiotic and biotic conditions.

Biofilms are communities of bacteria spread throughout an extracellular matrix, composed of secreted substances, particularly polysaccharides, proteins, and DNA [15, 14]. Biofilm formation can confer various forms of resilience to its occupants by creating spatially favourable distributions. This enables the bacteria to collectively share metabolic workload, and benefit from cooperative behaviours like quorum sensing [67, 37].

Free-living or planktonic bacteria move via diffusion or active swimming and readily adhere to surfaces or each other. These adhered bacteria become centres for more bacteria to attach, forming microcolonies that grow into mature biofilms through further attachment and internal replication. Researchers have examined the formation process using experiments [44], numerical simulations using hydrodynamic theory [68], ordinary and partial differential equation models[61, 10], and agent-based modelling [48]. These studies can provide valuable insights for several applications, including the prevention of biofouling [12, 38, 13], infections in medical devices [22], and applications to bioremediation [45]. By understanding biofilm formation dynamics under environmental stress, we also gain a crucial understanding of the wider ecological role of bacteria [8, 7, 21].

Differential equation models, where dynamics reflect expected in vivo behaviours, offer one way to model biofilm formation dynamics [4, 18]. In our study, we use these models to gain a deeper understanding of biofilm formation, particularly under external stressors. Our study addresses the types and amounts of data required to parameterise such models, as well as how existing forms of models can be adapted to capture long-term effects that are maintained through selection. Initially, we will focus on the model presented in [55], which utilises a set of differential equations to model a biofilm formation experiment under predation. The mathematical dynamics governing each population within such models follow established behaviours [47]. While these models are intuitively appealing due to their mechanistic interactions, they also make strong assumptions and introduce numerous potentially unknown parameters [20, 64]. As such, understanding how to measure these parameters and quantify their uncertainties is imperative [29, 58]. Therefore, we will determine a practical experimental setup that optimises the measurability of the parameters underlying our model. Additionally, this study enables us to explore how varying external stressors for such systems can lead to adaptations within the population through selection. The methodologies applied in this study allow us to develop a long-term model from a short-term model to study selection dynamics.

Our first aim is to evaluate this model using structural identifiability, practical identifiability, and sensitivity analysis, and use that information to propose an optimal data collection schedule [1]. Assuming that the model parameters can be categorised into a priori known and unknown categories, we determine the values and associated uncertainties of the unknown parameters. We validate the inverse problem on synthetic datasets generated from a known baseline parameter set, a common approach in such studies [6, 39]. The resulting re-estimated parameters and their confidence intervals are calculated from the synthetic dataset. Using relative practical identifiability, we use a genetic algorithm to propose optimal data sampling schedules. Further, we employ sensitivity analysis to determine which parameters drove significant differences in population dynamics. Finally, informed by the sensitivity and identifiability analysis results, we propose a simpler model that produces similar dynamics.

Our second aim is to model bacterial adaptation in a controlled setting. Researchers often use long-term evolutionary experiments (LTEE) to study these phenomena [33, 34]. LTEE involves the culture of numerous, initially identical, bacterial populations. Periodically subcultured populations in the LTEE allow the bacteria to proliferate continuously, resulting in many generations. By exposing subsets of these isolated cultures to fixed stressors, their development is guided by those selection pressures. The bacteria, in turn, will evolve different dynamics to offset the stressors. In a differential equation model, adaptation across generations can be thought of as bacteria changing the values of their fixed model parameters, or through generic functions that interpolate between conditions [35]. Assuming our model’s inverse problem is well-posed, experimental data could be used to perform the inverse problem at different generations of the LTEE, and then measure the various parameter values. However, how these parameter values change in this scenario is spontaneous, and there is no measurement of how they transition from their initial to final values. Since the initial model lacks a mechanism that enables these long-term changes to occur continuously, we will introduce a structured variant of this model, where the bacterial community is comprised of individuals with intrinsic differences in their dynamics. We will demonstrate using a systematic resetting of the population, mirroring an LTEE, that selection is then possible for such a structured model.

Section 2 introduced the two models used throughout the study. In Section 2.1 a short-term model is outlined, and in Section 2.2 we outline a modified long-term model. Section 3 introduced all the mathematical techniques applied throughout the results section. This section introduces structural and practical identifiability, the genetic algorithms techniques, and the sensitivity analysis. Section 4 discusses our results. We first outline the results of the short-term model’s structural and practical identifiability, before proposing a genetic algorithm approach to improve the model’s identifiability. Next, we introduce the long-term model in Section 4.4, and detail the sensitivity analysis of this model in Section 4.5. Finally, a discussion of the results is given in Section 5.

2 Population Models

2.1 Short-term Biofilm Formation Model

In the first part of our study, we focus on the dynamics of a population within one self-contained experimental setting. To represent the flow of nutrients through our biofilm-forming system, we use a stoichiometric mathematical model proposed in [55], see Figure 1A. The different populations are measured using state variables, which indicate the carbon content contained in these populations, access to which is limited by the initial content present [5, 25]. Our equations, which model the state variables, are given by

[TABLE]

Here, the state variables represent the compartments present in the system: the media source ($C$), planktonic bacteria ($P$), biofilm bacteria ($B$), planktonic predator ($S$), and biofilm predator ($A$). The carbon source (Equation 1) is consumed by both planktonic and biofilm cells. Equation 2 and Equation 3 represent planktonic and biofilm bacteria, respectively. Their similar structure includes (a) carbon-limited growth, (b) predation by their respective predators, (c) biofilm attachment, and (d) biofilm detachment. The attachment of planktonic bacteria to biofilm depends on the size of the biofilm: the bacteria attach at a maximum rate $\chi_{PB}^{\max}$ for small biofilm population sizes, but when biofilms are large, the rate reduces to $\chi_{PB}^{\min}$. The parameter $a$ determines the biofilm size at which the attachment rate transitions between these values. The detachment rate is proportional to the size of the biofilm, with a constant rate $\chi_{BP}$. Equation 4 and Equation 5 represent the planktonic and biofilm predators, which grow from consuming their respective food sources. In Equation 1-Equation 5, all growth and consumption functions follow a type-II functional response [19]. A schematic of nutrient flow is given in Figure 1A.

Table 4 presents the nominal values of the parameters used throughout our results. These values were calculated by non-dimensionalising the values from [55] (see Section 6.1). The non-dimensionalisation was chosen to provide consistency when comparing the relative sizes of the populations, regardless of whether they are in liquid (media, planktonic, or planktonic predator) or on the surface (biofilm, biofilm predator). However, our non-dimensional parameters implicitly depend on the geometry of the container used in experiments; thus, they must be re-determined for other experimental setups. Throughout, when referencing the values in Table 4, we will use a tilde (e.g., the nominal value of $\chi_{BP}$ will be denoted as $\tilde{\chi}_{BP}$). Our primary focus is to use this model to investigate changes in biofilm formation. Thus, we fix all parameters in the model except those most relevant to biofilm formation: $\chi_{\max},\chi_{\min},a,$ and $\chi_{BP}$. For a given synthetic dataset, our analysis will focus on re-estimating the values of these parameters and determining their uncertainties.

2.2 Long-term Structured Biofilm Model

The parameters in the model of Section 2.1 are fixed constants, which means the bacteria cannot adopt different dynamics or adapt. To study the long-term evolutionary dynamics of bacterial phenotypes, we propose a parameter-structure population model (see the schematic in Figure 6A). In this updated model, we assume that instead of having a single population of bacteria (and thus a single value of the constant parameter describing the binding dynamics), we have a collection of similar but slightly different subpopulations. Within these populations, the value of their maximum binding rate, referred to in the short-term model as $\chi_{PB}^{\max}$, takes on a continuum of values. The modified set of equations governing the system is:

[TABLE]

This set of equations is structurally similar to Equations (1)-(5), with the key difference being that we now have a compartment structure in $P$ and $B$. Equations (7) and (8) describe the dynamics of the bacteria occupying the $i$th compartment. The total amount of planktonic and biofilm populations is the sum of the individual compartments (Equation 11 and Equation 12). As before, there is growth from media source consumption and elimination due to predation. The growth and death terms are scaled to be proportional to the size of the compartment. In this updated model, there is no biofilm detachment rate. The remaining term concerns the biofilm attachment rate, which remains dependent on the size of the underlying biofilm. However, the minimum attachment rate is now zero for large biofilms. Each of the $N_{C}$ ordered compartments has an associated biofilm maximum attachment rate, $\chi_{i}$, at a uniformly spaced set of values (i.e., $\chi_{1}<\chi_{2}=\chi_{1}+\Delta\chi<...<\chi_{N}=\chi_{1}(N-1)\Delta\chi$). The bacteria are assumed to retain their maximum attachment rate $\chi_{i}$ upon transitioning between the phenotypes; thus, $P_{i}$ and $B_{i}$ are coupled, but do not interact with the other compartments. For clarity, we use the notation $\chi_{\max}$ in figures and discussions to represent the maximum value of the compartment. We select our range of $\chi_{i}$ so that evolving populations never reach the endpoints $\chi_{1}$ or $\chi_{N}$ [53]. To represent a continuum of $\chi_{\max}$ values, we used $N_{c}=2000$ compartments in our results. The dynamics for the carbon media source (Equation 6) and both predators (Equation 9 and Equation 10) are analogous to the previous model.

To solve Equations (6)-(12), the initial conditions are the same for $C$, $S$, and $A$, but the initial condition for $P_{i}$ or $B_{i}$ must be a distribution. To solve our system, we take the initial distribution:

[TABLE]

where $A$ is a constant chosen so that the sum across all compartments $\sum_{i}^{N_{c}}P_{i}=1$. The distribution has some width, measured by the standard deviation $\sigma_{\chi}$, and a mean $i_{c}$ corresponding to a binding rate of $\chi_{i_{c}}$. For convenience, we use $B_{i}=0.001P_{i}$ as the initial condition for the biofilm compartment to avoid issues related to calculations involving the total biofilm size.

To mimic the LTEE setup, the structured model is repeatedly solved to $t=1$, representing one generation, and re-seeded. This choice is motivated by an envisioned LTEE experimental setup in which the supernatant (containing media, planktonic predators, and biofilm predators) is removed from the system. Following this, the preserved biofilm is mechanically broken up to become the new planktonic population (e.g., initial conditions $P_{i}=B_{i}/\sum_{i}B_{i}$ and $B_{i}=0.001P_{i}$). New media and an appropriate level of predators are then reintroduced (e.g., $C$, $S$, and $A$ are reset to their initial values). By iterating this process, simulating the system to $t=1$ and then resetting the system using only the biofilm population, selection is enabled across multiple generations.

3 Methods

The dynamics of biofilm formation under predation are explored using the above population models. In this section, we outline the details of the techniques used to simulate these dynamics and describe the tools used to analyse model outputs.

3.1 Model Numerical Solutions and Synthetic Data Construction

The mathematical models of the populations (see Section 2.1 and Section 2.2) are solved numerically using the Julia package DifferentialEquations [50], specifically the Tsit5 implementation of the Runge-Kutta 4(5) [62].

We utilise noisy synthetic datasets to assess our ability to measure model parameters at various levels of sparsity. To generate our synthetic data for compartment $i$, $Y_{i}(t)$, we use the following statistical model:

[TABLE]

For each data point $X_{i}(t)$, we modify the exact numerical solution according to multiplicative noise drawn from a Log-normal distribution [46]. A synthetic dataset example is overlaid on the numerical solution in Figure 1B. We select the log-normal distribution to preserve the positivity of the synthetic data. The value of our Log-normal is given by $\mathcal{E}_{i}=\exp(\mu_{l}+\sigma_{l}\epsilon_{i})$ where $\mu_{l}$ and $\sigma_{l}$ are some real numbers and $\epsilon_{i}$ is a standard normally distributed random variable (i.e., mean zero and standard deviation one). These values $\mu_{l}$ and $\sigma_{l}$ tell us the mean and standard deviation of $\log(\mathcal{E}_{i})$, not of $\mathcal{E}_{i}$. We represent our distribution’s mean and standard deviation as $\mu$ and $\sigma$, respectively. The relationship between these values is as follows,

[TABLE]

Using the Julia packages Distributions and Random, the values $\mu_{l}$ and $\sigma_{l}$ were used to generate samples from the Log-normal distributions [3]. We assume that such multiplicative noise is appropriate to capture the errors present in measurements of each compartment, and that the values of $\mu$ and $\sigma$ are the same for each compartment. While we assume we know the value of $\sigma$ a priori, it can easily be estimated by extending our parameter estimation approach if it’s not known in advance [57].

In experiments, the measurement of bacterial population size through optical density or Colony-Forming Units is often a balance between speed and accuracy. Modern techniques report coefficients of variation in the range of 1-5% (corresponding to values as low as $\sigma=1\times 10^{-2}$) [51, 41]. In Section 4.1, we use an artificial value of $\sigma=2.5\times 10^{-3}$ to compare the relative difference between parameter estimation bounds, particularly in cases with very limited datasets. This noise represents an absolute best-case scenario value for these parameter estimate confidence bounds. At this noise level, we observe a range of practical identifiability across our parameters of interest, with some parameters estimated with very narrow 95% confidence intervals, while others have wide intervals.

3.2 Model Structural Identifiability

The first step in model evaluation is determining if the system is structurally identifiable. This intrinsic model property determines whether, given its noise-free outputs, the parameters within the model can be re-determined uniquely [64]. To determine whether the models presented are structurally identifiable, we used the online Maplecloud Structural Identifiability Toolbox [24]. This program uses the Structural Identifiability Analyser (SIAN) [23]. We determined whether model parameters were globally identifiable by providing different information availability. For example, in real-world experiments, the concentration of planktonic bacteria, the concentration in the biofilm, and the number of predators are relatively straightforward to measure (using Colony Forming Unit measurements (CFU), optical density (OD), haemocytometry, etc.) [49, 2, 66]. However, measuring the concentration of nutrients in the media, which can be done using techniques like liquid chromatography mass spectrometry, is comparatively impractical or destructive [16]. Restricting access to information about the media in a dataset (i.e., omitting synthetic data regarding $C$ throughout simulations) may pose issues regarding what parameters can be determined. Thus, structural identifiability analysis allows us to determine whether information about parameters within the system is even possible.

3.3 Model Practical Identifiability

Throughout our study, we will assess for given data scenarios whether the model is practically identifiable. This analysis allows us to determine the parameters within the model and quantify the statistical evidence that these measurements are significant.

We use the Maximum Likelihood Estimator (MLE) for a given set of synthetic data to re-estimate the baseline parameters used to generate the synthetic data [56]. To do this, suppose we have data $Y_{i,j}(t_{k})$ at a specific time $t_{k}$, where $i$ represents the state variable and $j$ represents the biological replicate. For a given initial guess at the baseline parameters $\bm{\theta}$, we can use our mathematical model to calculate what the parameter guess predicts the state variables should be, $X_{i}(t_{k}|\bm{\theta})$. The probability of the measurement $Y_{i}(t)$, $p(Y_{i,j}(t_{k})/X_{i}(t_{k};\bm{\theta})~|~\text{Log-normal}(\mu,\sigma))$, can be calculated since the distribution of the synthetic data is known. Working with the log-probability allows us to write the log-probability of the whole dataset as

[TABLE]

We have dropped the explicit error distribution label inside $p(\cdot)$ for simplicity. The best parameter set $\bm{\theta}_{opt}$ is the one which maximises the value of $\ell$ (equivalently, it is the parameter set which maximises the probability of measuring the data under our statistical model). We will refer to the optimal value as $\ell_{p}$. Hence,

[TABLE]

To calculate this, we used the NLopt Julia package [28]. Using an initial guess of the values of the parameters, the Nelder-Mead algorithm allows the optimal parameter values of Equation 18 to be found [17].

Without proper context, these estimated values can be somewhat misleading. The model is not practically identifiable with a given dataset if two distinct parameter sets provide the same log likelihood values. To quantify whether our parameters are practically identifiable, we calculated profile likelihood by fixing one parameter of interest at a time and then re-estimating the remaining parameters [52, 56]. For example, suppose we have measured the optimal values for our dataset $\bm{\theta}_{opt}=(\bar{\theta}_{1},\bar{\theta}_{2})$. We can, for some small value $\delta$, construct a new parameter set $\bm{\theta}^{\prime}=((1+\delta)\bar{\theta}_{1},\theta_{2})$. Where $(1+\delta)\bar{\theta}_{1}$ is fixed, and $\theta_{2}$ is re-estimated to maximise the log likelihood of our dataset. In that case, the log likelihood $\ell_{p}(\delta)$ must necessarily be less than the true optimal $\ell_{p}$ if $|\delta|>0$. By examining the value of $\bar{\ell}_{p}=\ell_{p}-\ell_{p}(\delta)$ across a range of values of $\delta$, we can determine how sensitive the estimation is to univariate variation. $\bar{\ell}_{p}$ is $\chi^{2}$-distributed [36] and, therefore, we can determine a one-degree-of-freedom quantile (for our univariate parameter) above which the probability of the data is equal to or greater than a desired level. This quantile can be implemented using the Distributions Julia package and has a value $\ell_{crit}=1.921$ for a 95% confidence level. A confidence interval can be constructed by interpolating $\bar{\ell}_{p}+\alpha$ and finding its roots. We refer to this interval as the confidence interval width $CIW_{95\%}$. The function CubicSplineInterpolation from the Interpolations Julia package, alongside the Roots package’s function find_roots, were used to perform these calculations [30, 63].

3.4 Genetic Algorithm

In our log-likelihood parameter estimation in Section 3.3, Equation 17, no pre-supposition is made on the sampling times $t_{k}$. One obvious schedule choice is uniform; however, different parameters can drive dynamics at different population sizes. For example, as in Figure 1B, when there is a tiny biofilm population, the attachment rate $\chi_{PB}^{\max}$ dominates the biofilm population size changes. In contrast, $\chi_{PB}^{\min}$ barely has any impact. To measure $\chi_{PB}^{\max}$ as accurately as possible, it is better to have more data samples when the parameter is more impactful. We used a genetic algorithm to systematically construct schedules that allow us to measure a given parameter most effectively [32]. A flowchart of the algorithm is given in Figure 3A.

As described in Section 3.3, for a given synthetic dataset, a 95% confidence interval can be calculated for the best fit parameters. The genetic algorithm we employed closely follows previous work [32], but is explained here for clarity. The genetic algorithm first generates $M=100$ random schedules to determine the best possible sampling schedule. These schedules consist of a fixed number of points; we generated 11, 21, and 51 sample point schedules in our case. The sample time points are generated by sampling a uniform distribution on $(0,1]$, with the initial condition always provided at t=0. Each new sample is initially checked to be greater than $\Delta t=0.01$ (14.4 minutes) away from other sample times. For each schedule, a synthetic dataset is generated, and the value of $CIW_{95\%}$ is found. The schedules can then be ranked according to the value of $CIW_{95\%}$, with the smallest $CIW_{95\%}$ values ranked highest. For the $M$ cases, we use a cloning process to generate a new set of schedules [32]. The best schedule is preserved unchanged. The new $i$th schedule in the cloned set is a copy of the $\text{ceil}(i/14)^{2}$th schedule in the previous set. This cloning method keeps many copies of the most highly ranked schedules and fewer, less highly ranked ones, up to the 51st best. Each cloned schedule time point is then modified by a value drawn from a normal distribution with a standard deviation $\sigma=9\Delta t/N_{epochs}$. The values were constrained to remain in the range $(\Delta t,1-\Delta t)$. Because of the relatively large standard deviation of the cloning noise, points could move freely within this interval, exchange positions, and the minimum distance between scheduled points was not enforced from one epoch to the next. The example from Figure 3B and Figure 3C used 11 schedule sample points (and three biological replicates) across five epochs. In later examples, three epochs were used to generate the optimal schedules, as in our example, relatively rapid settling of values motivated this approach.

3.5 Short-term Model Metrics and Sensitivity

We perform sensitivity analysis to assess the relationship between variations in biofilm dynamics parameters and their resulting population. Since we are interested in the dynamics of biofilm formation, we chose the final biofilm size as our quantity of interest, which we will refer to as $B^{*}$. We used Saltelli Sampling to generate a collection of 40,000 parameter sets, each of which the value of the parameters $\chi_{\max},~\chi_{\min},~a,~\chi_{BP}$ are varied over the range 10% to 190% of their nominal values (given in Table 4) [54]. In our results, we scale the parameters by their nominal value to allow for the best comparison between parameters.

We use a variety of metrics to analyse the relationship between the parameters and $B^{*}$. To display the general trend between parameters and $B^{*}$, we present univariate scatter plots (see Figure 5A-D) and their best-fit lines calculated via Julia’s CurveFit package. To measure linear relationships between the parameters and $B^{*}$, we calculate the Pearson and Partial correlations, the latter of which controls for the variation of other parameters. We also measure Spearman and Partial Rank Correlations, which allow us to remove the effects driven by outliers and measure the presence of monotonic relationships across the entire dataset. For a more in-depth discussion of these metrics, see [40]. The above metrics were calculated using the Julia package GlobalSensitivity [11].

Finally, we use Sobol Sensitivity indices to decompose the variance in the full dataset [59]. Sobol indices examine the magnitude of the present variance, even if the relationships between parameters are non-linear. Additionally, Sobol sensitivity provides a method for examining the collective effects of parameters, assessing potential relationships between them by measuring second-order indices.

3.6 Long-term Model Metrics and Sensitivity

In the long-term adaptation model, we apply sensitivity analysis methods to examine the relationship between model parameters and the change in biofilm binding rate. The average binding value can be found using

[TABLE]

Applying Equation (19) to the initial condition for the biofilm, we can get the initial biofilm average binding $\bar{\chi}_{\max}^{i}$ (which corresponds to $\chi_{i_{C}}$ from Section 2.2). We then recalculate the average binding rate after 50 generations, $\bar{\chi}_{\max}^{f}$. Our quantity of interest is the ratio of these rates $\bar{\chi}_{\max}^{f}/\bar{\chi}_{\max}^{i}$, which measures the relative impact of selection in the population. The parameters examined in the sensitivity analysis are the amount of resetting predators $S_{0}$ and $A_{0}$, the initial average binding rate $\chi_{i_{c}}$, and the variability of binding rates in the initial condition $\sigma_{\chi}$.

Sobol sensitivity was performed for a Saltelli sample using these methods (see Figure 7D), with 40,960 samples. The samples in this analysis were generated using the Python Library SALib sample.sobol’s function sample. The simulations were carried out in parallel using MATLAB to optimise the code runtime, the function ODE45 replacing the DifferentialEquations Tsit5 algorithm [26]. Since each sample’s execution took a relatively similar time, splitting the workload evenly between a MATLAB parallel pool of ‘threads’ workers heuristically resulted in the fastest execution. We saved the outputs for each parameter set, then analysed them separately using the Python Library SALib [27]. This library offers convenient functionality for importing parameters and datasets (note that results in Julia’s GlobalSensitivity must be generated locally).

4 Results

Our results fall into two main sections. In Sections 4.1 to 4.3, we will explore the short-term model presented in Section 2.1. We aim to determine whether the model parameters can be defined under specific data availability scenarios and assess the uncertainty in these measurements before examining the model’s sensitivity to variations in these parameters. In doing so, we show differences between the identifiability of the parameters and demonstrate that sampling schedule and frequency have a substantial effect on the parameters’ confidence intervals. Next, in Sections 4.4 to 4.5, we propose a new model for our system. This new model captures the impacts of selection on populations by introducing population structure within the bacterial compartments. Finally, we will measure this impact under specific conditions through a sensitivity analysis, showing that variability in bacterial parameters and predation type give the most significant changes in the selection dynamics.

4.1 Measuring the impacts of data sampling frequency on model practical identifiability

Our first aim is to explore the structural and practical identifiability of the short-term model described in Section 2.1. To determine the global structural identifiability of our initial mathematical model (Equation 1-Equation 5), we used the SIAN toolbox (see Section 3.2 for details). When information on the state variables $[C,P,B,S,A]$ is available, the complete set of parameters and initial conditions present in Table 4 is globally identifiable. In general, we will focus on the identifiability of the parameters $(\chi_{PB}^{\max},\chi_{PB}^{\min},a,\chi_{BP})$, since these parameters determine the dynamics through which the bacteria transition between planktonic and biofilm phenotypes. We assume all other parameters are known exactly a priori. When limiting our interest to this subset of parameters, only the biofilm $B$ needs to be measured for global identifiability.

To provide the best-case scenario for measuring these parameters practically, we will assume that we have access to measurements of all state variables. We generate a synthetic dataset with three biological replicates (described in Section 3.1) based on a uniform sample schedule as a proxy for experimental measurements. Additionally, the samplings include a measurement of the system’s initial condition. Figure 1B displays a synthetic data example using the nominal values from Table 4 for a single biological replicate with $\sigma=0.05$, comprising 21 equispaced time points (including the initial condition). Another example with all three biological replicates, the same noise level, and 51 uniformly spaced measurements is presented in Figure 1C.

We use MLE to estimate the parameters (see Section 3.3) from synthetic data. An example of the model with the best-fit parameters is displayed in Figure 1C. Qualitatively, we observe a good fit when visually comparing the curves in Figure 1C to the numerical solution in Figure 1B. Table 1 compares the recovered parameters for various noise levels up to $\sigma=0.05$. Despite visual agreement with the numerical solution, it is clear that the recovered parameter values sometimes differ from the nominal values. The estimation accuracy improves as the noise level decreases. For the remainder of this section, we will focus on the very-low noise case, $\sigma=2.5\times 10^{-3}$, since it consistently yields a very tight 95% confidence interval for parameters thought to be practically identifiable. This noise level also provides a significant, but measurable, confidence interval at some sampling levels for less identifiable parameters.

We next turn to profile likelihoods to explore parameter identifiability based on equispaced sampling schedules. Figure 2 displays the univariate profile likelihood (see Section 3.3 for details) to determine confidence intervals for each parameter under 11 (blue dotted), 21 (red dashes), and 51 (green line) equispaced time points. $\bar{\ell}_{p}$ is $\chi^{2}$ distributed, thus values greater than $\bar{\ell}_{p}$ (displayed as a horizontal yellow line) all have a 95% probability of being responsible for the synthetic dataset from which they are measured. Thus, the points where the blue, red, and green curves intersect this line represent the 95% confidence interval bounds of this estimation. We normalise using the nominal values on the $x$-axis in each plot to allow comparison of the relative widths between parameters. For all four parameters, the confidence intervals always contain the nominal value, and increasing the number of time points sampled consistently reduces the width of the confidence intervals, indicating improvement in measuring parameters. For the green curve (51 data points), the peak value and confidence intervals give a very tight interval around the nominal value for $\chi_{PB}^{max}$, Figure 2A. Parameters $a$ and $\chi_{BP}$, displayed in Figure 2C and Figure 2D, respectively, are relatively close ($\pm\sim 10\%$) to the nominal value. The range of values on the x-axis is extended for $\chi_{PB}^{min}$ (Figure 2B) to allow the confidence interval of the green curves to be measured (white and grey regions are given for visual reference to compare to the other plot’s ranges). The sparser datasets (11 and 21 measurements) in Figure 2A and Figure 2C produce reasonable estimated value measurements, with only Figure 2A giving a measurable interval in the range examined in 11 measurement datasets. In Figure 2B, the sparser datasets do not generate confidence intervals despite the larger domain, while the dense dataset can eventually give a value on a much larger range.

4.2 Applying Genetic Algorithms to produce optimal sampling schedules

Our previous results indicate that dense sampling is needed to identify parameters accurately when using uniform sampling. We now investigate whether alternative sampling schedules can improve parameter identifiability. We use a genetic algorithm to systematically evaluate numerous schedules and determine their corresponding confidence interval widths (see Section 3.4). A flowchart of the algorithm, along with examples, is presented in Figure 3. Keeping the same total number of data points as in Figure 2, a corresponding optimal schedule was produced, using three replicates in each synthetic dataset, based on its ability to minimise the corresponding parameter’s 95% confidence interval.

We compare the relative confidence interval widths between uniform and optimal sampling in Figure 4 and Table 2. For each parameter (separated by black dashed lines), there is a section representing uniform sampling (green background) and optimal sampling (grey background). Within each section, blue bars represent 11 time points, orange bars represent 21 time points, and green bars represent 51 time points. X’s represent a confidence interval larger than 0.7. As with uniform schedules, an increase in the available data (e.g., the number of time points) reduces the confidence interval, resulting in a more accurate parameter estimate. The trend in identifiability of the parameters is similar to uniform sampling, with $\chi_{PB}^{\max}$ and $a$ showing the smallest interval widths, followed by $\chi_{BP}$, and then $\chi_{PB}^{\min}$, which was not identifiable under any circumstances. By optimally choosing our sample schedules, we achieve a confidence interval width for 11 optimally spaced measurements that is comparable to, and in some cases narrower than, that of 51 uniformly spaced measurements. However, there are diminishing returns for optimal sampling compared to the improvement seen in the uniform sampling case. For example, when comparing $\chi_{PB}^{\max}$ measured with 11 vs. 51 sampling points, the confidence interval for uniform and optimal sampling was roughly one-sixth and one-half as wide, respectively. Even with optimal sampling, $\chi_{PB}^{\min}$ was highly non-identifiable, even in the highest data availability scenarios.

We observe stark differences in the identifiability of the parameters given these data scenarios. Overall, the maximal speed of binding $\chi_{PB}^{\max}$ is most readily measured and has the lowest associated uncertainty in all scenarios, suggesting that its value is essential in the overall dynamics of the system. The transition between fast $\chi_{PB}^{\max}$ and slow $\chi_{PB}^{\min}$ binding, governed by the value of $a$, as well as the binding off rate $\chi_{BP}$, were the following most readily measured parameters. This suggests they have some impact on the overall dynamics, but a lesser one than the maximal binding rate. Finally, we saw that the slow binding rate $\chi_{PB}^{\min}$ was not practically identifiable in all of the data scenarios. In contrast to $\chi_{PB}^{\max}$, the specific value didn’t readily impact the overall likelihoods, likely because the values were all small. As such, very few bacteria bind at this low rate, and so moderate changes to this tiny population are not easy to distinguish in light of the faster modifications dictated by the other parameters.

4.3 Characterising biofilm size sensitivity to formation parameters

Having examined the identifiability of the different parameters, we next turn towards sensitivity analysis to determine the impact of varying these parameters on the final size of the biofilm, denoted $B^{*}$ (see details in Section 3.5). Figure 5A-D displays the univariate dependence of $B^{*}$ on the normalised underlying parameter. Parameters are normalised to allow for a comparison of how proportional changes affect the output.

Table 3 depicts several sensitivity indices of the univariate relationships captured in Figure 5A-D. The strongest linear relationships, measured by either the gradient of the linear fit or the Pearson correlation coefficient, are present in $\chi_{PB}^{\max}$ and $a$. The remaining parameters only have weak linear dependence. The Partial Correlation for $\chi_{PB}^{\max}$ and $a$ reinforces our observation, providing a clearer relationship when controlling for the other parameters. Interestingly, the magnitude of the Partial Correlation for $\chi_{PB}^{\min}$ and $\chi_{BP}$ increases in both cases, suggesting the influence of these parameters is being suppressed by variation in other parameters. In all cases, the Spearman and Partial Rank correlations (PRCC) are small, indicating that outliers primarily drive the linear relationships suggested by our previous metrics. In Figure 5E, we can see this distribution of values, with the bulk of the values spread across $B^{*}=(0.2,0.8)$. The mean across the whole sample is $E(B^{*})=0.435$, which is not far from the nominal value $B^{*}_{0}=0.484$. Finally, to decompose the relative impact of the parameters on $B^{*}$ regardless of their underlying relationship, we calculated the Sobol indices via a variance-based decomposition (Figure 5F). Supporting what has been seen up to this point, $\chi_{PB}^{\max}$ and $a$ are dominant, accounting for 0.95 of the total variance. Additionally, we can see that their total-order Sobol index is very close to their first-order Sobol index. Thus, the influence of these two parameters has an independent impact on the quantity of interest.

Following our concluding remarks from Section 4.2, the results seen here further suggest that the overall dynamics of the bacteria are strongly impacted by the value of $\chi_{PB}^{\max}$. That is, with planktonic bacteria quickly able to bind to biofilms, the result is a significantly larger difference in final biofilm size $B^{*}$. The value of $a$ also has a relatively strong impact on this final biofilm size, likely since this value sets the duration for which the $\chi_{PB}^{\max}$ rate dictates the planktonic bacteria binding. As such, we see that bacteria that can maintain their fast-binding dynamics result in far larger final sizes $B^{*}$. The values of $\chi_{PB}^{\min}$ and $\chi_{BP}$ have relatively little impact on the final biofilm size, likely because their values are small. Therefore, bacteria that detach more readily or attach to large biofilms (i.e., $B>>a$) have relatively little impact on the overall final size of the biofilm.

4.4 Capturing bacterial evolution by selection with structured models

After evaluating the parameters that determine biofilm formation in the short term, we propose a structured model that includes a mechanism for long-term population changes through evolution (see Section 2.2). The key difference in our new model is that our system has multiple bacterial populations, each with distinct binding dynamics. An example schematic is presented in Figure 6A. This example features three compartments with different maximal binding rates, the values of which in the previous model were denoted $\chi_{PB}^{\max}$. Based on the results of the identifiability and sensitivity studies of the previous model, we removed the parameters $\chi_{BP}$ and $\chi_{PB}^{\min}$, since their impact on the system was minimal. Additionally, instead of a scalar initial condition for $P$ and $B$, the new populations $P_{i}$ and $B_{i}$ are given some initial distribution, namely normally distributed across some range of $\chi_{PB}^{\max}$ values, with the $i$th compartment having a value $\chi_{i}$.

As detailed in Section 2.2, we no longer focus only on simulations of single time intervals from $t=0$ to $t=1$, but instead take the state of the system at $t=1$ to reset the system periodically. In this way, it simulates multiple generations of the population subjected to a consistent and resetting environment. Figure 6B presents an example of the time interval $t=0$ to $t=1$ for these compartment populations. At $t=1$, our biofilm population is highlighted by a red curve; this population is used to reseed the planktonic population in the next generation.

We demonstrate the presence of selection by measuring the average binding rate $\bar{\chi}_{\max}$ present in our model both under no predation (Figure 6C) and with predation (Figure 6D). These figures display the probability distribution of the compartments $\chi_{\max}$ values. Without predators, the distribution shifts to the right as the system undergoes many generations and narrows. This right shift indicates that bacteria with a higher $\chi_{\max}$ value (i.e., those that more readily attach to the biofilm) are more represented at the end of each generation and are therefore being reseeded into the next generation. The increased peak height suggests that the distribution is narrowing over time, slowing the rate at which the curve moves to the right. Under no predation, we are naturally selecting for bacteria that can quickly attach to the biofilm, as these are the ones that make it to the next generation. The evolutionary effects of predator stress are displayed in Figure 6D. The main figure features biofilm predators, and the inset figure has planktonic predators, both of which are periodically refreshed. The biofilm predator drives opposite selection from the no-predator case, with the curve shifting to the left and favouring smaller values of $\chi_{\max}$ where bacteria transition to the biofilm more slowly. As before, this curve is narrowing, indicating that the process slows down as the bacteria undergo more generations. The planktonic predator dynamics look similar to the no-predator case, with bacteria with a higher $\chi_{\max}$ being favoured; however, it shifts right at a faster rate. The presence of predators results in very different types of selection within this model, which drives different outcomes in the bacterial population.

4.5 Characterising evolutionary dynamics’ sensitivity

To further understand the relationship between predation and bacterial populations, we can examine how the number of predators within the system affects their evolution. In Figure 7A and Figure 7B, we measure the average binding rate in the population, $\bar{\chi}_{\max}$, at the end of each generation. Figure 7A shows that over 500 generations, the population subjected to planktonic predators tends toward higher values of $\bar{\chi}_{\max}$, indicating that biofilms are forming quicker. These predators compound the selection effects in the no-predator system, further increasing the value of $\bar{\chi}_{\max}$ in the population. In contrast, Figure 7B shows a stronger relationship between the predator presence and the binding rate. Interestingly, Figure 7B transitions from an increasing $\bar{\chi}_{\max}$ value to a decreasing one as the predator quantity increases, as the biofilm predators overcome the selection from the experimental design. All the curves in Figure 7A and Figure 7B show varying slowing levels, with the gradient over time decreasing, suggesting each predator level may drive some asymptotic binding value on very-long timescales.

To understand the most significant drivers in selection for biofilm formation, we undertook a sensitivity study on the structured model (described in Section 3.6). We examined the effects of (1) the reset population of the planktonic predator $S_{0}$, (2) the biofilm predator $A_{0}$, (3) the initial population distribution mean $\bar{\chi}_{\max}$, and (4) the population initial standard deviation $\sigma_{\chi}$. The quantity of interest measured for our sampling was the ratio of the final population average binding rate to the initial. We used measurements at 50 generations based on the results presented in Figure 7A and B. These figures clearly distinguish their population’s average binding rates by 50 generations. Since these simulations are relatively computationally expensive, this choice balances execution times with simulations that are long enough to allow a clear distinction to result from the selection impacts. We examined a sample across the parameter space using the Sobol sampling method. Figure 7C displays the observed values distribution. We can see that the distribution is skewed to the left, with most of the samples reducing the final average binding rate compared to the initial. After measuring the overall variance in this distribution, we calculated the Sobol sensitivity indices to determine which parameters drive these variances (Figure 7D). Most of the variation is attributed to the variance of the parameters $A_{0}$ and $\sigma_{\chi}$, the resetting biofilm predator population size, and the initial bacterial binding variance, respectively. Unlike our previous sensitivity analysis, the first-order and total Sobol indices were unequal, indicating that the parameters are no longer independent. Instead, higher-order effects are also at play. The second-order indices have no significant effect, except the index between $A_{0}$ and $\sigma_{\chi}$, which had a value $0.274\pm 0.055$. Given first-order index values of $A_{0}$ and $\sigma_{\chi}$ of 0.381 and 0.317, respectively, the combination of these first-order values and the measured second-order between them explains $97.2\%$ of the overall variance in the data. This suggests that populations with solely high parameter variance, or solely under high levels of predation, are insufficient to drive the most significant levels of change. Instead, these two factors cooperate, driving the most considerable deviations from the nominal.

These results suggest that bacterial populations subject to such biofilm predators should be expected to display signs of selection in their biofilm formation dynamics most rapidly. While planktonic predation will drive some changes in biofilm formation dynamics, these changes occur more slowly. For less diverse populations of bacteria, as is common in lab strains, this lack of diversity likely limits the rates of selection since there is limited capacity of sub-populations to outperform others.

5 Conclusions

This study investigated models of bacteria switching between planktonic and biofilm-forming phenotypes, with results presented in two main areas. First, we explored the short-term population dynamics of bacteria under predation stress, assessing the real-world applicability of a mathematical model through structural and practical identifiability, optimal design, and sensitivity analyses. Recognising the short-term model’s limitation in reflecting long-term population changes, we then proposed a new structured long-term model. Our minimal approach divided phenotypic populations into compartments with varying behaviours, allowing some to outperform others. By periodically resetting our system and preserving the population distribution of this transition behaviour, mirroring the Long-Term Evolution Experiment (LTEE) design often conducted in laboratory settings, we enabled selection to occur within our model.

In the first part of our study, we analysed an existing stoichiometric model that describes bacteria transitioning between planktonic and biofilm phenotypes and their associated predators. We first confirmed the system’s global structural identifiability when information from all state variables was available. Furthermore, we demonstrated that restricting the parameters to those involved in biofilm dynamics ($\chi_{PB}^{\max}$, a, $\chi_{PB}^{\min}$, and $\chi_{BP}$), only information about the biofilm state variable ($B$), is required for global structural identifiability.

We then generated synthetic datasets from numerical solutions to estimate these biofilm dynamics parameters. We used Maximum Likelihood Estimation to estimate parameters from synthetic data. We found that, while the fitted solutions accurately reproduced the overall trends in the data, our practical identifiability study revealed that the accuracy of parameter reconstruction varied with sampling frequency and scheduling. In a best-case scenario, across a low-noise scenario with 11 to 51 measurements, the maximal binding rate of planktonic cells ($\chi_{PB}^{\max}$) and the biofilm size at which the binding rate decreased ($a$) exhibited the narrowest confidence intervals. This was likely due to initial conditions featuring a large planktonic population and a small biofilm, as is common in experimental setups. Conversely, $\chi_{PB}^{\min}$ and $\chi_{BP}$ were more challenging to measure accurately, likely because of their limited activity periods and smaller magnitudes. We propose that modifying experimental protocols, such as pre-growing biofilms, could enhance the measurement of these less identifiable parameters by shifting initial conditions to emphasise their roles in the dynamics. However, this approach might compromise the accuracy of other parameter measurements. Finally, our application of a genetic algorithm revealed that optimised sampling schedules could achieve data quality comparable to, or even exceeding, that of significantly higher uniform sampling rates, suggesting a limited inherent informativeness of state variables regarding underlying parameters.

Also, we conducted a sensitivity analysis on the initial model, which clearly demonstrated that $\chi_{PB}^{\max}$ and $a$ dominated the system’s overall variance. These results, combined with our parameter identifiability findings, led us to simplify the model for initially small biofilms. Thus, in our long-term structured model (Section (2.2)), we omitted the biofilm-to-planktonic transition rate ($\chi_{BP}$) and the large-biofilm binding rate ($\chi_{PB}^{\min}$). Our long-term model organises planktonic and biofilm populations into compartments, each defined by a unique maximum transition rate ($\chi_{PB}^{\max}$) from planktonic to biofilm. Initially, $\chi_{PB}^{\max}$ values are normally distributed within the planktonic bacteria. While the short-term evolution of this system qualitatively resembles the initial model, periodic reseeding using the biofilm’s $\chi_{PB}^{\max}$ distribution generates emergent differences over many generations. Selection was demonstrated to occur even in the absence of predators; this is likely driven by the reseeding mechanism, which acts as a selective pressure.

We observed that the presence of predators significantly alters the dynamics of selection. The biofilm-specific predator ($A$) had a significantly greater impact, resulting in faster adaptation. In contrast, a planktonic predator ($S$) had a much smaller effect compared to the inherent selection in the no-predator scenario. Interestingly, while the biofilm predator strongly influenced adaptation, the rate of adaptation slowed down more rapidly than in the other cases. This suggests a trade-off: organisms must transition quickly enough to form biofilm by the end of the generation, but not so quickly that they risk prolonged exposure and consumption within the biofilm. Consequently, the distribution of $\chi_{\max}$ shifts towards an optimal value that balances this trade-off.

We performed sensitivity analysis to understand how the ratio between final and initial mean binding rates is impacted by key model inputs: initial planktonic predator population size ($S_{0}$), initial biofilm predator population size ($A_{0}$), initial bacterial mean transition rate ($\bar{\chi}_{\max}$), and initial distribution mean ($\sigma_{\chi}$). This analysis revealed that variance is dominated by biofilm predator ($A_{0}$) and initial system variance ($\sigma_{\chi}$). The interaction between $A_{0}$ and $\sigma_{\chi}$ accounted for most of the unexplained variance. This is because $A_{0}$ directly impacts the biofilm population, exerting a strong and clear selective pressure. Similarly, a larger $\sigma_{\chi}$ enhances the system’s capacity to respond rapidly to selection, as it provides a greater initial diversity of competitively advantageous populations.

Our current models offer several avenues for future exploration. Although a comprehensive investigation of genetic algorithms was beyond the scope of this study, we found that optimal sample points often clustered. Future research could address this by either removing these clustered points or replacing them with new, randomly drawn points that maintain a minimum pairwise distance. Additionally, minimising the sum of all confidence intervals, as suggested by [32], could yield more robust sampling schedules. Finally, data sets incorporating multiple initial conditions could be employed to yield more accurate parameter set measurements, especially in models with regime-dependent identifiability.

Our structured model effectively captures selection dynamics; however, there is a need for further refinement to reflect real-world complexities, particularly the inherent limitations in bacterial cellular processes. Currently, our model simplifies these by characterising compartment differences solely based on $\chi_{\max}$ values, which could represent variations in attachment ability. However, selection likely also involves other parameters, such as bacterial reproduction rates ($r_{P}$ and $r_{B}$). Furthermore, real biological systems may impose discrete or bounded constraints on $\chi_{\max}$ values, and multiple underlying parameters often vary simultaneously under stress, potentially incurring energetic costs for traits like increased $\chi_{\max}$. Our model also assumed that all possible $\chi_{\max}$ values were initially present. To create more realistic models of systems, future studies should explore constrained and simultaneous parameter changes, incorporate mutation dynamics between compartments, and crucially, align these model additions with experimental observations to avoid unjustified theoretical assumptions.

Statements and Declarations

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

Data Availability: All data generated and codes used in analysis during this study are included in the published article’s associated GitHub at https://github.com/StephWilleniams/BiofilmModelAnalysis

Author Contributions

All authors contributed to the study conception and design. Stephen Williams performed numerics, data collection, and analysis. All authors read and approved the final manuscript.

Acknowledgements

This material is based upon work supported by the National Science Foundation under Grant No DBI-2214038.

6 Supplementary Materials

6.1 Nondimensionalisation of the Model

In the original model, as presented in [55], the carbon, planktonic, and planktonic predators have units of carbon per volume. The biofilm and its associated predator have units of carbon per area. These units are generally convenient, giving the closest description to what would be measured experimentally. However, the mixed-compartment units in the model can provide a somewhat obscured comparison of the relative sizes of the populations, making them challenging to compare. As such, the equations describing the populations were multiplied by their associated dimensions to avoid inconsistent comparisons, either a volume of $3\mu$l or an area of $8$cm2. With this change, the populations represent the total quantity of carbon stored in each state variable. For consistency, several parameters must also be handled this way, for example, $H_{C}$. Upon this dimensionalisation, the area/volume value in the original model equations was absorbed into other parameter values, where necessary. As such, the values of $\chi_{PB}^{\max}$ and $\chi_{PB}^{\min}$ have been scaled by the original system’s area/volume ratio, and $a$ by the system area. In addition, the parameters and timescales were also normalised so that 1 time unit corresponds to 1 day, since we entirely focus on simulated runs of this length. The efficiency calculation is performed to enforce model structure consistency, and the values of $r_{p}$ and $r_{b}$ are scaled by $1/e_{b}$. The efficiency value $e_{b}$ is applied in Equations (2) and (3), to be consistent with how it is used for Equations (4) and (5).

We choose our initial conditions for the system to be $C=1$ and $P=1$, corresponding to $0.3\dot{3}\mu$gC/ml in the dimensionalised units; the remaining populations are then chosen to be some orders of magnitude smaller. This choice provides a convenient way to consider how these populations evolve relative to their initial condition, which becomes an essential consideration in Section 2.2 to facilitate a simple reseeding process.

A limitation of this nondimensionalization approach is that the system’s parameter values are then tailored to the specific geometry, since the underlying parameters in [55] were scaled to fit the given experimental container. Another potential approach was scaling the state variables and parameters to the initial carbon content sum across compartments. This normalisation would allow a convenient way to view how the carbon is distributed proportionally. Additionally, that would clarify what carbon is lost through efficiency ($e_{b}$, $e_{S}$, $e_{A}$). However, the initial conditions favouring scaling were implemented, since we are primarily interested in between-parameter comparisons.

6.2 Code availability

The code used throughout is available on GitHub at:

https://github.com/StephWilleniams/BiofilmModelAnalysis

Where possible, fixed random seeding was used for all synthetic data in the figures to ensure reproducibility of all the data. Where practical, Julia was used to execute the calculations. The ModelAnalysis module was developed during the work and provides much of the functionality used for practical identifiability and genetic algorithm development. For convenience, some tasks required execution using MATLAB or Python. These tasks utilised the simple vectorisation and parallelisation capabilities that MATLAB allows.