Mathematical Modeling of Bacteria Communication in Continuous Cultures

Abstract

Quorum sensing is a bacterial cell-to-cell communication mechanism and is based on gene regulatory networks, which control and regulate the production of signaling molecules in the environment. In the past years, mathematical modeling of quorum sensing has provided an understanding of key components of such networks, including several feedback loops involved. This paper presents a simple system of delay differential equations (DDEs) for quorum sensing of Pseudomonas putida with one positive feedback plus one (delayed) negative feedback mechanism. Results are shown concerning fundamental properties of solutions, such as existence, uniqueness, and non-negativity; the last feature is crucial for mathematical models in biology and is often violated when working with DDEs. The qualitative behavior of solutions is investigated, especially the stationary states and their stability. It is shown that for a certain choice of parameter values, the system presents stability switches with respect to the delay. On the other hand, when the delay is set to zero, a Hopf bifurcation might occur with respect to one of the negative feedback parameters. Model parameters are fitted to experimental data, indicating that the delay system is sufficient to explain and predict the biological observations.

1. Background

More than twenty years ago it was first discovered that even primitive single-celled organisms such as bacteria are able to communicate with each other and coordinate their behavior [1,2]. Bacterial communication is based on the exchange of signaling molecules, or autoinducers, which are produced and released in the surrounding space. At the same time, bacteria are able to measure the autoinducer concentration in the environment, and according to this, they can coordinate and even switch their behavior, adapting to environmental changes. The term “quorum sensing” was coined to summarize the cell-to-cell communication mechanism thanks to which single bacteria cells are able to measure (“sense”) the whole population density [3]. Quorum sensing was first observed for the species Vibrio fischeri [2], which uses such a mechanism to regulate its bioluminescence. Nowadays, it is known that many bacterial species are able to use similar regulation systems, controlling biofilm formation, swarming motility, and the production of antibiotics or virulence factors [4,5,6].

The basis for cell-to-cell communication is a gene regulatory network that not only controls certain target genes, but often also their own production, resulting in a positive feedback loop. Gram-positive bacteria use so-called two-component systems (see e.g., [7]), whereas Gram-negative bacteria produce autoinducers directly in the cells, release them to and take them up from the extracellular space without any further modification or transformation.

In the following, we focus on the architecture of a quorum sensing system in Gram-negative bacteria, which mainly communicate via N-Acyl homoserine lactones (AHLs) [3,8], typically produced by a synthase. AHL molecules bind to receptors, which control the transcription of target genes. The receptor–AHL complex usually induces the expression of AHL synthases in a positive feedback loop.

We restrict our considerations to the bacteria species Pseudomonas putida, a root colonizing, plant growth-promoting organism [9]. Nevertheless, these basic principles may be easily transferred to related bacterial species.

Mathematical modeling of quorum sensing systems has developed in the last decade. Basic principles for a mathematical approach can be found, for example, in [10], where quasi-steady state assumptions for mRNA and corresponding protein in Pseudomonas aeruginosa were introduced, or in [11], which focuses on the basic feedback system of Vibrio fischeri and the resulting bistability. Alternative approaches for Gram-negative bacteria can be found in [12] (focussing on population dynamics) and [13] (including a further feedback loop). Classical mathematical models for Gram-positive bacteria were introduced, for example, for Staphylococcus aureus in [7,14,15].

Several model approaches have also been proposed for Pseudomonas putida, in closed systems (batch) as well as in continuous cultures (chemostat) [16,17,18]. The goal of this manuscript is to review such models, investigating mathematical properties and principles underlying the equations. The interesting component of quorum sensing models of Pseudomonas putida is that beside a positive feedback for the autoinducer one also finds a negative feedback via an autoinducer-degrading enzyme, a Lactonase. This is initialized with a certain time lag, leading to a system of delay differential equations (DDEs).

The paper is organized as follows. In Section 2 we provide a short overview of previous modeling approaches for quorum sensing of Pseudomonas putida. Starting from ordinary differential equations (ODEs) for the regulatory network in one single cell, in a second step we extend to quorum sensing in populations, including signal exchange among cells and Lactonase activity. The latter component introduces delays into the system. The delay represents the activation time of the Lactonase-dependent negative feedback. Bacteria population might be considered in batch as well as in continuous cultures. It is our purpose to investigate the long term behavior of the presented dynamical systems, and this can be achieved via a reduced model of two delay equations. We explain in great detail how to obtain the two-equation system, maintaining key properties of the gene regulatory network.

In Section 3 we present results concerning the existence and uniqueness of solutions to the reduced model. Moreover, we show that non-negative initial data yield non-negative solutions, a fundamental property of models in biology that is often violated when working with delay differential equations (cf. [19]). We compute stationary states of the dynamical system and investigate local stability properties. To this purpose, we compare the DDE system (with a constant delay, $\mathsf{\tau }>0$) to the associated ODE system ($\mathsf{\tau }=0$), studying delay-induced stability switches. In the last part of Section 3, model parameters are fitted to experimental data from [18], indicating that in the long run the reduced model is sufficient to explain and predict the general behavior of the system.

Everywhere in this manuscript, if not otherwise specified, we shall denote variables dependent on time by x or $x\left(t\right)$. First derivatives with respect to time are denoted by $\dot{x}$, respectively by $\dot{x}\left(t\right)$.

2. Methods

2.1. Compartmental Models

We present in the following compartmental models for quorum sensing of bacteria in a continuous culture. One compartment represents either bacterial population density, the nutrient concentration in the medium, or the concentration of a certain protein/enzyme/signaling substance in a single cell or in the medium.

2.1.1. Regulatory Pathway in One Cell

Let us start to consider the gene regulatory system for a single Pseudomonas putida cell. We follow a standard approach for modeling the quorum sensing system in Pseudomonas putida (ppu), analogous to the lux system in Vibrio fischeri [11], where polymers of the receptor–AHL complex initiate a positive feedback loop. The autoinducer concentration in Pseudomonas putida is regulated by a (self-induced) positive feedback as well as by a negative feedback via the AHL-degrading enzyme Lactonase. Transcriptional activators PpuR bind to AHLs, forming a PpuR-AHL complex which polymerizes. PpuR-AHL n-mers bind to the AHL-dependent quorum sensing locus (ppu-box) and synthesize PpuI. This protein is finally responsible for AHL synthesis. We neglect possible feedbacks (cf. [13,20]) on the transcription of PpuR, as these seem to be of minor influence [16]. Thus, just a constant basic production of the receptor PpuR is considered, as in [10,11,17]. Further, we assume as in [16,17,18] that PpuR-AHL n-mers induce synthesis of Lactonase molecules. A schematic representation of this regulatory pathway is given in Figure 1.

Everywhere in this work, mRNA equations are assumed to be in quasi-steady state. This assumption is justified by the evidence that many proteins are more stable than their own mRNA code (cf. [10] and references thereof). Let us denote the intracellular concentrations of PpuI and PpuR at time t by $I\left(t\right)$ and $R\left(t\right)$, respectively. The variables C and $C_i$ indicate the concentration of the PpuR–AHL complex and of its i-mer in one bacterial cell, respectively. It is assumed that the formation of complex i-mers takes place via the combination of an $(i-1)$-mer with a single PpuR-AHL complex, cf. [11]. AHL concentration will be denoted by x; in some cases, it might be convenient to distinguish between intracellular ($x_{int}$) and extracellular concentration ($x_{ext}$). The variable y denotes Lactonase concentration.

To begin with, we consider only the positive feedback which regulates AHL. The Lactonase-degrading activity shall be included in a separate step. The positive feedback loop of the regulatory pathway on the protein level described in Figure 1 can be written in the form of an ODE system (cf. [21]):

$$\begin{array}{cc}\hfill \dot{I}& =\underset{\begin{array}{c}\text{basic}\\ \text{production}\end{array}}{\underbrace{{\alpha }_I}}+\underset{\begin{array}{c}\text{feedback-regulated}\\ \text{production}\end{array}}{\underbrace{{\beta }_I\frac{C_n}{I_{th}+C_n}}}-\underset{\begin{array}{c}\text{natural}\\ \text{decay}\end{array}}{\underbrace{{\gamma }_{I}I}}\hfill \\ \hfill {\dot{x}}_{int}& =\underset{\text{production}}{\underbrace{\widehat{\alpha }I}}-\underset{\begin{array}{c}\text{natural}\\ \text{decay}\end{array}}{\underbrace{{\gamma }_A{x}_{int}}}-\underset{\begin{array}{c}\text{complex}\\ \text{formation}\end{array}}{\underbrace{{\pi }_1^{+}{x}_{int}R}}+\underset{\begin{array}{c}\text{complex}\\ \text{degradation}\end{array}}{\underbrace{{\pi }_1^{-}C}}+\underset{\begin{array}{c}\text{exchange}\\ \text{with}\text{medium}\end{array}}{\underbrace{d(x_{ext}-x_{int})}}\hfill \\ \hfill {\dot{x}}_{ext}& =\underset{\begin{array}{c}\text{exchange}\\ \text{with}\text{medium}\end{array}}{\underbrace{d(x_{int}-x_{ext})}}-\underset{\begin{array}{c}\text{natural}\\ \text{decay}\end{array}}{\underbrace{{\gamma }_A{x}_{ext}}}\hfill \\ \hfill \dot{R}& =\underset{\begin{array}{c}\text{basic}\\ \text{production}\end{array}}{\underbrace{{\alpha }_R}}-\underset{\begin{array}{c}\text{complex}\\ \text{formation}\end{array}}{\underbrace{{\pi }_1^{+}{x}_{int}R}}+\underset{\begin{array}{c}\text{complex}\\ \text{degradation}\end{array}}{\underbrace{{\pi }_1^{-}C}}-\underset{\begin{array}{c}\text{natural}\\ \text{decay}\end{array}}{\underbrace{{\gamma }_{R}R}}\hfill \\ \hfill \dot{C}& =\underset{\begin{array}{c}\text{complex}\\ \text{formation}\end{array}}{\underbrace{{\pi }_1^{+}{x}_{int}R}}-\underset{\begin{array}{c}\text{complex}\\ \text{degradation}\end{array}}{\underbrace{{\pi }_1^{-}C}}+\underset{\begin{array}{c}\text{dimer}\\ \text{degradation}\end{array}}{\underbrace{2{\pi }_2^{-}{C}_2}}-\underset{\begin{array}{c}\text{dimer}\\ \text{formation}\end{array}}{\underbrace{2{\pi }_2^{+}{C}^2}}+\underset{\begin{array}{c}\text{j-mer}\\ \text{degradation}\end{array}}{\underbrace{{\displaystyle \sum _{j=3}^n}{\pi }_j^{-}{C}_j}}-\underset{\begin{array}{c}\text{j-mer}\\ \text{formation}\end{array}}{\underbrace{{\displaystyle \sum _{j=3}^n}{\pi }_j^{+}C{C}_{j-1}}}\hfill \\ \hfill {\dot{C}}_i& =\underset{\begin{array}{c}\text{i-mer}\\ \text{formation}\end{array}}{\underbrace{{\pi }_i^{+}C{C}_{i-1}}}-\underset{\begin{array}{c}\text{i-mer}\\ \text{degradation}\end{array}}{\underbrace{{\pi }_i^{-}{C}_i}}+\underset{\begin{array}{c}\text{(i+1)-mer}\\ \text{degradation}\end{array}}{\underbrace{{\pi }_{i+1}^{-}{C}_{i+1}}}-\underset{\begin{array}{c}\text{(i+1)-mer}\\ \text{formation}\end{array}}{\underbrace{{\pi }_{i+1}^{+}C{C}_i}}\mathrm{for} 2\le i\le n-1\hfill \\ \hfill {\dot{C}}_n& =\underset{\begin{array}{c}\text{n-mer}\\ \text{formation}\end{array}}{\underbrace{{\pi }_n^{+}C{C}_{n-1}}}-\underset{\begin{array}{c}\text{n-mer}\\ \text{degradation}\end{array}}{\underbrace{{\pi }_n^{-}{C}_n}}.\hfill \end{array}$$

(1)

Although this regulatory pathway seems to be well understood, experimental settings cannot provide information on the dynamics of all components described in system (1). Typically only data for the time course of AHL (and for the population dynamics of the bacteria, which will be introduced in the next step) are available. For this reason, one is interested in a model reduction, decreasing the number of variables and parameters in the system of equations. In a first step, we assume the formation of complexes and its polymers to take place on a fast time scale. Quasi steady state assumptions ($\epsilon \to 0$) yield for the n-mer,

$$\epsilon {\dot{C}}_n={\pi }_n^{+}C{C}_{n-1}-{\pi }_n^{-}{C}_n\underset{\epsilon \to 0}{\underbrace{\to }}{C}_n=\left(\frac{{\pi }_n^{+}}{{\pi }_n^{-}}\right)C{C}_{n-1}$$

Consider now the $(n-1)$-mer for $\epsilon \to 0$ and substitute the last expression. We find

$$0={\pi }_{n-1}^{+}C{C}_{n-2}-{\pi }_{n-1}^{-}{C}_{n-1}+\underset{=0}{\underbrace{{\pi }_n^{-}{C}_n-{\pi }_n^{+}C{C}_n}}$$

It follows that

$$C_n=\left(\frac{{\pi }_n^{+}{\pi }_{n-1}^{+}}{{\pi }_n^{-}{\pi }_{n-1}^{-}}\right)C^2{C}_{n-2}$$

and, recursively,

$$C_n=\left(\prod _{j=2}^n\frac{{\pi }_j^{+}}{{\pi }_j^{-}}\right)C^n$$

We denote $p_I:={\prod }_{j=2}^n\frac{{\pi }_j^{+}}{{\pi }_j^{-}}$ and substitute the result of the quasi-steady state assumption into the I-equation in (1), obtaining

$$\dot{I}={\alpha }_I+{\beta }_I\frac{p_I{C}^n}{I_{th}+p_I{C}^n}-{\gamma }_{I}I.$$

Observe that the Hill coefficient n covers the fact that polymers (n-mers) of the complex PpuR-AHL are relevant for the positive feedback loop (see also [17]).

To reduce the system further, we also assume that PpuI is in quasi steady state, as in [11,13,17], for example, resulting in

$$I=\frac{{\alpha }_I}{{\gamma }_I}+\frac{{\beta }_I}{{\gamma }_I}\frac{C^n}{I_{th}/p_I+C^n}$$

Let $C_{th}:=\sqrt[n]{I_{th}/p_I}$, then the modified equation for $x_{int}$ reads

$${\dot{x}}_{int}={\alpha }_A+{\beta }_A\frac{C^n}{C_{th}^n+C^n}-{\gamma }_A{x}_{int}-{\pi }_1^{+}{x}_{int}R+{\pi }_1^{-}C+d(x_{ext}-x_{int}),$$

where ${\alpha }_A:=\widehat{\alpha }{\alpha }_I/{\gamma }_I$ and ${\beta }_A:=\widehat{\alpha }{\beta }_I/{\gamma }_I$

Diffusion through the cell membrane plays an important role in regulation processes. Nevertheless, AHL diffusion into and out of the cytoplasm does not require any transport mechanisms and the whole diffusion process goes rather fast, compared to the time scale chosen for the experimental measurements (1 h) [17,22]. This allows us to assume that $x_{int}$ and $x_{ext}$ are in equilibrium. Via steady state assumption, we get

$$x_{ext}=\frac{d{x}_{int}}{d+{\gamma }_A}\approx x_{int}$$

as $d\gg {\gamma }_A$. Taken together, the resulting AHL concentration (now simply denoted by x) follows

$$\dot{x}={\alpha }_A+{\beta }_A\frac{C^n}{C_{th}^n+C^n}-{\gamma }_{A}x-{\pi }_1^{+}xR+{\pi }_i^{-}C$$

and the simplified version of the single cell model (1) reads

$$\begin{array}{cc}\hfill \dot{x}& ={\alpha }_A+{\beta }_A\frac{C^n}{C_{th}^n+C^n}-{\gamma }_{A}x-{\pi }_1^{+}xR+{\pi }_i^{-}C\hfill \\ \hfill \dot{R}& ={\alpha }_R-{\pi }_1^{+}xR+{\pi }_1^{-}C-{\gamma }_{R}R\hfill \\ \hfill \dot{C}& ={\pi }_1^{+}xR-{\pi }_1^{-}C\hfill \end{array}$$

(2)

2.1.2. Population Dynamics

In the next step, the model is adapted for a bacterial population, including its growth in the classical experimental situation of a batch culture [17]. We denote the bacteria density in the medium at time t by $N\left(t\right)$. The dynamics of the bacterial population is classically described by logistic growth,

$$\dot{N}=rN\left(1-\frac{N}{K}\right)$$

where r is the bacterial growth rate and K the carrying capacity of the batch culture system.

Alternatively, one can consider the situation in a continuous culture, also called chemostat, with a continuous inflow of water and nutrient substrate for the bacteria and an outflow for all extracellular players. In this setting, one introduces a separate variable (S) for the available substrate concentration, which limits the bacterial growth. Consumption of nutrients is usually assumed to lead directly to a proportional increase of the biomass (N). The consumption term includes a saturation with the possibility of a further nonlinearity via the Hill coefficient $n_s$. Standard equations for nutrient–bacteria dynamics in a chemostat, with dilution rate $D>0$, are given by [23]

$$\begin{array}{cc}\hfill \dot{S}& =\underset{\mathrm{inflow}}{\underbrace{D{S}_0}}-\underset{\mathrm{consumption}}{\underbrace{{\gamma }_{S}N\frac{S^{n_s}}{K_m^{n_s}+S^{n_s}}}}-\underset{\mathrm{washout}}{\underbrace{DS}}\hfill \\ \hfill \dot{N}& =\underset{\mathrm{growth}}{\underbrace{aN\frac{S^{n_s}}{K_m^{n_s}+S^{n_s}}}}-\underset{\mathrm{washout}}{\underbrace{DN}}\hfill \end{array}$$

(3)

2.1.3. Lactonase Regulates AHL Degradation

It turned out by experimental observations [16,18] that a further process plays a major role in the AHL dynamics. In both the batch [16] and the continuous culture experiments [18], maximum concentrations of detected AHLs were followed by a rapid degradation of AHLs to Homoserines, indicating the presence of extracellular enzymatic activity. It is reasonable to assume that the AHL-degrading enzyme is a Lactonase [16], whose production or activation could also be initiated by polymers of the PpuR–AHL complex. Experiments in [16] suggested that Lactonases are activated with a certain delay (about 2 h) compared to the up-regulation of AHL production. From a mathematical point of view, this time lag can be included in the model via a delay differential equation [17].

2.1.4. Full Model

Let us now see how the regulatory pathway model (1), respectively the simplified system (2), can be adapted for a bacteria population. It can be convenient to distinguish between intracellular and extracellular components, and different assumptions are reasonable. For example, whereas in [17] the PpuR concentration was thought for the whole population, we consider here a system where the intracellular components (like PpuR) are interpreted per single (typical) cell.

In [18], to keep the model simple and at the same time to cover some details in the dynamics, equations for the concentrations of AHL (x) and Lactonase (y) in the medium, as well as one equation for the intracellular concentration of PpuR–AHL (C) were added to (3). At the same time, the total amount of PpuR (either free or in the PpuR–AHL complex) in one cell was assumed to be constant.

This does not correspond exactly to reality, but covers the idea that a cell typically maintains the number of receptors within a certain range. This simplification is justified by the still realistic resulting AHL-dynamics (see [18] for details).

The result is the following system of equations:

$$\begin{array}{cc}\hfill \dot{S}\left(t\right)& =D{S}_0-{\gamma }_{S}N\left(t\right)\frac{S{\left(t\right)}^{n_s}}{K_m^{n_s}+S{\left(t\right)}^{n_s}}-DS\left(t\right)\hfill \\ \hfill \dot{N}\left(t\right)& =aN\left(t\right)\frac{S{\left(t\right)}^{n_s}}{K_m^{n_s}+S{\left(t\right)}^{n_s}}-DN\left(t\right)\hfill \\ \hfill \dot{x}\left(t\right)& =\underset{\begin{array}{c}\text{total}\\ \text{AHL}\text{production}\end{array}}{\underbrace{\left({\alpha }_A+{\beta }_A\frac{C{\left(t\right)}^n}{C_{th}^n+C{\left(t\right)}^n}\right)N\left(t\right)}}-\underset{\begin{array}{c}\text{natural}\\ \text{decay}\end{array}}{\underbrace{{\gamma }_{A}x\left(t\right)}}-\underset{\begin{array}{c}\text{complex}\\ \text{formation}\end{array}}{\underbrace{{\pi }_1^{+}(R_{const}-C\left(t\right))x\left(t\right)}}\hfill \\ & +\underset{\begin{array}{c}\text{complex}\\ \text{degradation}\end{array}}{\underbrace{{\pi }_1^{-}C\left(t\right)}}-\underset{\begin{array}{c}\text{washout}\end{array}}{\underbrace{Dx\left(t\right)}}-\underset{\begin{array}{c}\text{Lactonase-regulated}\\ \text{degradation}\end{array}}{\underbrace{\delta x\left(t\right)y\left(t\right)}}\hfill \\ \hfill \dot{C}\left(t\right)& ={\pi }_1^{+}(R_{const}-C\left(t\right))x\left(t\right)-{\pi }_1^{-}C\left(t\right)\hfill \\ \hfill \dot{y}\left(t\right)& =\underset{\begin{array}{c}\text{total}\\ \text{Lactonase}\text{production}\end{array}}{\underbrace{{\alpha }_L\frac{C{(t-\mathsf{\tau })}^m}{C_{th2}^m+C{(t-\mathsf{\tau })}^m}N\left(t\right)}}-\underset{\begin{array}{c}\text{natural}\\ \text{decay}\end{array}}{\underbrace{{\gamma }_{L}y\left(t\right)}}-\underset{\text{washout}}{\underbrace{Dy\left(t\right)}}\hfill \end{array}$$

(4)

where $m,C_{th2}$ are the Hill coefficient and the threshold for Lactonase activation, respectively, and δ is the Lactonase-dependent degradation rate of AHLs. Observe that there is no outflow term in the complex equation, as PpuR–AHL is considered to be intracellular.

The model (4) can be extended by adding one equation for PpuR dynamics in one cell, as in [17] or in system (1). Then the system reads

$$\begin{array}{cc}\hfill \dot{S}\left(t\right)& =D{S}_0-{\gamma }_{S}N\left(t\right)\frac{S{\left(t\right)}^{n_s}}{K_m^{n_s}+S{\left(t\right)}^{n_s}}-DS\left(t\right)\hfill \\ \hfill \dot{N}\left(t\right)& =aN\left(t\right)\frac{S{\left(t\right)}^{n_s}}{K_m^{n_s}+S{\left(t\right)}^{n_s}}-DN\left(t\right)\hfill \\ \hfill \dot{x}\left(t\right)& =\left({\alpha }_A+{\beta }_A\frac{C{\left(t\right)}^n}{C_{th}^n+C{\left(t\right)}^n}\right)N\left(t\right)-{\gamma }_{A}x\left(t\right)-{\pi }_1^{+}R\left(t\right)x\left(t\right)+{\pi }_1^{-}C\left(t\right)-Dx\left(t\right)-\delta x\left(t\right)y\left(t\right)\hfill \\ \hfill \dot{C}\left(t\right)& ={\pi }_1^{+}R\left(t\right)x\left(t\right)-{\pi }_1^{-}C\left(t\right)\hfill \\ \hfill \dot{R}\left(t\right)& ={\alpha }_R+{\pi }_1^{-}C\left(t\right)-{\pi }_1^{+}R\left(t\right)x\left(t\right)-{\gamma }_{R}R\left(t\right)\hfill \\ \hfill \dot{y}\left(t\right)& ={\alpha }_L\frac{C{(t-\mathsf{\tau })}^m}{C_{th2}^m+C{(t-\mathsf{\tau })}^m}N\left(t\right)-{\gamma }_{L}y\left(t\right)-Dy\left(t\right)\hfill \end{array}$$

(5)

2.1.5. Reduced Model

When being interested in the long term behavior of regulatory systems in the chemostat, one can assume that substrate concentration and bacterial density have approximately assumed a stationary state $(N^{*},S^{*})$. We consider the system (5) for large values of t and impose quasi-steady state conditions for PpuR and complex. In other words, we assume that when bacteria stay at their saturation level, the dynamics of R and C is slow compared to those of AHL and Lactonase. The equilibrium conditions are given by

$$R^{*}=\frac{{\alpha }_R}{{\gamma }_R},C^{*}=\underset{=:\tilde{\gamma }}{\underbrace{\frac{{\pi }_1^{+}{\alpha }_R}{{\pi }_1^{-}{\gamma }_R}}}x=\tilde{\gamma }x$$

(6)

Define the parameters

$$\begin{array}{cccc}\alpha ={\alpha }_A{N}^{*},& \beta ={\beta }_A{N}^{*},& x_{th}=C_{th}/\tilde{\gamma },& \omega ={\gamma }_L+D\\ \gamma ={\gamma }_A+D,& \rho ={\alpha }_L{N}^{*},& y_{th}=C_{th2}/\tilde{\gamma }\end{array}$$

(7)

Substituting the equilibrium conditions (6) into (5), we obtain the system

$$\begin{array}{cc}\hfill \dot{x}\left(t\right)& =\alpha -\gamma x\left(t\right)-\delta x\left(t\right)y\left(t\right)+\beta \frac{x{\left(t\right)}^n}{x_{th}^n+x{\left(t\right)}^n}\hfill \\ \hfill \dot{y}\left(t\right)& =\rho \frac{x{(t-\mathsf{\tau })}^m}{y_{th}^m+x{(t-\mathsf{\tau })}^m}-\omega y\left(t\right)\hfill \end{array}$$

(8)

Observe that all parameter values are non-negative. Their meaning is summarized in

Table 1. Variables and parameters in model (8), with values used for data fit in Figure 3.

Symbol	Description	Value (Unit)	Comments/Source
$N^{*}$	Cell density at equilibrium	$4.5929\times {10}^{11}$ (cells/lit)	[18]
α	Basic AHL production rate	$1.0564\times {10}^{-7}$ (mol/(lit${}^2·$ h))	$={\alpha }_A*N_{equi}$, [18]
γ	AHL decay rate (includes washout)	0.105 (1/h)	$={\gamma }_A+D$, [18]
δ	Lactonase-dependent degradation rate	$1.5000\times {10}^{-4}$ (lit/(mol · h))	[18]
β	Feedback-regulated AHL production rate	$1.0564\times {10}^{-6}$ (mol/(lit${}^2·$ h))	$={\beta }_A*N_{equi}$, [18]
n	Hill coefficient for x	2.3 (dimensionless)	[18]
$x_{th}$	Critical threshold for positive-feedback in x	$3.597\times {10}^{-13}$ (mol/lit)	estimated
ω	Lactonase decay rate (includes washout)	0.105 (1/h)	$={\gamma }_e+D$, [18]
ρ	Lactonase production rate	$5.0521\times {10}^3$ (mol/(lit${}^2·$ h))	$={\alpha }_e*N_{equi}$, [18]
τ	Delay in the release of y	2 (h)	[18]
m	Hill coefficient for x	2.5 (dimensionless)	[18]
$y_{th}$	Critical threshold for positive-feedback in y	$3.597\times {10}^{-13}$ (mol/lit)	estimated
$x_0\left(t\right)$	AHL concentration (initial data) $t\in [18,20]$	$5.4044\times {10}^{-7}$ (mol/lit)	mean of exp. data
$y_0\left(t\right)$	Lactonase (initial data) $t\in [18,20]$	$5.2\times {10}^3$ (mol/lit)	estimated

.

2.2. Experimental Data

We report experimental data as published in the previous publication [18]. Pseudomonas putida IsoF was cultivated and grown in a continuous culture with a working volume of $2\mathrm{L}$, under controlled conditions at $30{}^{\circ }\mathrm{C}$, enabling the reproducible establishment of defined environmental conditions.

AHL molecules and their degradation products were identified and quantified via two different methods. The first one is the so-called ultra-high-performance liquid chromatography (UHPLC), a technique used to separate different components in a mixture. The second method, the enzyme-linked immunosorbent assay (ELISA), allows the rapid detection and quantification of AHLs and Homoserines directly in biological samples with the help of antibodies.

2.3. Parameter Estimation

In [18], the model (4) was fitted to a first set of experimental data using a mean square error algorithm and the simplex search algorithm in MATLAB${}^{\text{®}}$ (Version 2013b, The Mathworks, Natick, MA, USA, 2013). Obtained parameter values were used to validate further data sets with minor adaptations for some initial values, which increased the quality of the fit.

Starting from these estimated parameter values, we fit the reduced system (8) to the same experimental data published in [18]. The fit was performed using curve fitting tools in MATLAB${}^{\text{®}}$ and Wolfram Mathematica${}^{\text{®}}$ (Version 10, Wolfram Research, Champaign, IL, USA, 2014). The reduced model (8) is obtained assuming the cell population to be in equilibrium; that is, it holds only for times $t>t_{ec}$, where $t_{ec}$ is the time at which the cell population has reached its saturation level.

3. Results

In this section we present analytical results concerning qualitative properties of the solution of the reduced model (8), as well as numerical simulations and data fit.

3.1. Existence of Solutions

Theorem 1. 

Let the system (8) hold for $t\ge t_0$, and let initial data $x\left(t\right)=x_0\left(t\right),y\left(t\right)=y_0\left(t\right)$ be given for $t\in [t_0-\mathsf{\tau },t_0]$, $\mathsf{\tau }>0$, with $x_0,y_0$ Lipschitz continuous. Then there is a unique solution to (8) in $[t_0,\infty )$. Moreover, if $x_0,y_0$ are non-negative, the solution is also non-negative.

Proof. 

The proof follows from basic principles of DDE theory, cf. [19,24,25]. We provide here a sketch of the proof steps. For simplicity, we shall denote the right-hand side of the system (8) by $f(u,v)$, where $u=\left(x\right(t),y(t\left)\right)$ and $v=\left(x\right(t-\mathsf{\tau }),y(t-\mathsf{\tau }\left)\right)$.

Local existence. 

For the construction of a local (maximal) solution on an interval $[t_0,t_0+\Delta ),\Delta >0$, it is sufficient to guarantee Lipschitz continuity of the initial data, as well as of f with respect to both arguments, cf. [25] (Thm. 2.2.1). It is easy to verify that the right-hand side of (8) is continuously differentiable with respect to the delayed, as well as to the non-delayed argument, and that the partial derivatives are bounded (computation not shown).

Non-negativity. 

Preservation of positivity is due to the fact that the delay only appears in the positive feedback term. Indeed, if for some $\overline{t}>t_0$, $x\left(\overline{t}\right)=0$ then $\dot{x}\left(\overline{t}\right)=\alpha >0$, and $x\left(t\right)$ remains non-negative. With this result it follows that also y stays non-negative. If for some $\overline{t}>t_0$, $y\left(\overline{t}\right)=0$, then $\dot{y}\left(\overline{t}\right)=\rho \frac{x{(\overline{t}-\mathsf{\tau })}^m}{y_{th}^m+x{(\overline{t}-\mathsf{\tau })}^m}\ge 0$.

Global existence. 

We show that the maximal solution is bounded. This follows with estimates on the right-hand side. Observe that

$$\begin{array}{cc}\hfill \dot{y}\left(t\right)& =\rho \frac{x{(t-\mathsf{\tau })}^m}{y_{th}^m+x{(t-\mathsf{\tau })}^m}-\omega y\left(t\right)\hfill \\ & \le \rho -\omega y\left(t\right)\hfill \end{array}$$

hence for all $t\ge t_0$ we have $0\le y\left(t\right)\le \widehat{y}$, where $\widehat{y}:=\left(y_0\left(t_0\right)-\frac{\rho }{\omega }\right)e^{-\omega t}+\frac{\rho }{\omega }.$

Similarly,

$$\begin{array}{cc}\hfill \dot{x}\left(t\right)& =\alpha -\gamma x\left(t\right)-\delta x\left(t\right)y\left(t\right)+\beta \frac{x{\left(t\right)}^n}{x_{th}^n+x{\left(t\right)}^n}\hfill \\ & \le (\alpha +\beta )-(\gamma +\delta y(t\left)\right)x\left(t\right)\hfill \\ & \le (\alpha +\beta )-\gamma x\left(t\right)\hfill \end{array}$$

Thus for all $t\ge t_0$ we have $0\le x\left(t\right)\le \widehat{x}$, with $\widehat{x}:=\left(x_0\left(t_0\right)-\frac{\alpha +\beta }{\gamma }\right)e^{-\gamma t}+\frac{\alpha +\beta }{\gamma }.$ The maximal solution is bounded, hence it exists on $[t_0,\infty )$, cf. [25] (Thm. 2.2.2).

3.2. Fixed Points

Fixed points of (8) are given by the solutions of

$$\left\{\begin{array}{cc}0\hfill & =\alpha -\gamma \overline{x}-\delta \overline{x}\overline{y}+\beta \frac{{\overline{x}}^n}{x_{th}^n+{\overline{x}}^n}\hfill \\ 0\hfill & =\rho \frac{{\overline{x}}^m}{y_{th}^m+{\overline{x}}^m}-\omega \overline{y}.\hfill \end{array}\right.$$

So we have

$$\overline{y}=\frac{\rho }{\omega }\frac{{\overline{x}}^m}{y_{th}^m+{\overline{x}}^m}$$

where $\overline{x}$ is given by the solutions of

$$\alpha -\gamma \overline{x}-\frac{\delta \rho }{\omega }\frac{{\overline{x}}^{m+1}}{y_{th}^m+{\overline{x}}^m}+\beta \frac{{\overline{x}}^n}{x_{th}^n+{\overline{x}}^n}=0$$

(9)

Recall that for the biological motivation of the model, we are only interested in non-negative $\overline{x}$. In the following, for simplicity of notation, we shall omit the bars from $\overline{x}$.

In the general case $n\ne m$ and $x_{th}\ne y_{th}$, solutions of (9) are the zeros of the polynomial

$$a_0{x}^{n+m+1}+a_1{x}^{n+m}+a_2{x}^{n+1}+a_3{x}^{m+1}+a_4{x}^n+a_5{x}^m+a_{6}x+a_7=0$$

where

$$\begin{array}{cc}\hfill a_0=-(\gamma \omega +\rho \delta )<0,& a_1=\omega (\alpha +\beta )>0,\hfill \\ \hfill a_2=-\gamma y_{th}^m\omega <0,& a_3=-x_{th}^n(\gamma \omega +\rho \delta )<0,\hfill \\ \hfill a_4=y_{th}^m\omega (\alpha +\beta )>0,& a_5=\alpha \omega x_{th}^n>0,\hfill \\ \hfill a_6=-\omega \gamma x_{th}^n{y}_{th}^m<0,& a_7=\alpha \omega x_{th}^n{y}_{th}^m<0.\hfill \end{array}$$

Let us consider a special case which is relevant for our application, and assume $n=m=2$ and $x_{th}=y_{th}$. Then, fixed points $\left(\overline{x},\overline{y}\right)$ satisfy

$$\overline{y}=\frac{\rho }{\omega }\frac{{\overline{x}}^2}{x_{th}^2+{\overline{x}}^2}$$

(10)

with $\overline{x}$ given by the solutions of a cubic equation

$$(\delta \rho +\omega \gamma ){\overline{x}}^3-\omega (\beta +\alpha ){\overline{x}}^2-\omega \gamma x_{th}^2\overline{x}-\alpha \omega x_{th}^2=0$$

(11)

which has either three real zeros, or one real and two complex solutions. Thus, we might have up to three biologically-relevant fixed points.

3.3. The case $\mathsf{\tau }=0$

Consider the ODE system obtained from (8) by setting $\mathsf{\tau }=0$:

$$\begin{array}{cc}\hfill \dot{x}\left(t\right)& =\alpha -\gamma x\left(t\right)-\delta x\left(t\right)y\left(t\right)+\beta \frac{x{\left(t\right)}^n}{x_{th}^n+x{\left(t\right)}^n}\hfill \\ \hfill \dot{y}\left(t\right)& =\rho \frac{x{\left(t\right)}^m}{y_{th}^m+x{\left(t\right)}^m}-\omega y\left(t\right),\hfill \end{array}$$

(12)

It is important to know the dynamics of (12), because for small delays ($\mathsf{\tau }>0$), the DDE system (8) will very likely behave as the ODE system (12), cf. [24].

Observe that the ODE system (12) and the DDE system (8) have exactly the same equilibrium points. In general, a DDE system and the associated ODE system have the same number of fixed points, but if the delay appears in the coefficients, the fixed points of the DDE system could be shifted with respect to those of the ODE system.

The presence of a negative feedback in (12) leads to the hypothesis that oscillatory solutions might show up. We investigate local properties of the steady states, looking for Hopf-bifurcations. For linear (local) stability analysis, we compute the Jacobian matrix of system (12),

$$J=\left(\begin{array}{cc}-\gamma -\delta \overline{y}+\beta \frac{n{x}_{th}^n{\overline{x}}^{n-1}}{{(x_{th}^n+{\overline{x}}^n)}^2}& -\delta \overline{x}\\ \rho \frac{m{y}_{th}^m{\overline{x}}^{m-1}}{{(y_{th}^m+{\overline{x}}^m)}^2}& -\omega \end{array}\right).$$

In the special case $n=m=2$ and $x_{th}=y_{th}$, we have

$$J=\left(\begin{array}{cc}-\gamma -\delta \overline{y}+\beta \frac{2x_{th}^2\overline{x}}{{(x_{th}^2+{\overline{x}}^2)}^2}& -\delta \overline{x}\\ \rho \frac{2x_{th}^2\overline{x}}{{(x_{th}^2+{\overline{x}}^2)}^2}& -\omega \end{array}\right)$$

(13)

The trace and the determinant of (13) at a stationary point at $(\overline{x},\overline{y})$, with $\overline{y}$ in (10), are given by

$$\begin{array}{ccc}\hfill Tr\left(J\right)& =& \frac{2\beta x_{th}^2\overline{x}}{{(x_{th}^2+{\overline{x}}^2)}^2}-\gamma -\delta \frac{\rho }{\omega }\frac{{\overline{x}}^2}{x_{th}^2+{\overline{x}}^2}-\omega \hfill \\ \hfill det\left(J\right)& =& -\omega \left(\frac{2\beta x_{th}^2\overline{x}}{{(x_{th}^2+{\overline{x}}^2)}^2}-\gamma -\delta \frac{\rho }{\omega }\frac{{\overline{x}}^2}{x_{th}^2+{\overline{x}}^2}\right)+\delta \overline{x}\frac{2\rho x_{th}^2\overline{x}}{{(x_{th}^2+{\overline{x}}^2)}^2}.\hfill \end{array}$$

If there was only one stationary point and this one is a repellor, one can use the fact that all solutions are bounded and stay positive, then the Poincare–Bendixson theorem yields the existence of periodic solutions. For a Hopf-bifurcation, necessary conditions are $Tr\left(J\right)=0$ and $\Delta \left(J\right)=Tr{\left(J\right)}^2-4det\left(J\right)<0$. We choose δ, the Lactonase activity, as bifurcation parameter. From the trace condition, we get

$$\delta \left[\frac{\rho }{\omega }{\overline{x}}^2(x_{th}^2+{\overline{x}}^2)\right]-2\beta x_{th}^2\overline{x}+\gamma {(x_{th}^2+{\overline{x}}^2)}^2+\omega {(x_{th}^2+{\overline{x}}^2)}^2=0.$$

We solve for δ and obtain

$$\begin{array}{ccc}\hfill \delta =\delta \left(\overline{x}\right)& =& \frac{2\beta x_{th}^2\overline{x}-(\gamma +\omega ){(x_{th}^2+{\overline{x}}^2)}^2}{\frac{\rho }{\omega }{\overline{x}}^2(x_{th}^2+{\overline{x}}^2)}\hfill \\ \hfill & =& \frac{-\omega (\gamma +\omega )(x_{th}^2+{\overline{x}}^2)}{\rho {\overline{x}}^2}+\frac{2\beta x_{th}^2\omega }{\rho \overline{x}(x_{th}^2+{\overline{x}}^2)}.\hfill \end{array}$$

(14)

Note that, in turn, $\overline{x}$ also depends on δ, cf. Equation (11). Neglecting this for a minute, we observe that $\underset{\overline{x}\to \infty }{lim}\delta \left(\overline{x}\right)=-(\gamma +\omega )\frac{\omega }{\rho }<0$, whereas $\underset{\overline{x}\to 0}{lim}\delta \left(\overline{x}\right)\to +\infty$. Due to the intermediate value theorem, there exists a $\tilde{x}>0$, such that $\delta \left(\overline{x}\right)>0$ for $\overline{x}>\tilde{x}$. It is possible to choose $\tilde{x}$ as the smallest positive solution of $(\gamma +\omega ){(x_{th}^2+{\overline{x}}^2)}^2>2\beta x_{th}^2\overline{x}$. If a $\overline{x}=\overline{x}\left(\delta \right)>0$ satisfies this condition, then $Tr\left(J\right)=0$ at $\left(\overline{x},\overline{y}\right)$.

In the next step, we check the discriminant condition ($\Delta \left(J\right)<0$), or equivalently, $det\left(J\right)>0$, as for the Hopf-bifurcation we need simultaneously $Tr\left(J\right)=0$.

$$\begin{array}{ccc}\hfill det\left(J\right)& =& -\omega (Tr\left(J\right)+\omega )+\delta \frac{2\rho x_{th}^2{\overline{x}}^2}{{(x_{th}^2+{\overline{x}}^2)}^2}\hfill \\ & =& -{\omega }^2+\delta \frac{2\rho x_{th}^2{\overline{x}}^2}{{(x_{th}^2+{\overline{x}}^2)}^2}\hfill \\ & =& ...\hfill \\ & =& {\omega }^2\left[\frac{-3x_{th}^2-{\overline{x}}^2}{x_{th}^2+{\overline{x}}^2}\right]+\omega \frac{2{\overline{x}}^2}{{(x_{th}^2+{\overline{x}}^2)}^3}\left[2\beta \overline{x}{x}_{th}^2-\gamma {(x_{th}^2+{\overline{x}}^2)}^2\right].\hfill \end{array}$$

We solve $det\left(J\right)=0$ in dependence of the Lactonase decay rate $\omega >0$, and find the roots ${\omega }_1=0$ (which does not provide further information), and

$${\omega }_2=\frac{\frac{2x_{th}^2}{{(x_{th}^2+{\overline{x}}^2)}^2}\left[2\beta \overline{x}{x}_{th}^2-\gamma {(x_{th}^2+{\overline{x}}^2)}^2\right]}{3x_{th}^2+{\overline{x}}^2}.$$

Hence $det\left(J\right)>0$ when ${\omega }_2>0$. We need to distinguish between two cases. If $\overline{x}>x_{th}$, i.e., a stationary state with high AHL concentration and activated bacteria, then we get

$$2\beta \overline{x}{x}_{th}^2-\gamma {(x_{th}^2+{\overline{x}}^2)}^2>2\beta x_{th}^3-\gamma {\left(2{\overline{x}}^2\right)}^2=2\beta x_{th}^3-4\gamma {\overline{x}}^4.$$

Thus, if $2\beta x_{th}^3-4\gamma {\overline{x}}^4>0$, then ${\omega }_2>0$. Analogously, we get ${\omega }_2>0$ if $2\beta {\overline{x}}^3-4\gamma x_{th}^4>0$, in case of $\overline{x}<x_{th}$, i.e., with bacteria in a non-activated quorum sensing state. All in all, if the model parameters and the resulting stationary point satisfy the last condition yielding ${\omega }_2>0$, and simultaneously $\delta >0$ according to (14), then a Hopf-bifurcation takes place.

3.4. The Case $\mathsf{\tau }>0$

We are interested in stability switches due to the presence of a delay $\mathsf{\tau }>0$ in (8). Consider the case $n=m=2$ and $x_{th}=y_{th}$, and let $(\overline{x},\overline{y})$ be one equilibrium point of (8). The linearized system about $(\overline{x},\overline{y})$ is given by

$$\dot{Z}\left(t\right)=AZ\left(t\right)+BZ(t-\mathsf{\tau }),$$

(15)

with

$$Z\left(t\right)=\left(\begin{array}{c}{z}_1\left(t\right)\\ z_2\left(t\right)\end{array}\right),A=\left(\begin{array}{cc}a& b\\ 0& d\end{array}\right),B=\left(\begin{array}{cc}0& 0\\ c& 0\end{array}\right),$$

and

$$\begin{array}{cc}\hfill a& =-\gamma -\delta \overline{y}+2\beta x_{th}^2\frac{\overline{x}}{\left(x_{th}^2+{\overline{x}}^2\right)},\hfill \\ \hfill b& =-\delta \overline{x}\le 0,\hfill \\ \hfill c& =2\rho x_{th}^2\frac{\overline{x}}{\left(x_{th}^2+{\overline{x}}^2\right)}\ge 0\hfill \\ \hfill d& =-\omega <0.\hfill \end{array}$$

(16)

The characteristic equation corresponding to (15) is given by

$$det\left(\lambda \mathbb{I}-A-B{e}^{-\lambda \mathsf{\tau }}\right)=0,$$

or equivalently,

$${\lambda }^2-\lambda (a+d)+ad-bc{e}^{-\lambda \mathsf{\tau }}=0.$$

(17)

Characteristic equations of this and more general type have been studied in [24]. In the following, we report results from [24], adapting them to our specific example. We apply standard methods for the analysis of characteristic equations and switches with respect to increasing delays, hence we consider purely imaginary roots, $\lambda =i\phi ,\phi >0$. Separating real and imaginary parts in (17) we obtain

$$\begin{array}{cc}\hfill {\phi }^2-ad& =-bccos\left(\phi t\right)\hfill \\ \hfill \phi (a+d)& =bcsin\left(\phi t\right).\hfill \end{array}$$

Now we square left- and right-hand sides and sum up the two equations, obtaining

$${\phi }^4+{\phi }^2(a^2+d^2)+a^2{d}^2=b^2{c}^2.$$

(18)

Its roots are

$${\phi }_{\pm }^2=\frac{1}{2}\left(-(a^2+d^2)\pm \sqrt{{(a^2-d^2)}^2+4b^2{c}^2}\right).$$

It can be seen from (18) that the parabola in ${\phi }^2$ is open upwards and it has:

• no positive intercept with the horizontal axis, if $a^2{d}^2-b^2{c}^2>0$, i.e., if $\left|ad\right|>-bc$;
• one positive intercept (${\phi }_{+}$) with the horizontal axis, if $\left|ad\right|<-bc$.

In the first case, there is no stability switch with respect to τ; that is, the stability of the equilibrium point $(\overline{x},\overline{y})$ remains the same for any $\mathsf{\tau }\ge 0$, and it is sufficient to study the ODE case (12). In the case $\left|ad\right|<-bc$, there is one root (${\phi }_{+}$) with positive imaginary part, hence one stability switch. In order to find out in which direction the stability switch occurs, we study the sign of the real part $\Re \lambda \left(\mathsf{\tau }\right)$ in $\lambda =i{\phi }_{+}$, for $\mathsf{\tau }>0$. From (17) we have

$$\left\{2\lambda -(a+d)+\mathsf{\tau }bc{e}^{-\lambda \mathsf{\tau }}\right\}\frac{d\lambda \left(\mathsf{\tau }\right)}{d\mathsf{\tau }}=\lambda bc{e}^{-\lambda \mathsf{\tau }}.$$

It follows

$$\begin{array}{cc}\hfill sign{\left\{\frac{d\Re \lambda \left(\mathsf{\tau }\right)}{d\mathsf{\tau }}\right\}}_{\lambda =i{\phi }_{+}}& =sign{\left\{\Re {\left(\frac{d\lambda }{d\mathsf{\tau }}\right)}^{-1}\right\}}_{\lambda =i{\phi }_{+}}\hfill \\ & =sign{\left\{\Re \left(\frac{2\lambda -(a+d)}{\lambda ({\lambda }^2-(a+d)\lambda +ad)}\right)\right\}}_{\lambda =i{\phi }_{+}}\hfill \\ & =sign\left\{\frac{2({\phi }_{+}^2-ad)+{(a+d)}^2}{({\phi }_{+}^2{(a+d)}^2+{(ad-{\phi }_{+}^2)}^2}\right\}\hfill \\ & =sign\left\{2({\phi }_{+}^2-ad)+{(a+d)}^2\right\}\hfill \\ & =sign\left\{2{\phi }_{+}^2+a^2+d^2\right\}\hfill \\ & =+1.\hfill \end{array}$$

Roots cross the imaginary axis from the left to the right, indicating stability loss. If the solution $(\overline{x},\overline{y})$ is asymptotically stable for $\mathsf{\tau }=0$, then it is uniformly asymptotically stable for all $\mathsf{\tau }<{\mathsf{\tau }}_c$ and unstable for $\mathsf{\tau }>{\mathsf{\tau }}_c$, where

$${\mathsf{\tau }}_c=\frac{{\theta }_c}{{\phi }_{+}},$$

(19)

with ${\theta }_c$ implicitly defined by

$$arctan\left({\theta }_c\right)=\frac{(a+d){\phi }_{+}}{ad-{\phi }_{+}^2}.$$

(20)

All in all, we have shown the following result.

Theorem 2. 

Let $(\overline{x},\overline{y})$ be one equilibrium point of (8), with $\mathsf{\tau }>0$, $n=m=2$ and $x_{th}=y_{th}$. Assume that $\left|ad\right|<bc$, with $a,b,c,d$ given in (16). Then, the equilibrium point is uniformly asymptotically stable for all $0<\mathsf{\tau }<{\mathsf{\tau }}_c$ and unstable for $\mathsf{\tau }>{\mathsf{\tau }}_c$, with ${\mathsf{\tau }}_c$ defined by (19)–(20).

3.5. Numerical Simulations and Data Fitting

We consider experimental data published in [18] and perform numerical simulations in MATLAB${}^{\text{®}}$ and Wolfram Mathematica${}^{\text{®}}$. The reduced model (8) is obtained assuming the cell population to be in equilibrium, that is, for times $t>t_{ec}$, where $t_{ec}$ is the time at which bacteria have reached the saturation level. In Figure 2, we read from experimental data that the cell population reaches the equilibrium after ca. 20 h from the beginning of the experiment. Hence, we take $t_{ec}=20$ as the starting time point for numerical simulations of the reduced system (8), and define initial data on the time interval $[t_{ec}-\mathsf{\tau },t_{ec}]$.

For simplicity, we assume that initial data are constant functions on the definition interval, see also [17,18]. We fix the value of the delay, $\mathsf{\tau }=2$ h, as in [18]. Then we take $x\left(t\right)={\widehat{x}}_{19}$, for $t\in [18,20]$, ${\widehat{x}}_{19}$ being the mean value of ELISA and UHPLC measurements at 19 h from the beginning of the experiment. Initial data for the Lactonase are estimated from simulations of the full model (4) in [18]. To date, there is no experimental data available for Lactonase concentration, thus parameters associated with Lactonase production (ρ), decay (ω), and activity (δ) can be only estimated from AHL experimental data. This means in turn that there are several plausible solutions for the estimation of $\rho ,\omega$, and δ. We choose to maintain parameter values as estimated in [18].

It can be seen from the model reduction assumptions, as well as from the simplified parameters (7) that we lack information on the receptor production (${\alpha }_R$) and decay (${\gamma }_R$); indeed, there is no equation for PpuR in (4). These parameters play a role for the critical threshold value, $x_{th}$, in the complex-regulated processes. We fit $x_{th}$ and $y\left(t\right)=y_0,t\in [18,20]$, fixing all other parameter values as in [18]. The results are summed up in Table 1, with parameter values as provided by the fitting procedure, without rounding. In Figure 3 we show a comparison between the numerical solution of model (4), that of the simplified model (8) and experimental data for AHL time series.

With the estimated parameter values in Table 1, we consider the analytical results in Section 3.2 and Section 3.4. The system has three equilibrium points ${\overline{x}}_1,{\overline{x}}_2,{\overline{x}}_3$, but we only consider the stability properties of the largest one $(x_3,y_3)=(1.593\times {10}^{-7},4.809\times {10}^4)$, which corresponds to high AHL level and to an activated state of the bacteria population. Parameter values satisfy $a^2{d}^2>b^2{c}^2$, thus there is no stability switch with respect to τ, and the system behaves as in the case $\mathsf{\tau }=0$. We go back to Section 3.3 and consider the Jacobian matrix (13), obtaining $tr\left(J\right)=-7.4243$ and $det\left(J\right)=0.7685$. Hence, with the estimated parameter values, the system (8) has a locally asymptotically stable equilibrium $(x_3,y_3)$ in which bacteria are activated.

4. Discussion

In this paper we have introduced a system of two delay differential Equation (8) for quorum sensing of Pseudomonas putida in a continuous culture. Motivated by experimental data, a more detailed mathematical model (4) was previously proposed in [18]. Though the system (4) describes the regulatory network in greater detail, in the long run, bacteria reach a saturation level and the model can be reduced to two governing Equation (8), as we have shown in Section 2.1.

Surprisingly, even a simple model such as (8) can be used to explain experimental data (Figure 3), maintaining parameter values from a previous fit [18] for almost all model parameters. However, one should take into account that this is valid only from the moment the bacteria population has reached its saturation level. If one is interested in understanding quorum sensing in the initial phases (lag and exponential phase) of bacterial population, then it is convenient to use a more detailed model, such as (4) or (5).

The advantage of system (8) is that it can be investigated thoroughly thanks to well-established methods. We have shown existence and uniqueness of solutions and, more importantly, we could guarantee preservation of positivity. This property is often violated in systems of delay differential equations. We have studied linearized stability of non-negative equilibria and proved that the delay system (8) might show stability switches as the delay increases. On the other side, the Lactonase activity (δ) can induce Hopf-bifurcation in the associated ODE system (12).

For simplicity of computation in the analysis of the system, we have considered only the case of small Hill coefficients ($n=m=2$), which corresponds to a maximum of three biologically-relevant stationary states. This assumption is, however, not as restrictive as it seems. Three stationary states, with an intermediate unstable one, are the basis for the bistability situation already discovered in analogous regulatory networks [11,17]. Moreover, similar small values for the Hill coefficients were found to fit experimental data (Table 1) and correspond well to the biological assumption of a dimer being relevant for the positive feedback in the quorum sensing system of Pseudomonas putida [16]. With the fitted parameter values, we do not find delay-induced stability switches. This is not a hint of the delay being not relevant. Though a positive time lag might not change the main qualitative behavior of the system, the DDE model still describes the experimentally-determined data and their time course quantitatively better than the associated ODE system, in particular when the bacteria population is in the lag or exponential phase. Stability switches and periodic oscillatory behavior might appear for a different choice of the parameters in system (8). As the main focus of this work was to provide a description for a real biological process, we decided to omit further numerical investigation on the qualitative behavior of (8).

Delay equations have been previously used in mathematical models of continuous cultures. Commonly, a time lag was included to describe the time necessary for the bacteria to convert nutrients in new biomass [26,27]. Being interested in the long term dynamics with bacteria being at an equilibrium, we have chosen not to consider such reproduction lags in our model. In our case, the time lag arises from the dynamics of the regulatory network, in particular from the initialization processes of the AHL-degrading enzyme.

Taken together, the presented simplified delay equation system is a good compromise between refined modeling for a well-known gene regulatory network with several players, and a system of equations which still allows explicit analysis of the basic qualitative behavior as well as parameter determination from few experimental data.