Dynamical behaviors of quorum sensing network mediated by combinatorial perturbation

1.   Introduction

Quorum sensing (QS) is a signaling system that senses bacterial population density and coordinates gene expression by secreting small diffusible autoinducer signaling molecules ^[1,2,3,4,5]. Autoinducers are also well-known signaling molecules, which can diffuse into the extracellular space and provide a communication mechanism for cell aggregation ^[6,7]. Acyl-homoserine lactone (AHL) is a commonplace autoinducer chemical signaling molecule in the QS system, which consists of a homoserine lactone (HSL) loop carrying acyl chains of length C4 to C18 ^[8] and can be produced by LuxI enzymatically ^[9]. Moreover, there are many differences in the production and diffusion of AHL autoinducer between low cell density (LCD) and high cell density (HCD). At LCD, LuxI is expressed at a basal level. The autoinducer AHL synthesized as a result of LuxI expression does not accumulate in the cell, but instead rapidly diffuses out. When population density rises (HCD), AHL accumulates in the cell owing to the lowered diffusion gradient across the cell membrane ^[10]. The detection of the minimum stimulation concentration threshold of AHL will cause changes in gene expression. As we know, bacteria control the gene expression in response to changes in cell density and species complexity through the QS mechanism. In many cases, bacterial communities can recognize themselves and foreign objects, which is very beneficial, especially in symbiosis, niche adaptation, and production of secondary metabolites ^[11]. When the bacterial colony perform these beneficial activities, they will work cooperatively and synchronously ^[12].

Dynamical behaviors are commonplace in QS system, which controls various aspects of biological events ^{[13,14,15,16]}. Extensive evidences show that the urease PH value of plants and microorganisms and drug delivery depend on the oscillation and bistable states of the QS system. Quorum sensing of plants and microorganisms is accomplished by the diffusion of enzyme particles. The diffusion of different enzyme loads corresponds to different states of quorum sensing. The PH exhibits an s-shaped switch when high enzyme load, oscillations when intermediate enzyme load, and a hysteresis switch when low enzyme load ^[17]. Furthermore, many scientists are committed to building a QS circuit for purpose to prompt bacterial populations to carry medicine to treat disease ^[18,19]. Research has shown that minor changes in the amplitude and period of the QS circuit can severely affect the dose and time of the bacteria spilling the drug on the target ^[20].

Previously, some Repressilator models (TetR CI LacI) have been proposed ^{[21,22,23,24,25,26,27,28]}, which showed some oscillatory dynamics from a mathematical point of view. However, none of these models explain which phenomena may be affected by changes in the concentration of some important substances in the QS system. In addition, some evidences indicate that LuxI, AiiA and H$_2$O$_2$ are key substances involved in quorum sensing ^{[29,30,31,32]}. In addition, a problem often encountered is that there may be a wide range of time intervals in the process of gene expression. However, in previous studies, for the sake of simplicity, many models implicitly assumed that the gene expression during the QS regulation process was transient. In this regard, it is important to note that most gene expression time delays are estimated to be about 15–27 minutes, including about 10–20 minutes for transcription, about 4 minutes for transport out of the nucleus, and about 1–3 minutes for translation ^[33,34]. Therefore, ignoring the impact of such a time delay in the QS module is disadvantageous ^[35]. However, considering the effects of combinatorial perturbations to LuxI, AiiA, H$_2$O$_2$, and delay on the system dynamics is not easily achieved due to the complexity of the network. Moreover, most studies fix the number of multi-cells and analyze the regulation effect of internal elements on a certain protein ^{[23,24,25,28]}. Thus, here we chose a research approach similar to ^[25] to consider the dynamical effects of the co-regulation of these three substances and delays on single cell. Therefore, a new single cell computational model containing these substances and delay is established on the basis of the original model ^[28]. To simplify the analysis, all time delays are assumed to be the same. When they are different, similar analysis can be performed.

To sum up, the novelty of this paper is mainly reflected in the following aspects. First, the individual effect of time delay on the QS system is investigated. The stability and oscillations of the QS system are studied for the system with and without delay. Secondly, the individual effects of the three substances AiiA, LuxI, and H$_2$O$_2$ on the QS system dynamics are discussed. Furthermore, based on two-parameter bifurcation analysis, the impact of the combinatorial perturbation of AiiA, LuxI and H$_2$O$_2$ on oscillatory behaviors of the QS system is performed. Finally, the influence of combinatorial of perturbation of AiiA, LuxI, H$_2$O$_2$ and delay on QS system is analyzed. That is, the sensitivity of the delay bifurcation value $\tau_{0}$ to AiiA, LuxI and H$_2$O$_2$ is systematically discussed. Interestingly, the effects of changes in the concentration of these internal substances on the drug release and the PH value is also revealed.

2.   Dynamic modelling

In this paper, we use the dimensional repressilator model (LacI, CI, TetR) and AHL as the original model ^[28]. To capture the important impact of the transcriptional processes on the QS system, we convert them into a time delay. Crucially, the time delay of protein synthesis is objective, which is a complex and time-consuming process ^[36]. In addition, AiiA and LuxI are mainly used to inhibit and generate AHLs, which are signal molecules for communication between bacterial populations ^[29,30]. In comparison, H$_2$O$_2$ is relatively rare in quorum sensing interaction networks. However, studies have shown that H$_2$O$_2$ has a significant effect on inducing the production of AiiA and LuxI ^[32]. To reveal their key roles in dynamics and regulatory mechanisms, we incorporate them into the QS network and analyze their effects on dynamics and biological phenomena. The new QS circuit is shown in Figure 1. Based on these descriptions, the mathematical model is construed as follows:

Here, $x$, $y$, $z$, $w$, $v$, $p$, and $s$ represent the concentrations of TetR, CI, LacI, AiiA, LuxI, AHL, and H$_2$O$_2$ at time $t$, respectively. First, protein levels of TetR, CI, and LacI are explained. The nonlinear functions $\frac{eK_m^n}{\delta_m(K_m^n+(z(t-\tau))^n)}$, $\frac{eK_p^n}{\delta_p(K_p^n+(x(t-\tau))^n)}$ and $\frac{eK_f^n}{\delta_f(K_f^n+(y(t-\tau))^n)}$ respectively indicate that the TetR is inhibited by LacI, CI is inhibited by TetR, and LacI is inhibited by CI, which constitutes a negative feedback loop. These terms are consistent with the original model. Studies have shown that AHL degradation enzymes AiiA and LuxI have a certain promotion effect on LacI in the QS system ^[37]. The promotion of AiiA on LacI is reflected by the term $\frac{o(w(t-\tau))^n}{J^n+(w(t-\tau))^n}$ in the third equation. These two terms represent protein-protein interactions, and we choose the Hill function modeling method similar to the original model. The fourth term $\frac{Q(v(t-\tau))^n}{M^n+(v(t-\tau))^n}$ represents that the promotion of LuxI on LacI. The term $\frac{k p(t)}{N+p(t)}$ represents LacI is promoted by AHL, which is consistent with the original model.

Next, protein levels of AiiA and LuxI are explained. The term $\frac{u(z(t-\tau))^n}{q^n+(z(t-\tau))^n}$ represents that AiiA is promoted by LacI with rate $u$. The term $\frac{d(z(t-\tau))^n}{\delta^n+(z(t-\tau))^n}$ reflects that LuxI is promoted by LacI with rate $d$. However, H$_2$O$_2$ promotes AiiA and LuxI mainly in an unique way. Under normal conditions, ArcA partially inhibits Lux promoter. However, it is inactivated under oxidative conditions triggered by H$_2$O$_2$, alleviating the inhibition of Lux and increasing AiiA and LuxI ^[32]. In addition, the evidence shows that the production of H$_2$O$_2$ on AiiA can be controlled by $f_i(Z) = (\delta_i+Z)/(c_i+Z)$ ^[38]. Based on these chemical reactions, the promotion effects of H$_2$O$_2$ on AiiA and LuxI are represented by two functions, i.e., $\frac{\xi (\delta_i+s(t))}{V+s(t)}$ and $\frac{\mu (\delta_l+s(t))}{H+s(t)}$, with rates $\xi$ and $\mu$, respectively.

Finally, the levels of AHL and H$_2$O$_2$ are explained. The term $-\frac{r w(t-\tau)p(t)}{h+p(t)}$ represents that AHL is hydrolyzed by AiiA. The term $\frac{\beta v(t-\tau)p(t)}{\alpha+p(t)}$ represents that AHL is promoted by LuxI. Study has shown that AiiA and LuxI have a binding process with AHL ^[32]. Therefore, we treat the process as a form of multiplication, which is similar to most kinetic modeling. $l$ refers to the diffusion coefficient of AHL out of the cell membrane. Here, $U$ refers to the extracellular concentration of AHL. In fact, $U = \frac{N}{Q}\underset{i = 1}{\overset{N}{\sum\limits}}p_i$, $N$ is the total cell number, $Q$ is the intensity of communication. However, here we consider the regulation of the above three molecules and delay on the single cell dynamics, so choose the approach similar to ^[25], i.e., $U = 0$. In other words, we mainly study the single cell dynamics and omit the effects of cell-cell communication. Studies have shown that LuxI can promote the production of H$_2$O$_2$, and LuxI also has a binding process with H$_2$O$_2$ ^[32]. Therefore, the promotion effect of LuxI to H$_2$O$_2$ is represented by $\frac{\rho v(t-\tau)s(t)}{\sigma+s(t)}$.

The complete list of parameter descriptions and default values is given in Table 1. We adopt some parameter values within the ranges, which are used in other similar quorum sensing studies ^{[28,32,37,38]}. Other non-reference parameter values are selected within a certain range, which can be referred to ^[39]. The reason for choosing these parameters accurately is to approach the production rate and degradation rate of most proteins in the QS network. In the model, positive regulation, i.e., promotion, and negative regulation including degradation or inhibition are all represented by Hill functions. In addition, here, the time delay $\tau$ of each protein synthesis in this network is considered to be the same. The formation of a protein requires a certain time $\tau$, when another protein or substance interacts with it, the starting time is $t-\tau$, which is similar to the introduction principle of time delay in p38 and p53 network ^[34,40].

3.   Theoretical and numerical results

In the section, we will focus on the influence the protein synthesis time delay and related parameters of AiiA, LuxI, and H$_2$O$_2$ on the QS network. All theoretical analyses and computer simulations are given by using XPPAUT, Mathematica and Matlab. All parameter descriptions and default values are given in Table 1, and are used in all calculations and simulations.

3.1.   Theoretical analysis

From the perspective of dynamical system, changes in protein synthesis delays in most systems usually cause oscillations, which enrich dynamic behavior and increase mathematical complexity. Therefore, we need to clearly know the critical value of protein synthesis delay, which is the beginning of the oscillation behavior of the QS system. To this end, oscillatory dynamics of the system (2.1) is analyzed by using Hopf bifurcation theory. From the viewpoint of biology, only positive equilibria are of interest. It is assumed that system (2.1) has a unique positive equilibrium $E^{*} = (X^{*}, Y^{*}, Z^{*}, W^{*}, V^{*}, P^{*}, S^{*})$, which satisfies the following equation:

From Eq (3.1), replace $Y^{*}$, $W^{*}$, $V^{*}$, $P^{*}$ with the equations for $Z^{*}$ and substitute them into the third equation in Eq (3.1).

Let

where

where

then

where

Therefore, $H(0) > 0$ if and only if $-A_2 \pm \sqrt{A_2^2-4A_1A_3} < 0$. And $lim_{Z^{*}\rightarrow +\infty}H(Z^{*}) = -\infty$.

Thus, there is at least one positive equilibrium point for the system (2.1) under this condition.

Moreover

Let

Therefore, the system (2.1) has a unique positive equilibrium point if and only if $O_0 < 1$.

Let $\bar{x}(t) = x(t-\tau)-X^{*}$, $\bar{y}(t) = y(t-\tau)-Y^{*}$, $\bar{z}(t) = z(t-\tau)-Z^{*}$, $\bar{w}(t) = w(t-\tau)-W^{*}$, $\bar{v}(t) = v(t-\tau)-V^{*}$, $\bar{p}(t) = p(t)-P^{*}$, and $\bar{s}(t) = s(t)-S^{*}$, $\bar{x}$, $\bar{y}$, $\bar{z}$, $\bar{w}$, $\bar{v}$, $\bar{p}$, $\bar{s}$ are still denoted by $x$, $y$, $z$, $w$, $v$, $p$, $s$, respectively. Therefore, the linearized equation of system (2.1) at the positive equilibrium point $E^{*} = (X^{*}, Y^{*}, Z^{*}, W^{*}, V^{*}, P^{*}, S^{*})$ is given:

where

The characteristic equation of the linearized system (3.2) is

where

Then the following exponential polynomial equation is obtained

where

It is known that $E^{*} = (X^{*}, Y^{*}, Z^{*}, W^{*}, V^{*}, P^{*}, S^{*})$ is locally asymptotically stable if and only if all roots of Eq (3.4) have strictly negative real parts. But the equilibrium will lose its stability if a pair of purely imaginary roots appear. To study the roots of Eq (3.4) with the time delay, we consider two cases, i.e., the system with and without the time delay.

Case I : In the absence of the delay, i.e., $\tau = 0$, Eq (3.4) can be simplified into the form

where

By the well-known Routh-Hurwitz criterion, a set of necessary and sufficient conditions for all roots of Eq (3.5) having negative real parts can be expressed as:

The equilibrium point $E^{*} = (X^{*}, Y^{*}, Z^{*}, W^{*}, V^{*}, P^{*}, S^{*})$ is locally asymptotically stable if and only if (H1) $\Delta_{i} > 0 \; (i = 1, \cdots, 7)$ holds. In order to verify whether the above inequality holds under the parameter values of Table 1, we substitute them into the above formulas. We get $\Delta_{1} = 10.2466$, $\Delta_{2} = 43.2115$, $\Delta_{3} = 114.8$, $\Delta_{4} = 140.696$, $\Delta_{5} = 157.22$, $\Delta_{6} = 3.47169$, $\Delta_{7} = 0.026455$. Therefore, when $\tau = 0$, the system (2.1) is asymptotically stable under the parameters in Table 1. The theoretical analysis is consistent with the numerical simulation in Figure 2.

As mentioned earlier, the quorum-sensing oscillatory behavior can start the bacterial population to lyse the drug and release it. This means that when other conditions are fixed and all protein synthesis does not require a specific time, the phenomenon of bacteria carrying drugs to the target cannot be initiated. In addition, the protein synthesis time is regarded as instantaneous, which violates the actual biological meaning.

Case II : In the presence of the delay with $\tau > 0$, Eq (3.4) has the form

It is assumed that $\pm i\omega \; (\omega > 0)$ are a pair of pure imaginary roots of Eq (3.6), which implies that $\omega$ must satisfy the following equation

Separating real and imaginary parts, we get

where

By simple calculation, the following equations are obtained

where

Employing $\cos^{2}(\omega\tau)+\sin^{2}(\omega\tau) = 1$, one can obtain the following equation

where

If Eq (3.10) has a positive real root $\omega_{0}$, the Eq (3.6) has an imaginary root $\lambda = i \omega_{0}$. Therefore, then we obtain the following critical value of the time delay

At $\tau = \tau_{0}$, the characteristic equation (3.7) has a pair of purely imaginary roots $\pm i\omega_{0}$. Obviously, the system (2.1) undergoes a Hopf bifurcation at $\tau = \tau_{0}$. The critical value $\tau_{0}$ is essential to determine if an oscillation can occur. Under the given parameter values in Table 1, we can figure out the critical value $\tau_0 = 4.4378$.

By taking derivative with respect to $\tau$ in Eq (3.4), we obtain

where

and

Then we can easily find

where

Obviously, according to the data in Table 1, Eqs (3.10) and (3.11), we obtain $GM > 0$ and $VK > 0$, and thus

These results demonstrate that the roots of characteristic Eq (3.4) cross the imaginary axis at $\lambda(\tau_{0}) = \pm i\omega_{0}$ from left to right. Clearly, when $\tau = \tau_{0}$, the system (2.1) loses its stability and undergoes a supercritical Hopf bifurcation at the equilibrium point $E^{*} = (X^{*}, Y^{*}, Z^{*}, W^{*}, V^{*}, P^{*}, S^{*})$, according to Hopf bifurcation theorem. This indicates that when the protein synthesis time reaches a certain critical value, the bacterial population will initiate lysis and release the drug. In addition, the quantitative control of drugs will be given in the next section. How to actually regulate protein synthesis time is also considered below.

3.2.   Numerical results

3.2.1.   The Individual perturbation of time delay

To determine the important mechanism of time delay in the QS system, it is crucial to detect the dynamic behavior without time delay. Figure 2 shows the dynamics of QS network without the delay, i.e., $\tau = 0$, evolves toward an equilibrium and asymptotically stable. Conversely, when the time delay is introduced into the QS circuit, the critical value of the time delay that induces the oscillation of the QS system can be calculated, i.e., $\tau_{0}$ = 4.4378 min. When $\tau < \tau_{0}$, the system still converges to an equilibrium and asymptotically stable, as shown in Figure 3. As $\tau$ increases and passes through the critical value $\tau_{0}$, the system will reach a regime of sustained oscillations. The corresponding periodic oscillation at $\tau$ = 5.5 min is shown in Figure 4. The results show that the QS system changes from asymptotic stability to continuous periodic oscillation with the increase of delay $\tau$.

To more specifically study the effectiveness of the delay on oscillations, we next further study how the variation of $\tau$ affects the amplitudes and periods of the oscillations. The time evolution diagrams of various time delays are presented in Figure 5. Obviously, with the increase of the $\tau$, the amplitudes and periods of the oscillations increase significantly. In general, these results suggest that the amplitudes and periods of these oscillations are sensitive to the variation of $\tau$. In other words, the time delay leads to Hopf bifurcation, which drives the QS oscillation and controls its amplitude and period. These results indicate that the time delay is necessary for the oscillation behavior of the QS system, and the critical value of the time delay is a prerequisite for oscillation. As noted above, the amplitude and period of the QS oscillations are crucial in terms of the dose and time of the bacterial synthesis and release of drugs ^[18,19]. These results indicate that perhaps the amount and timing of drug release can be controlled by individually regulating the delay time in the QS network.

3.2.2.   The influence of parameters related to AiiA, LuxI and H$_2$O$_2$ on QS dynamics

The interaction principle of AiiA, LuxI, H$_2$O$_2$ to the QS system has been widely involved. However, the impact of the changes of these three substances on the dynamics and actual biology of the QS system is still an unsolved mystery. Therefore, in this part, the effects of parameters related to AiiA, LuxI, and H$_2$O$_2$ on dynamics of system (2.1) is studied by perturbing single parameter and fixing all other parameters.

First, the impact of the parameters related to protein AiiA, i.e., $o$, $r$, and $u$, on the dynamics behaviors is studied. The bifurcation diagram of the system (2.1) regard the parameter $o$ as a control parameter is shown in Figure 6(a). It can exhibit four distinct dynamical regimes with the variation of $o$. When AiiA has a low induction rate to LacI, a subcritical Hopf bifurcation occurs, i.e., $o_{HB}$ = 0.3545. When the induction rate of AiiA to LacI is low, i.e., $o < o_{HB}$, the system maintains a continuous period oscillation, and the amplitudes of these oscillations decrease with the increasing of $o$. Studies have shown that the amplitude and period of oscillations in the QS system determine the dose and time interval of the bacteria released drugs. Therefore, perhaps the dose and time interval for the bacteria to release the drug can be reduced by increasing the production rate of AiiA to LacI. When the induction rate of LacI by AiiA reaches 0.9858, i.e., $o_{HB} < o < o_{LP2}$ = 0.9858, the system is monostable with a higher TetR level. When the induction of LacI by AiiA is at an intermediate level, $o_{LP2} < o < o_{LP1} = 1.425$, the system becomes bistable. When the induction of LacI by AiiA is at a high level, i.e., $o > o_{LP1}$, the system is also monostable with a lower TetR level. This means that as the induction of LacI is increased by AiiA, the QS system will switch from oscillation to bistable. Figure 6(b), (c) shows the similar dynamical properties of the other two parameters $r$, $u$ related to AiiA, which are opposite to the parameter $o$. Particularly, $r$ and $u$ occur supercritical Hopf bifurcation occurs at $r_{HB}, u_{HB}$ and two saddle-node bifurcations occur at $r_{LP2}, u_{LP2}$ and $r_{LP1}, u_{LP1}$. When $r, u < r_{LP1}, u_{LP1}$ or $r_{LP2}, u_{LP2} < r, u < r_{HB}, u_{HB}$, the system is monostable. The system becomes bistable at the region $r_{LP1}, u_{LP1} < r, u < r_{LP2}, u_{LP2}$. While oscillations occur when $r, u > r_{HB}, u_{HB}$. Furthermore, the amplitudes of these oscillations increase with the increasing of $r$, $u$ in a wide range. Therefore, maybe the dose and time interval of the drug released by the bacteria in QS can be increased by increasing $r$ or $u$. These dynamics and applications are exactly opposite to effects of increasing $o$.

Second, the regulation of parameters related to LuxI, i.e., $Q$, $d$, $\beta$, on the on dynamics behavior is analyzed. From the results of Figure 7(a)–(c), it can be seen that the induction rate of LuxI to LacI $Q$, the induction rate of LacI to LuxI $d$ and the induction rate of LuxI to AiiA $\beta$ have similar kinetic characteristics. As $Q$, $d$, $\beta$ increases, the level of TetR shows a downward trend. Specifically, as the three parameters increase, TetR initially undergoes a Hopf bifurcation and maintains a continuous oscillation. However, both the amplitude and period of the oscillation are reduced. Therefore, the dose and time interval of the bacteria releasing the drug can be reduced by controlling the production rate of LacI induced by LuxI, LuxI induced by LacI, or AHL induced by LuxI. As they continue to increase, TetR jumps from a high-level monostable state to LP2, then continues to increase and jumps to LP1, which reflects bistable switch, and finally returns to a low-level monostable state. Therefore, perhaps changing the parameters related to LuxI can also induce changes in the PH of plants and microorganisms.

Finally, the bifurcation diagrams that regard H$_2$O$_2$ related parameters $\xi$ and $\mu$ as control parameters are shown in Figure 8(a), (b), respectively. In Figure 8(a), the entire evolution process goes through from the low-level stationary state to the bistable to the high-level stationary state, and finally produces continuous oscillations. However, it can be seen from Figure 8(b) that the induction rate of LuxI induced by H$_2$O$_2$ is similar to the qualitative properties of several other parameters. In this case, when LuxI is adjusted by H$_2$O$_2$ appropriately, the system is in an oscillating state. Therefore, the production rate of LuxI induced by H$_2$O$_2 \; (\mu)$ and other parameters have similar roles in the release of drugs from bacteria and the regulation of PH in plants and microorganisms.

3.2.3.   The combinatorial interference of parameters related to AiiA, LuxI and H$_2$O$_2$ on QS dynamics

In this section, two-parameter bifurcation diagrams are drawn to address the dynamics mechanism of the interaction of AiiA, LuxI, H$_2$O$_2$ by simultaneously perturbing two parameters. Regions of different dynamical behaviors in $u$ and $o$ parameter plane is shown in Figure 9(a), which indicates the effects of combinatorial perturbations of AiiA and LacI on system (2.1). The plane is divided into three regions at which monostability, bistability, and oscillations may occur. It can be seen that when the induction rate of LacI to AiiA, $u$, and the induction rate of AiiA to LacI, $o$, are high, the system maintains an oscillating state. When the induction rate of LacI to AiiA, $u$, are low, and induction rate of AiiA to LacI, $o$, increases appropriately, the system is in bistable state. However, when the induction rate of LacI to AiiA, $u$, and the induction rate of AiiA to LacI, $o$, are higher level or lower level, and the system is transformed into monostable state. Especially, different stable states may coexist, i.e., coexistence of equilibrium and periodic oscillations, the state to which the system evolves depends on initial conditions. As shown in Figure 9(b), when the induction rate of LuxI to LacI, $d$, and the induction rate of LacI to LuxI, $Q$, are within a certain range, the system appears a coexistence of monostability and oscillation. The main reason for this phenomenon is that when $d$ and $Q$ are in (Monostability, oscillation), the LP2 point is located on the left side of the HB point, so there is a coexistence of the oscillation and monostable. Similar dynamic simulations appear, as shown in Figure 10(a), (b), which shows the combinatorial effects of LuxI and H$_2$O$_2$, AiiA and H$_2$O$_2$, respectively. Figure 11(a), (b) both show the combinatorial effect of LuxI and LacI, LacI and AiiA, respectively.

These results show that the simultaneous perturbation of any two parameters related to AiiA, LuxI, or H$_2$O$_2$ can induce occurrence of oscillations and bistability, and even the coexistence of monostability and oscillation. In short, the combinatorial perturbations of AiiA, LuxI and H$_2$O$_2$ are needed to control the dynamics of the QS system. As mentioned earlier, theoretical studies have shown that the dynamical behaviors induced of the QS system are essential for controlling the PH value of microorganisms and plants and bacterial delivery of drugs. Therefore, this means that the simultaneous changes of important molecules in the QS system may be used to control these biological functions.

3.2.4.   The influence of parameters related to AiiA, LuxI, H$_2$O$_2$ on the critical value of delay

The effects of combinatorial regulation of AiiA, LuxI, H$_2$O$_2$, and the delay are reflected by the sensibility of the critical value $\tau_{0}$ to the parameters. It is clearly elucidated that the delay is critical to induce oscillations, which indicates that the critical value $\tau_{0}$ is one of prerequisites for the oscillations. In addition, the size of the time delay has great flexibility in adjusting the amplitude and period of these oscillations. However, the sensibility of the critical value $\tau_{0}$ to the parameters related to AiiA, LuxI, H$_2$O$_2$ is not fully understood, which is reflected by the effect of combinatorial perturbations. Next, the influence of important parameters related to AiiA, LuxI, H$_2$O$_2$ on the critical value $\tau_{0}$ is studied.

First, how the critical value $\tau_{0}$ depends on the parameters $o$, $r$, $u$ related to AiiA has been carefully probed. The relationship between $\tau_{0}$ and $o$ is shown in Figure 12(a). It can be observed that with slightly increase of $o$, the trend and degree of the critical value $\tau_{0}$ keeps increasing and the rate rises rapidly, which indicates that the critical value $\tau_{0}$ is sensitive to the variation of $o$. Similarly, Figure 12(b), (c) show the sensitivity of the critical value $\tau_{0}$ to $r$ and $u$. It can be seen that as $r$ and $u$ increase, the critical value $\tau_{0}$ value continues to decrease, but the change is extremely imperceptible, which indicates that $\tau_{0}$ is insensitive to the increase of $r$ and $u$. Next, the sensitivity of $\tau_{0}$ to the parameters $Q$, $d$, $\beta$ related to LuxI is analyzed. The Figure 13(a) shows the sensitivity analysis of $Q$ and $\tau_{0}$, which is similar to the relationship between $o$ and $\tau_{0}$, which shows that $\tau_{0}$ is sensitive to $Q$. Instead, $\tau_{0}$ is insensitive to the variation of $d$, $\beta$ parameters, as shown in Figure 13(b), (c). Finally, the sensitivity of $\tau_{0}$ to the parameters $\xi$ and $\mu$ related to H$_2$O$_2$ is considered. The Figure 14(a) shows that $\tau_{0}$ monotonically decreases and the speed is very slow with the variation of the parameters $\xi$. While with variation of $\mu$, $\tau_{0}$ also increases but the rate is also very slow, as shown in Figure 14(b). In summary, these results means that the critical value $\tau_{0}$ is insensitive to $\xi$, $\mu$ and $\beta$. The critical value $\tau_{0}$ is the most sensitive to $o$ and $Q$.

As mentioned earlier, the threshold $\tau_{0}$ is a critical condition for initiating oscillations. When the time delay $\tau$ exceeds the critical value $\tau_{0}$, a larger delay induces a larger amplitude and period. In addition, when $\tau$ is less than the critical value $\tau_{0}$, we can change the value of $o$ and $Q$ so as to compensate the inefficiency of $\tau$ and induce the occurrence of oscillations. In other words, by fine-tuning the most sensitive parameters, the QS system can quickly realize the transition from steady state to oscillating state.

4.   Discussion and conclusions

In this paper, a computational model is constructed by incorporating existing biological knowledge and biochemical reaction modeling methods, with the aim of capturing the dynamics mechanism of QS system. To this end, first, by considering the delay required for protein synthesis as a bifurcation parameter, its individual effect on the oscillatory behavior of the QS system is investigated. The research results find that $\tau_{0}$ is the critical indicator to start the oscillation of the QS system. In addition, the amplitudes and periods are sensitive to the slight variation of the delay. Therefore, this reveals that the size of the time delay is an advantageous condition for controlling the period and amplitude of QS circuit. Second, the influence of combinatorial perturbation of AiiA, LuxI and H$_2$O$_2$ on the oscillatory behavior of QS system is studied. It is indicated that the dynamics region of the QS system is generated by simultaneously perturbing any two parameters related to AiiA, LuxI, H$_2$O$_2$. Third, the dependence of the dynamic properties of QS system on AiiA, LuxI and H$_2$O$_2$ is analyzed. The result indicates that concentration change of AiiA, LuxI and H$_2$O$_2$ can induce the bistability and oscillatory behavior of QS system. Finally, the influence of the combinatorial perturbation of AiiA, LuxI, H$_2$O$_2$ and delay on the QS system is analyzed. It is explained by the sensitivity analysis of the critical value $\tau_{0}$ to the parameters related to AiiA, LuxI, and H$_2$O$_2$. The results show that when the induction rate of AiiA to LacI, $o$, and the induction rate of LuxI to LacI, $Q$, are increased, the change of critical value $(\tau_0)$ of delay is extremely obvious. Therefore, the critical value $\tau_{0}$ is the most sensitive to the parameters $o$ and $Q$.

Extensive experimental research conclusions indicate that the dynamical behaviors of the QS system will have an important impact on certain bacterial populations. For instance, the bistable switching of the QS system in plants and microorganisms can be used to control the conversion of PH value ^[17]. In addition, the fine adjustment of the period and amplitude of the QS system also has important applications in the dosage and time of the scattered medicine ^[20]. However, the amplitude and period of these stability and oscillations may be interrupted through minor change of certain components in the QS system. Interestingly, we find that changing the length of the time delay and adjusting the concentration of AiiA, LuxI, and H$_2$O$_2$ have become important issues to provide verifiable predictions for experiments. In addition, how to interfere with the protein synthesis time and these important substances through actual external means are also considered. The length of time delay may be changed by transcription inhibitors (Rif) and translation inhibitors (Cam, Ksg) ^{[41,42,43,44]}. Moreover, AiiA and LuxI can be cloned from certain bacteria by PCR ^[45,46,47]. Finally, H$_2$O$_2$ in the QS system can be produced by adjusting the concentration of SOD ^[48]. Therefore, focus on the dynamic characteristics of the QS system and its important components will help to further understand some mechanisms of QS, and provide better ideas for practical biological applications.

However, our study also has some restrictive assumptions and unresolved issues. The choice of model parameters is always a top priority for dynamics. Therefore, it is more interesting to determine the kinetic properties through statistical random parameter selection within a certain range. In addition, although the cell density has an important effect on the cell fate, due to the complexity of the regulatory network, we only consider the single cell dynamics which is the first step toward the multicellular dynamics. It would be an interesting topic to further analyze the dynamical effects induced by cell density. Finally, in the QS system, the diffusion of the signal molecule AHL is very important. In future, the diffusion term of AHL can be considered in partial differential equations, and the dynamical effects of time delay and diffusion terms on the system can also be analyzed by referring to the modeling and theoretical methods of Li et al. ^[49]. Additionally, it is expected that the spatial diffusion kinetics of these AHLs can guide rational therapeutic strategies and improve agricultural practices ^[50,51].

Acknowledgement

This research is supported by the National Natural Science Foundation of China (Grant No. 11971297).

Conflict of interest

All authors declare no conflicts of interest in this paper