Using G -Functions to Investigate the Evolutionary Stability of Bacterial Quorum Sensing

: In ecology, G -functions can be employed to deﬁne a growth function G for a population b , which can then be universally applied to all individuals or groups b i within this population. We can further deﬁne a strategy v i for every group b i . Examples for strategies include diverse behaviour such as number of offspring, habitat choice, and time of nesting for birds. In this work, we employ G -functions to investigate the evolutionary stability of the bacterial cooperation process known as quorum sensing. We employ the G -function ansatz to model both the population dynamics and the resulting evolutionary pressure in order to ﬁnd evolutionary stable states. This results in a semi-linear parabolic system of equations, where cost and beneﬁt are taken into account separately. Depending on different biological assumptions, we analyse a variety of typical model functions. These translate into different long-term scenarios for different functional responses, ranging from single-strategy states to coexistence. As a special feature, we distinguish between the production of public goods, available for all subpopulations, and private goods, from which only the producers can beneﬁt.


Introduction
Since its discovery in 1977 [1], quorum sensing (QS) has received increasing attention as a mechanism of bacterial communication. Nowadays, we know that QS coordinates a range of behaviours at the population level [2]. One of those is the production of iron-scavenging molecules, so-called siderophores, in Pseudomonas aeruginosa. The bacteria secrete the siderophores and other virulence factors into the surrounding medium, where their benefit can be shared by all cells in the local population, leading to the term public goods [3,4], e.g., with the focus on kin selection in [5] or the conflict between individual and group interests in [6]. However, such a cooperative system is vulnerable to exploitation by non-cooperative cheaters. Gaining all the benefits of QS without incurring the metabolic costs, these cheater mutants have a growth advantage compared to the wild-type. Studies have shown that such mutants arise in laboratory conditions as well as in infection models [7][8][9][10][11][12][13] and act as social cheats. A few examples of these social cheats can be found in [3,14,15] (all considering the Gram-negative P. aeruginosa), while Pollitt et al. [16] provided an example from the Gram-positive S. aereus.
This gives rise to questions about the evolutionary stability of QS. Once cheaters arise (e.g., through loss-of-function mutation), they should theoretically outgrow the producers as they have more resources available to invest in cell division. Cheaters have been shown to outcompete producers both in vitro and in vivo, but QS seems to be evolutionary stable in natural systems nevertheless. In contrast to the aforementioned public goods, private goods are only accessible to the producing cell itself. They are hence innately protected from cheaters and provide their benefit exclusively for cells with a functioning QS system. Apart from extracellular molecules, QS also controls the production of proteins which act within the cell. In P. aeruginosa, one such protein is nucleoside hydrolase (Nuh). As Nuh is involved in metabolising adenosine, only bacteria with intact signal receptors can digest this carbon source. In this way, cooperation via QS provides a private fitness benefit to cooperating cells if adenosine is available as carbon source [17]. Schuster et al. [18] discussed the benefits attained from private goods. Specifically, on the question why public goods are ultimately stable, they suggested that, in fact, the private good is, directly or indirectly, involved in cooperative behaviour.
Extracellular molecules do not provide benefit indiscriminately, but are limited by diffusion and habitat structure. Kümmerli et al. [19] found a negative correlation between habitat structure and water solubility of siderophores, a class of secreted enzymes under QS control in a wide range of bacteria. For highly structured environments such as animal tissues, water solubility of siderophores is high, while microstructures in the environment naturally limit the resulting diffusion. Conversely, water solubility of siderophores is low in unstructured environments such as water habitats. This leads to siderophores clinging to each other as well as to lipid membranes. In this way, a fraction of the siderophores stay with their producer and provide some private benefit.
Several mechanisms, such as kin selection [20] and policing [21], have been described that could explain the evolutionary stability of cooperation and QS despite the advantages cheaters have in such a system.
In addition to experimental approaches, there have been many different mathematical modelling approaches to study QS, from stochastic cellular automata [22] to classic individual-based models [23], including differential equations [24]. In this paper, we develop a mathematical model to investigate the evolutionary stability of QS, which results in a semi-linear parabolic system of equation, following Frank [25] and Mund et al. [26], who used an age-dependent model as basis to check evolutionary stability. We explore the population dynamics between secretors and cheaters and are also interested in the evolutionary pressure that arises through those population dynamics. As such, we employ the G-function ansatz that was introduced by Cohen et al. [27] and further developed since then [28][29][30][31]. With the G-function ansatz, one finds a function G which describes the growth rate of a subpopulation. This function is dependent on the strategy of this subpopulation and both population numbers and strategies of other existing subpopulations (see Figure 1). A strategy can be taken to be any behaviour of interest which can be assigned a real number.
In our case, we take the QS-intensity of the subpopulation as the strategy. A strategy value of zero therefore denotes a QS cheater, while we assume a value of one to be the normal output of QS active secretors. If we take v i as the strategy of the subpopulation, v is the vector of strategy values for all subpopulations. Accordingly, b i and b represent the population numbers of subpopulation i and the complete vector, respectively. The population dynamics for subpopulation i thus follows the equatioṅ At the same time, the per-capita growth function G also governs the evolutionary dynamics of the population. If we calculate the derivative of G with respect to the subpopulation strategy, we find the direction of higher per-capita growth. This is the direction evolutionary forces will drive the population, leading tov (2) To simplify notation from here on, we define Figure 1. Schematic representation of influences on the growth rate of a population: the growth rate G of a subpopulation will be mainly influenced by three factors: its own strategy v i , the environment defined by the complete strategy set v and the whole population b.
This framework allows us to find strategy values which are resistant against others in such a way that, under the influence of natural selection, no other strategy can invade the population. These strategies are called evolutionary stable strategies (ESS) [32].
Many experiments have shown the importance of spatial arrangements when considering bacterial interaction [33,34]. As such, we aim to include spatial variables in our model. This can be done in various ways, but we concentrate on the one that was introduced by Kronik and Cohen [31]. This leads to a system of partial differential equations.
In [35], we showed the existence and uniqueness of weak and classical solutions to the associated semi-linear parabolic system in Equations (1) and (2) with Robin boundary conditions using the coupled upper-lower solution approach. In this paper, we analyse a variety of typical model functions. These translate into different long-term scenarios for different functional responses, ranging from single-strategy states to coexistence. In Section 2, we introduce the basic model and some assumptions we work with, before considering private goods models in Section 3. These models are then analysed numerically in Section 4. In Section 5, we summarise and interpret our findings.

The Basic Model: Using the G-Function Approach
In this case G takes the form: For a formulation with carrying capacity instead of population-dependent death rate, we can transform this equation into which gives us B(v, b)C(v i )/µ as the capacity for the population. We can now look at the derivative of G with respect to v i .
We make some regularity assumptions on G and its derivative to simplify the mathematical analysis. These do not restrict our biological application, as our model growth functions do not include discontinuous behaviour and are smooth.
where T denotes transpose. Additionally, we make some assumptions on G from the biological background.
(III) G(v i , v, b) has a (direct or indirect) negative feedback loop in b i . The population number cannot go towards infinity.
From a certain threshold onwards, higher strategy values no longer not a higher net growth.
Reducing the strategy further than the value u provides no higher growth rate.
These assumptions ensure that Equation (21) does not exhibit divergent behaviour (v → ∞ and/or b → ∞), which is not biologically plausible. If, e.g., It follows that v i ≥ c + λt and thus v i → ∞. This would mean ever-increasing strategy values without some kind of trade-off, which we do not find in nature.
In the special case of QS, we take u to be 0, requiring ∂ 1 G(0, v, b) = 0. This keeps v from leaving the biologically meaningful parameter range R + 0 (production cannot be lower than 0).

Cost and Benefit
To model the G-function for QS, one of the avenues we can take is to divide the impact of v i , v, and b into two parts: a growth term influenced by v i and a benefit provided by v, b, and possibly also v i . This reflects the fact that production of QS molecules is costly to the individual, while the resulting factors are public goods and therefore provide benefit to all bacteria. The additional dependence on v i can be seen as a form of private benefit and is discussed in detail when it occurs.
One important thing to note is that the growth term is actually reduced with rising v i , as increased public good production incurs increased metabolic costs. In this way, less energy is retained for reproduction. We denote the cost term by C : R → R + , and make the following assumptions: (VI) As public good production is costly, C is strictly monotonically decreasing for v i ∈ [0, ∞). (VII) When producing public goods, the growth rate is reduced by a certain factor: While Item (VI) is clear from our assumptions on QS, Item (VII) is not as immediately clear. We introduce Item (VII) because we use this factor multiplicatively for G. Thus, a value of 1 would signify unimpeded growth, while a value between 0 and 1 reduces growth. In this way, we assume that QS costs alone do not lead to negative growth rates.
The benefit of QS is provided by secreting extracellular proteins. We denote it by B : R m × R m → R + and make two main assumptions here: (VIII) There is a limit to how much benefit can be obtained, (IX) There is no benefit if no public goods are produced, Equation (8) models a saturation behaviour for the benefit-even if the cells were producing an infinite amount of extracellular protein, the benefit that can be derived is still capped through saturation of enzymes or similar phenomena. Equation (9) ensures that there is no benefit from QS when there are no living bacteria, or all of them have stopped participating in QS.

Analysis
We know from Assumption (VI) that C(v i ) is strictly monotonically decreasing for positive v i . Combined with Assumption (V), we have C (v i ) < 0 ∀v i > 0, C (0) = 0 and there exists only v i = 0 as an ESS. As the only stable stationary point, v i = 0 will attract all other solutions, leading to the decline of QS.

Models with Private Benefit
Following up on the private goods-hypothesis explained in Section 1, we now assume that there is a private benefit associated with producing the public goods, e.g. a small percentage of the produced enzymes may cling to the producing bacteria. We model this by adding a term B(v i ), which is a benefit term that is solely dependent on the strategy of the subpopulation itself. This B(v i ) should fulfil similar assumptions to those we made of B(v, b), which is why we choose to denote it by the same letter. Overall, we can write for G: It follows that hence for an ESS it must hold that To proceed with a more in-depth analysis, we choose a specific cost function which satisfies all our assumptions, and look at two specific types of model terms for B(v i ).

First Type of Model Terms with Monotonicity Property
The first type of model terms consists of those functions B(v i ), for which B (v i )/v i is monotonically decreasing in v i . This is the case for all possible concave functions (e.g., those with B(0) = 0, B(∞) =B > 0), as well as for Equation (13) allows us to write Equation (12) as forv i = 0. Since we know the left-hand side of this equation to be monotonically increasing, a monotonically decreasing right-hand side means that there can only be one solution to Equation (14). Whether or not a solution exists at all can be determined by looking at the limit behaviour.
We also take a look at the limiting values at 0:

Second Type of Model Terms: Hill Equations
For a second type of model terms, we consider Hill equations. These relatively simple functions are often used in mathematical modelling for cooperative binding of molecules. In the case of QS, Dockery and Keener [36] employed them in the modelling of QS signal regulation. We look at terms of the form with h ≥ 3, as the cases h = 1 and h = 2 fall into the first type of model terms. Then, we have Instead of being monotonically decreasing as before, R(v i ) has the following properties: For a better understanding of the ideas behind, we assume for a moment that , two positive stationary solutions v i to Equation (12) exist.
In that way, which unfortunately cannot be simplified much further, without assumptions on the parameters K and a. We remind ourselves that K is the cost of cooperation, while a denotes the amount of signalling necessary to gain half of the maximal (private) benefits. Equation (19) thus gives upper limits to these terms, dependent on the steepness of the activation curve h as well as the public benefit B(v, b). A larger public benefit is actually detrimental to the existence of positive stationary strategies v i . At the same time, there cannot be more than two intersections in [0, +∞) by Rolle's theorem (see, e.g., [37], p. 134), since where the term in parenthesis is strictly monotonically increasing and the equation therefore only has one ξ ∈ (0, +∞) for which (2KB In biological terms, our results indicate that the surrounding subpopulations have a profound impact on the evolution of cooperativity. If one of the subpopulations is cooperating (B(v, b) = 0), the others will experience bistable behaviour-they might cooperate themselves or not participate in QS, depending on starting strategy. If, on the other hand, all of the subpopulations do not cooperate at first (B(v, b) = 0), then the evolutionary pressure will drive them towards the positive ESS, meaning they pick up QS with time.

Numerical Simulations
The full set of equations we are working with therefore readṡ noting thatḃ i (x, t) = ∂b i ∂t . We used our own experimental data from [34] to derive values for the parameters in our model, which can be found in Table A1. Our aim here was to obtain rough estimates which replicate the general qualitative behaviour, as opposed to fit the experimental data perfectly. To solve the system numerically, we employed the method of lines [38] and used an explicit Runge-Kutta solver on the discretised system. The simulations were run in Matlab.
We assumed a closed system in which the existing subpopulations are mixed but not completely homogeneous in the beginning. This leads to Neumann boundary conditions and the initial distribution displayed in Figure 2. Visual representation of initial conditions for numerical simulations. For initial conditions, we used patches. The two subpopulations (dark and light patches) were mixed but not homogeneously distributed.
We also added in a feature of QS that we have neglected before: the costs of QS depend on the signal level, as QS products are only produced when the signal level is high enough. We can achieve this by modifying the cost term from Equation (13) to as B(v, b) is a measure for the level of QS. We have not done so before, since it does not change the analysis qualitatively, but leads to unnecessarily long and obfuscating terms.

Basic Model
We have seen in our calculations that v i = 0 is the only stable strategy in this situation. We should therefore have both v i → 0 and b i → 0. However, when looking at the simulation results for 3000 h as displayed in Figure 3, we notice that this convergence is very slow indeed. A quick calculation of ∂ 1 G and G for our parameter values confirms that both are close to zero. Thus, even though QS is theoretically unstable, both cheaters and producers persist alongside each other for a very long time, albeit at different densities.

Models with Private Benefit
In Section 3 we distinguish between two types of model functions. Depending on the Hill coefficient, Hill-terms can fall into either one. We investigated the long-term behaviour for h = 1, 2 or 3.
As predicted, the behaviour in these cases is quite different. For h ≤ 2, B(v i ) falls into the first type of model functions. For this type, we postulate that there can only be one positive stationary point for v and if it is stable, and then the zero solution is unstable. This is the case for the parameters we calculated. Additionally, we found that for h = 1 the zero solution is not a stationary point. Thus, both strategies converge towards the stable positive equilibrium, one from above and the other from below. Population 2, which started out as cheaters, gains QS functionality, while there is reduced production from Population 1. Since reducing production is slower to reach the stable point in this parameter constellation, Population 1 succumbs to the population pressure from Population 2 and dies out (see Figure 4). This happens on a rather short time frame of less than 200 h.  Figure 4. Evolution of two populations with different start strategies, using a G-function with private benefit and hill factors of h = 1 (line), 2 (dashed) and 3 (dash-dotted). For h = 1, Population 1 (dark colour) dies out rather quickly, while Population 2 (light colour) gains the QS functionality, albeit at a low level, and remains at a stable population level. If h = 2, Population 1 reduces its QS strategy to a lower, but stable value. Population 2 is unable to compete and dies out. For h = 3, Populations 1 and 2 coexistt in a stable way with similar numbers. One population is QS active while the other is not.
For h = 2, we can see in Figure 4 that once again the strategies converge to a positive value. However, this time the zero solution is a stationary point, although an unstable one. Hence, Population 2 cannot gain QS functionality by starting out with v 2 = 0. The end result is the extinction of Population 2, while Population 1 remains stable.
Overall, for h ≤ 2, we found that one population is driven to extinction, while the other stays at a stable level with QS intact at lower levels.
When h > 2, B(v i ) falls into the second type of model functions. That means there might be a stable positive strategy in addition to a stable zero solution. We found this to be the case for h = 3 for our parameters, as one can see in Figure 4. This means that Populations 1 and 2 remain at stable population and strategy levels, with Population 1 being QS active while Population 2 consists of cheaters. In this scenario, producers and non-producers can live side-by-side indefinitely.

Discussion
We have applied the G-function ansatz to the QS dynamics of P. aeruginosa and analysed the resulting model both analytically and numerically. We have found that the type of G-function considered changes the expected outcome drastically.
If we assume that there is no private benefit to QS at all, then the long-term outcome is preassigned and independent of parameters. Whether QS is inactive in all bacteria after some time, or the subpopulation of secretors dies out, the effect remains the same: QS is evolutionary unstable.
As soon as we assume that there is some kind of private benefit, the situation changes. Depending on parameter values, QS might still be unstable, but there can also be positive evolutionary stable secretion values. For the first type of model functions, there is only one stable secretion value. It follows that in the long term the population will unitise (Figure 4). The situation remains much the same if, instead of a direct benefit, we consider a reduction in death rate (results not shown). In contrast, private benefit functions of the second type admit bistable behaviour. This allows secretors and cheaters to coexist indefinitely in addition to single subpopulation states. We have also found that bacteria are driven towards QS if there is no pre-existing cooperation in their environment.
Since the left-hand-side of Equation (12) has a value of zero for v i = 0, whereas the right-hand-side has a value greater than zero for v i = 0, there could be any number of stable strategies (or none) if one chooses different functional responses.
Our model functions have simplified the actual process of QS a great deal, but one could also include direct dependencies on the signal molecules and enzymes produced. Such models with abiotic components were considered by Cohen et al. [30]. While their short-term dynamics might differ from the scenarios we have considered, their long-term outcome does not. They are also difficult to analyse without resorting to numerical methods, which is why we have concentrated on the simplified functions.
The G-function has allowed us to determine the long-term evolution of QS by considering the short-term population dynamics. By analysing different functional responses, we have explored some possible outcomes. One could also further refine the population dynamics for different real world scenarios. The G-function serves to always model the resulting evolutionary pressure, which tells us something about the way mutation and selection are going to drive the population. As both are happening on a fast time-scale in bacteria, they have a big impact on the interaction between bacterial populations as well as on their interaction with humans (see, e.g., [3]). The G-function can thus be a valuable tool in modelling bacterial interactions.
For completeness, in Appendix C, we show the spatial evolution of two populations with different start strategies, using a G-function without additions ( Figure A1) and with private benefit and a hill factor of h = 1 ( Figure A2).
Lastly, we would like to refer the interested reader to the following associated publications which complement this study: Mund et al. [34] and Mund et al. [35].

Conflicts of Interest:
The authors declare no conflict of interest. We first prove that R(v i ) has exactly one maximum as well as calculate its maximal value. To that end, we take the derivative of R

Appendix A. Variables and Parameters
and search for critical points For h > 3, there exists a critical point at v i = 0. However, since R(0) = 0 and R(v i ) > 0 for v i > 0, this critical point is clearly a minimum on the bounded space v i ∈ R + 0 . The other critical point satisfies and is thus unique in R + 0 . It remains to show that v * i is indeed a maximum. However, knowing that R(0) = 0 and lim v i →∞ R(v i ) = 0, we can conclude that it can only be a maximum.
We prove Property 5 by induction.
base case We can derive from the definition of R(v i ) that a h and thus f (0) = ha 2h (h − 2) > 0. inductive step Assuming that the statement holds for R (n) (v i ) and that n + 1 ≤ h − 2, we have and thus g(0) = a h (h − (2 + n)) f (0) > 0 since n < h − 2 and f (0) > 0.
It follows that R (n) (0) = 0 for n < h − 2 and R (h−2) (0) > 0. At the same time, a similar argument shows that Looking at the second derivative of G, we have from which we can gather thatv i = 0 is stable if B (0) < 2KB(v, b). We use Equation (14) and substitute for B(v, b) + B(v i ) As B (v i ) > 0, stability of positivev i immediately follows from Equation (A3).

Second Type of Model Terms
We start out by assuming B(v, b) > 0. We know that v i = 0 is stable if B (0) < 2KB(v, b) and that for h ≥ 3 it holds that B (0) = 0. Hence, v i = 0 is an ESS regardless of parameter values.
A positive stationary solutionv i is stable, if We can rearrange this inequality to The larger of both positive stationary solutions fulfilsv i > v * i = h h−2 h+2 · a, if it exists. It follows that Hence, the larger of both positive solutions is always stable, if it exists. It follows immediately from the stability of the zero solution that the smaller of both positive solutions must be unstable.
In the other case, we first note that both B (0) = 0 and B(v, b) = 0, which means we cannot gain insight into the stability of the zero solution directly through Equation (A4). However, we know that, for small enough ε > 0, R(ε) > 2KB(ε) and that there exists only one positive ξ for which 2KB (ξ) = R (ξ). Additionally, we know that ξ lies in the open interval between zero and the positive intersectionv i and as such is unequal tov i . It can thus follow that which is exactly the condition for stability from Equation (A4). For B(v, b) = 0, we thus have an unstable zero solution and a stable positive strategy.

Appendix C. Spatial Evolution of Two Populations with Different Start Strategies
In this section, we show the spatial evolution of two populations with different start strategies, using a G-function without additions ( Figure A1) and with private benefit and a hill factor of h = 1 ( Figure A2).  Figure A2. Evolution of two populations with different start strategies, using a G-function with private benefit and a hill factor of h = 1. Note the shorter timescale in this plot.