Rule–based regulatory and metabolic model for Quorum sensing in P. aeruginosa

Computational analysis

In this multi–level formalism each node can adopt several possible states (see Additional file 3: Table S3). An overview of the updating scheme and work flow is given in Figure 3. We consider a certain number of cells that share the same environment, which means that the external autoinducer concentrations are equal to each cell. Usually, the simulation starts with a single cell. In the results section, we only discuss the results for this first cell. During exponential growth, i.e. before completing six cell divisions, a cell divides into two cells after every 60 time steps. The simulations are iterated until time step 600. Each new cell is initialized in the same way as the first cell and has a delay time of ten steps during which only transport processes are possible. Thus, in this model, the first cell represents a biological cell at the end of the exponential growth phase.

whereby W=(w _{i j} ) is a weighting matrix and contains all relations between the nodes. Its entry w _{i j} is $\frac{1}{l}$ when node j is one of l nodes that are together required to activate node i. If node j deactivates node i then w _{i j} is −k where k is the number of reactions that activate node i. Otherwise, w _{i j} is zero. In contrast to some other Boolean models, the second case implies that one inhibiting reaction overrules the activating reactions of node i. Also, a node that is not explicitly activated at a certain time step is inactive. As our network only included deactivating nodes in cases where there is a good evidence that they inhibit certain activations, the weighting scheme makes sure that deactivations are stronger than activations. As the concentration of the autoinducers must reach a certain threshold ε to start the positive feedback loop, some reactions in our model are only activated if a critical threshold of certain metabolites is exceeded. Function K controls that the state of every node j that is related to node i has reached the required threshold at time t. Function H rounds the weighted sum to an integer valued output. For multi–level nodes with more than two possible states, i.e. autoinducers, RsaL, and virulence factors the level at the current time step x _i (t) is added via function G. Function M ensures that no node adopts a higher value than its maximum possible value. The explicit forms of the functions M, H, K, and G are given in Additional file 4.

After updating all nodes, the states of multi–level nodes are decreased by the number of reactions in which they were used. Let us assume, that for example, HHQ forms a complex with PqsR and HHQ is additionally required to yield PQS. Then, in this step the level of HHQ is decreased by two. Those multi–level nodes are subject to a further reaction that decreases the state value after every 20 time steps to account for degradation (see Section Choice of parameters). This degradation is realized equally for all multi–level nodes independent of their maximum possible states. A time step represents either the whole formation of a protein including transcription and translation or a faster enzymatic reaction. Hence, a degradation happens infrequently. However, a node has a relatively small number of possible states, such that a single degradation process consumes a comparably large part of its concentration. For example, in the case of RsaL a degradation decreases the level by two–thirds. Transport reactions (dotted lines in Figure 2) are neither considered in the logical updating scheme nor in the decreasing step afterwards. Additionally, we included some random reactions for processes that take place either infrequently or very slowly in comparison to the other processes, for example, the transformation of HHQ to PQS [18]. For this, we used a random reaction that happens about every fourth step. For simplicity, the enzymes PqsB, PqsC, and PqsD are considered as one node.

The communication between different cells is based on the exchange of the signaling molecules. For this, we modeled a simple transport process that occurs after updating all states inside the cells. When the internal concentration of an autoinducer is higher than the transport threshold, it is transported out of the cell except for the case where the external concentration relative to the number of cells is higher than the internal concentration. Otherwise, for external levels larger than zero, an autoinducer is transported into the cell. In the case that there is no autoinducer outside, there is also no transport. We do not distinguish between diffusion and efflux pumps or active transports and apply usually a transport threshold of three, i.e. if the autoinducer concentration is high.

Knock–out mutants, enzyme inhibitors and receptor antagonists

Besides the wild type also knock–out mutants and inhibitors were considered. For such mutants (PqsR ⁻, PqsA ⁻–PqsBCD ⁻, PqsE ⁻, PqsE ⁻–PqsR ⁻, PqsE ⁻–PqsA ⁻–PqsBCD ⁻) all respective formation edges were removed. Conversely, inhibitors (PqsR antagonist, PqsBCD inhibitor) block a certain fraction of the target. Therefore, inhibition of the target was set to happen with a certain inhibition level that indicates at what percentage the target is blocked (see Figure 4). An inhibitor is initialized as being active, but the time when the inhibitor blocks its target is chosen randomly.

Choice of parameters

For the sake of simplicity and interpretability, the model was designed as a minimal system. A synchronous, parameter–free Boolean logical analysis was applied where nodes can be either “on” or “off” except for some nodes and edges that are discussed below in detail. In a Boolean network with a synchronous updating scheme, there is no time scale so that all reactions take place at the same time with the same rate. This is a valid approach unless this alters the system behavior in a crucial waycompared to either considering different reaction time scales, or if certain reactions depend on special conditions. We will show in the following that the first assumption (same time scale) is not crucial for the system properties studied here. The second point, however, required a generalization of the Boolean model by introducing multi-level logical functions at various points of our model.

We will now discuss the elements of the systems where more than two states were found necessary. (i) The regulator RsaL blocks the transcription of the autoinducer synthase LasI. Hence, three different RsaL levels were used, namely “not available” (0), “available at low concentrations” (1), and “sufficient to block LasI” (2) are used. (ii) Besides PQS, HHQ is also able to bind to PqsR, and to activate PqsR, but PQS has a 100–fold higher affinity than HHQ [18]. Therefore, PQS was set to inactivate the complex C5 of HHQ and PqsR with a threshold of two representing a large concentration of PQS. (In Figure 2, all reactions with a threshold of two are colored in blue). (iii) Since the activation of gene expression (activation of Ci:Gi by Ci) is concentration–dependent, the complexation between autoinducer and receptor occurs only when the autoinducer concentration has reached the lower threshold of two. (iv) Moreover, the autoinducer of the las system AI–1 is able to form a weak complex with the second receptor RhlR as long as the autoinducer of the rhl system AI–2 has not reached a very high concentration. Thus, formation of the complex C4 between AI–1 and RhlR required at least four AI–1 and its inactivation at least four AI–2 (edges colored in green in Figure 2). (v) Due to the accumulation of autoinducers and virulence factors, six different levels were considered for the internal signaling molecules AI–1, AI–2, HHQ, as well as PQS and the virulence factors. The external autoinducer levels were allowed to adopt arbitrarily large numbers.

Conversion frequency

HHQ is formed by one cell, transported into another cell, and there converted into PQS [17]. Although this may yield to a time delay in reality, our model assumes that both steps take place in a single cell. Therefore, we designed this reaction (colored in gray in Figure 2) such that it occurs randomly with a certain probability. Figure 5 shows how the autoinducer and pyocyanin levels depend on the reaction rate of this process. Due to the activation of pyocyanin biosynthesis by PqsE and the production of PqsE by the complex C5 using HHQ, the pyocyanin level is independent of the reaction rate. Kesarwani et al. reported that the HHQ concentration is about 12% of the PQS concentration in the beginning of the stationary growth phase [38]. To match this experimental finding, we used a conversion frequency of 55% in the following.

Initialization and starting behavior

In this study, we are primarily interested in identifying attractors in which the QS systems are active. Obviously, when all initial states of the network are set to zero (inactive or not available), the whole system will remain inactive. The las system becomes active when (1) one or both of the global activators GacA and Vfr are initially set to one and when (2) the complex C1 between AI-1 and LasR binds to the las operon (C1:G1 is set to one). The rhl system is activated with a certain time delay. This delay can be avoided when GacA, C1:G1, and additionally C2:G2 (complex C2 between AI–2 and RhlR binds to the rhl operon) are active in the beginning. If in the beginning C5:G3 and C3:G3 are activated as well, the pqs system becomes active after a short delay, too. Also, switching on PqsR avoids the delay. In general, the initialization of Cx:Gy (complex Cx bound to operon Gy) to one is necessary since basal transcription rates are not included in the logical approach used here. In our system, the basal transcriptions could be included either with the same rate as the activated gene expression or in a stochastic way as often as degradation. The second case does not result in an initiation of the system. Alternatively, it is possible to start with elevated autoinducer levels or to begin with externally available autoinducers that are able to diffuse or can be transported into the cells.

When more nodes are initially active, the initialization phase of the system takes a different route because it is no longer required to start up the whole system. However, at later times, the behavior and the levels are the same as using minimal initial conditions. Additional file 7: Figure S3 shows that, after about 100 time steps, the metabolite levels reached in a computation with all nodes initially activated are very similar to those reached with minimal initial conditions. For all other further considerations, we used minimal initial conditions in which only Vfr, C1:G1, C3:G3, and C5:G3 are initially set to one.

Additional file 8: Table S4 shows an example trajectory where the system has reached a stable regime after time 100. Therefore, we can consider all time steps after 100 for computing equilibrated average values. For statistical reasons, we carried out all simulations until a time of 600. The differences of the levels using different growth rates were found to be very small. We compared growth processes using different division frequencies at time steps 60, 120, and 240 and modified the simulation time such that in each case about 200 time steps in stationary phase are possible. For most internal nodes, the levels were the same. Due to the stochastically modeled formation of PQS, the levels differ from each other in the same range as between different runs that used the same growth parameters but different random numbers for the pqs system. The external levels also display only negligible variations due to degradations. For example, the external level of autoinducer AI–1 divided by the total number of cells is 3.1 for a growth rate of 60, 3.3 for a growth rate of 120, and 2.9 for a growth rate of 240 averaged over the last 100 time steps. When comparing the behavior of a single cell over 600 time steps to that of 64 cells the internal and external levels are slightly changed. However, the qualitative behavior remains unchanged as well. For a single cell, the external level of AI–1 averaged over the last 100 time steps is 2.6. Hence, the behavior of the model is practically independent of the growth processes. In the following, we used a growth rate of 60 time steps and 6 division processes.

Figure 6 shows the behavior of the complexes at the start of the simulation. In general, we defined a delay time of 10 time steps for each new cell during which only transport processes are possible. At first, the las system becomes active. The first peak is at time 13 because three reactions are required to build complex C1 starting from C1:G1. Those reactions include the whole biosynthesis. After time 19, the las complex C1 is always active. After a short delay, the rhl and pqs systems are activated since C1 is necessary to activate C1:G2 and C1:G3. Three further reactions are required to build the rhl complex C2. However, C2 is not activated directly since the level of AI–2 is too low. Therefore, these reactions have to take place twice. Then, the first peak of C2 is at time 20. The first peak of the C3 and C5 complexes of the pqs system occurred in the time interval between 19 to 23 and at time 16. PqsR and PqsH are activated by C1:G3 at time 15 for the first time. Until time 15 the level of HHQ is two and the level of PQS is one, such that at time 16 the complex C5 between PqsR and HHQ is formed. Afterwards, the activation of C3 and C5 depends on the stochastically modeled formation of HHQ by PqsH. When C3 is active C5 is inactive since PQS inhibits C5.

Behavior of wild type and knock–out mutants

In the following, network initialization is performed as before using the described minimal system. Figure 7 summarizes the behavior of the wild type system after the starting phase. The virulence factors LasB, Rhm2, and pyocyanin are generated at maximum level. The level of AI–2 is higher than the levels of AI–1, HHQ, and PQS since there is no negative regulator in the rhl system whereas there is RsaL in the las system and PqsE in the pqs system. Due to the high level of AI–2, the complex C4 between the first autoinducer AI–1 and the second receptor RhlR is inhibited. Complex C5 that is formed from PqsR and HHQ is clearly less active than the complex C3 between PqsR and PQS since the averaged internal PQS level is almost two.

Now, we focused on the pqs system shown in Figure 8 in order to illustrate the behavior of the network topology in detail. All other parts of the network were still initialized and processed as described before. Additionally, we considered the effects of knock–out mutants and inhibitors of the PQS biosynthesis or antagonists of the receptor PqsR. As mentioned before, the complexation between autoinducer and receptor requires at least an autoinducer level of two in our model (reactions colored in blue). Further, when the level of PQS is two or higher, HHQ can no longer bind to PqsR (inhibition of C5 by PQS colored in blue). HHQ and PQS can only be transported outside of the cell when their level is three or higher (colored in orange). The red colored edges denote degradations that decrease the level of the respective node by one after every 20 time steps. The gray colored reaction happens stochastically with a conversion frequency of 55%.

Figure 9 shows the behavior of the knock–out mutants modeled here, namely PqsR ⁻, PqsA ⁻, PqsBCD ⁻, and PqsE ⁻, as well as the double mutants PqsA ⁻– PqsBCD ⁻, PqsE ⁻–PqsR ⁻ and PqsE ⁻–PqsA ⁻–PqsBCD ⁻ (PqsABCDE ⁻) in comparison to the wild type. Knock–out mutants within the pqsA–E operon and pqsR are known to produce zero or drastically reduced amounts of pyocyanin compared to wild type cells [17, 20, 39]. As reported in the literature deleting PqsR, PqsA, or PqsBCD in our model results in a deficiency in PQS and pyocyanin formation. The single PqsE ⁻ knock–out mutant reached a higher internal PQS level than the wild type due to the missing inhibition of C3:G3 and C5:G3 by PqsE and therefore produced pyocyanin near the maximum level. Yet, it was reported that PqsE ⁻ knock–out mutants produce clearly less pyocyanin than wild type cells whereas an increased concentration of PqsE enhances the pyocyanin concentration [20, 40]. We will return to this issue below where we will investigate the network topology around PqsE in detail. PqsE is able to regulate pyocyanin independent of PQS and PqsR, but subject to the las and rhl system [40]. The external level of PQS in a PqsE ⁻ mutant is strongly increased. In contrast, the external concentration of autoinducers was reported to be similar in the wild type and PqsE ⁻ knock–out mutants [40]. However, in PqsE ⁺ mutants the concentration of HHQ and PQS is reduced [20].

Behavior of enzyme inhibitors and receptor antagonists

Next, we analyzed the effect of blocking different target proteins on the levels of the autoinducers HHQ and PQS as well as the level of the virulence factor pyocyanin. We compared inhibitors of the signal molecule biosynthesis and antagonists of the PqsR receptor and studied the influence of the inhibition level. For example, an inhibition level of 40% means that 40% of all receptors or enzymes of the corresponding type are blocked. Figure 10 shows the external HHQ and PQS levels for the unblocked system (inhibition level of 0%) and for systems including inhibitors.

In the presence of PqsBCD inhibitors, the external HHQ and PQS levels are noticeably decreased with a high dependence of the inhibition level. For example, a PqsBCD inhibitor with inhibition level of 30% reduces the external HHQ level by about 80% and the external PQS level by about 55% (see Figure 10). Indeed, Storz and coworker showed that a PqsD inhibitor reducing the HHQ concentration by 77% also reduced the PQS concentration by 42% [41]. Figure 11 shows the effect of PqsBCD inhibitors with varying inhibition levels on the internal autoinducer and pyocyanin levels. The pyocyanin level is also clearly decreased. For PqsBCD inhibitors with an inhibition level of 30%, it is reduced by about 10%.

The external PQS level is strongly increased by PqsR antagonists due to missing usage for complex C5 (see Figure 10). Figure 12 shows the effect of weak and strong PqsR receptor antagonists on the internal HHQ and PQS levels and on the pyocyanin level. As long as less than 60% of PqsR is blocked, the pyocyanin level is only very slightly decreased. In contrast, Klein et al. presented PqsR antagonists with affinity in the low micromolar range that indeed reduce the pyocyanin concentration (IC₅₀: 87 μ M) [42]. Further, it was reported that PqsR antagonists, derivates of the natural effector HHQ, with K _D values in a low nanomolar range reduce the pyocyanin formation by about 75% at 3 μ M [43].

Updated network

While the behavior of the HHQ and PQS levels mostly fulfills experimental observations, the discrepancy for pyocyanin between predicted and observed behavior suggests that pyocyanin production is co–regulated by further functionally unknown proteins that may or may not be connected to the pqs system. Here, we only considered the connected case. Thus, we proceeded by studying networks with modified connectivities. For simplicity, we did not consider negative regulations or varying conditions such as reaction rates or modifications in the number of states for certain nodes. Instead of the feed–forward link pqs–member → X → pyocyanin we used a direct link from pqs–member → pyocyanin such that we do not have to include new hypothetical proteins. We further assumed that relations are only possible between a protein acting as receptor and a corresponding signaling molecule.

So far, the natural ligands of PqsE are unknown and it was shown that PQS does not bind to PqsE [44]. According to reaction R19, formation of a complex between PQS and PqsE would lead to up–regulation of pyocyanin. The corresponding modified network is illustrated in Figure 13 and denoted as N3. We postulate that the formation of complex C6 does not require a high concentration of PQS. The edge from C6:G4 to PQS represents the degradation of the complex.

In the networks N5, N6, N7, and N8, PqsA ⁻–PqsBCD ⁻ and PqsE ⁻ mutants produce pyocyanin at maximum level due to the postulated reaction R4, see Table 2. According to this or–linkage, PqsR is essential for the biosynthesis of pyocyanin, so that its level is not decreased in PqsA ⁻–PqsBCD ⁻ and PqsE ⁻ mutants. The levels of all other metabolites are the same as for the corresponding network without reaction R4, e.g. N5 behaves as N1. Besides N1 PqsE ⁻ mutants also produce pyocyanin at maximum level in networks N4, N11, and N12. In the case of N2 and N3, where pyocyanin is activated by a complex between PqsE and HHQ or PQS, no pyocyanin is formed by knock–out mutants. However, weak PqsR antagonists again do not decrease the pyocyanin level.

We consider the networks N9 and N10 as being the closest to literature. In our model, PqsBCD inhibitors are more suitable to reduce the HHQ and PQS level and PqsR antagonists reduce the pyocyanin level at higher rate. In general, the levels strongly decrease with increasing inhibition level and the level of HHQ is clearly stronger decreased than the PQS and pyocyanin levels for the same inhibition level. In the case of network N10, a PqsBCD inhibitor strongly reduces the external HHQ level, whereas the effect on the external PQS level is clearly not as emphasized as for comparable inhibitors of the other networks. In reality, reactions R21 and R22 may be represented in two ways. (1) Two complexes up–regulate the transcription of two different enzymes. The first enzyme catalyzes the formation of a metabolite that is used by the second enzyme to form pyocyanin. (2) A first complex up–regulates the transcription of an enzyme that catalyzes the formation of a certain metabolite and a second complex up–regulates a receptor protein that forms a third complex with this metabolite. Then, the third complex up–regulates the transcription of pyocyanin. In both ways, either the first or the second complex contains PqsR and the other PqsE and one of those receptors uses an autoinducer as substrate. The unknown proteins in our model may be directly or indirectly related to the proteins MvaT and MvaU that are required to form pyocyanin [45]. It was reported that the biosynthesis of pyocyanin is a complex pathway including two operons that share high sequence homology. These are responsible for forming a precursor as well as the enzymes PhzM and PhzS which are required to convert this precursor to pyocyanin [46]. Therefore, we suggest that PqsE might be a regulatory protein and that pyocyanin biosynthesis might be regulated by the two regulators PqsR and PqsE as well as an autoinducer as ligand.