Epidemiological dynamics of bacteriocin competition and antibiotic resistance

Bacteriocins, toxic peptides involved in the competition between bacterial strains, are extremely diverse. Previous work on bacteriocin dynamics has highlighted the role of non-transitive ‘rock–paper–scissors’ competition in maintaining the coexistence of different bacteriocin profiles. The focus to date has primarily been on bacteriocin interactions at the within-host scale (i.e. within a single bacterial population). Yet in species such as Streptococcus pneumoniae, with relatively short periods of colonization and limited within-host diversity, ecological outcomes are also shaped by processes at the epidemiological (between-host) scale. Here, we first investigate bacteriocin dynamics and diversity in epidemiological models. We find that in these models, bacteriocin diversity is more readily maintained than in within-host models, and with more possible combinations of coexisting bacteriocin profiles. Indeed, maintenance of diversity in epidemiological models does not require rock–paper–scissors dynamics; it can also occur through a competition–colonization trade-off. Second, we investigate the link between bacteriocin diversity and diversity at antibiotic resistance loci. Previous work has proposed that bacterial duration of colonization modulates the fitness of antibiotic resistance. Due to their inhibitory effects, bacteriocins are a plausible candidate for playing a role in the duration of colonization episodes. We extend the epidemiological model of bacteriocin dynamics to incorporate an antibiotic resistance locus and demonstrate that bacteriocin diversity can indeed maintain the coexistence of antibiotic-sensitive and -resistant strains.


Introduction
Bacteriocins are toxic peptides that allow bacteria to eliminate competitors. Bacteriocins systems are pervasive in bacterial species and are thought to play a significant role in competition within (and possibly between [1]) species [2]. A mechanistic understanding of the role of bacteriocins in competition is therefore important for characterizing the ecological dynamics of bacteria.
At a high level of abstraction, bacteriocin systems can be thought of as consisting of three components: a system to produce and secrete toxins, a system to achieve immunity against these toxins, and a regulatory system to control toxin and immunity production-which may involve release of signalling molecules ( pheromones) allowing quorum-sensing and kin recognition. Even at this level of abstraction, bacteriocin systems are highly diverse. A single species may have a number of distinct systems, with considerable diversity within each of these: coexistence of strains with different combinations of toxin, immunity and regulatory genes is pervasive [3].
For example, the bacterial species Streptococcus pneumoniae has multiple different bacteriocins systems. The two best described systems-the blp (bacteriocin-like peptide) locus [4][5][6][7] and the cib (competenceinduced bacteriocin) locus [8]-are both ubiquitous. The specific combination of genes that make up the blp locus is highly variable between strains, while the cib locus is more conserved. Both loci are associated with a pheromone signalling system, which regulates the expression of toxins and immunity. These signalling systems are also diverse: there are multiple distinct 'pherotypes' (i.e. specific pheromonereceptor pairs allowing targeted signalling to cells of the same pherotype [9,10]). By contrast, other pneumococcal bacteriocin systems, such as the pld locus [11] and the circular toxin pneumocyclicin [12] are found on only a subset of strains. Thus, to fully capture bacteriocin competition, models should be able to explain both the diversity in bacteriocin profiles and the variation in the specifics of this diversity between different bacteriocin systems.
Previous theoretical and experimental work has highlighted the role of non-transitive competition in maintaining bacteriocin diversity. Indeed, two-strain models of toxin-producing ( producer) and toxin-susceptible (non-producer) strains do not typically predict coexistence (with some exceptions, see [13]): depending on the effectiveness of the toxins and the cost associated with their production, either the producer or the non-producer strategy out-competes the other [14,15]. Coexistence can be achieved by inclusion of a third (immune) strain which does not produce toxins but does have immunity against them (or, alternatively, produces less effective and less costly toxins [15]). If both toxin and immunity incur a fitness cost, the total cost to the producer strain (which has both toxins and immunity), is greater than the total cost to the immune strain (which only incurs the cost of immunity, but not toxin production). As a consequence, in head-to-head competition, the non-producer out-competes the immune strain; the immune strain out-competes the producer; and the producer out-competes the non-producer. This type of non-transitive competitive structure is referred to as 'rock-paper-scissors' dynamics [16].
To provide additional clarity in the context of bacteriocin dynamics, we make a distinction between two types of rockpaper-scissors model. In the most general formulation of rock-paper-scissors dynamics (figure 1a) [16,18], strain interactions are considered in terms of head-to-head competition for the occupation and invasion of patches or sites. Strains are differentiated through their ability to invade occupied patches: each strain wins against one of the other strains and loses against the other. Apart from the outcomes of this head-to-head competition, there are no further ecological differences between the strains. Models building on this structure have been studied in a number of ecological contexts [16], including bacteriocin dynamics [18]-where, at the epidemiological scale, patches would represent hosts to colonize, or, at the within-host scale, space to occupy in the colonized niche. Such models predict oscillatory dynamics of the three strains in well-mixed environments and stable coexistence in spatially structured environments [16,18].
This general formulation of rock-paper-scissors competition contrasts with bacteriocin-specific rock-paper-scissors models (figure 1b) [15,17] that explicitly represent bacteriocin-related processes within a host (or within a liquid culture or Petri dish). There are two key differences between these models. The first is in the interaction between producer and non-producer strains: in the bacteriocin-specific formulation, toxins lead to the death of non-producer cells, rather than their direct replacement by the producer cells. The death of non-producer cells frees up resources (e.g. nutrients or space), which can then be used by any of the strains-or any of the locally present strains in a spatially structured model. As a result, the immune strain can gain the same benefit from the production of toxins as the producer strain, without paying the cost (a 'cheater' strategy). The second difference is how the cost of immunity and bacteriocin production are represented: these are not modelled in terms of strain interactions, but rather as a reduced growth rate or increased death rate for the immune and producer strains. These bacteriocin-specific models predict stable coexistence of the three strains in spatially structured environments, but dominance of a single strain in well-mixed systems [15]. This prediction has been verified in an experimental model of producer, non-producer and immune strain interactions in Escherichia coli [17]. More recent modelling suggests that when the benefits of immunity are shared, rather than specific to the immune strain (e.g. immunity involves secretion of toxindegrading compounds), coexistence of strains can also arise in well-mixed environments [19]-although others have suggested that coexistence in such models is not generally robust to invasion by a non-producer cheater strain [20].
These bacteriocin-specific models capture within-host interactions between different bacteriocin profiles and thus provide insights into dynamics at this scale. Yet, within-host interactions will also impact dynamics at the epidemiological scale by allowing strains to invade already colonized hosts and/or by preventing such invasion. In addition, any fitness  Figure 1. Strain interactions arising in different models of bacteriocin dynamics. P represents the producer strain (both toxins and immunity), M the immune strain (immunity only) and N the non-producer (no toxins, no immunity). (a) A general formulation of rock-paper-scissors competition. Strain interactions are modelled in terms of the outcome of head-to-head competition between two strains. The costs of toxins and immunity, as well as the killing action of the toxins, are modelled implicitly in the outcome of this head-to-head competition. At the epidemiological scale, assuming fitness differences between strains are only apparent when the strains are competing for the same host gives rise to this model structure. (b) Bacteriocin-specific within-host model [15,17]. The producer strains kill non-producer strains, leading to resources being freed up. The costs of immunity and bacteriocin production are modelled as decreased growth rate or increased death rate.
(c) Bacteriocin-specific epidemiological model. The killing of the non-producer strain by the producer strain leads to the host becoming colonized with the producer strain. In this epidemiological model, the costs associated with bacteriocin production and immunity can be represented in terms invasion probabilities and/or epidemiological parameters (i.e. transmission and/or clearance rates).
royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 289: 20221197 costs associated with bacteriocin production and immunity may also impact epidemiological parameters (i.e. transmission and clearance rates). In species such as Streptococcus pneumoniae, these effects at the epidemiological scale are likely to play an important role in shaping the ecology and evolution of bacteriocin systems, for two reasons. Firstly, within-host models assume all three species are co-colonizing the same host. However, multiple colonization is not necessarily the norm: for example, reported rates of co-colonization with multiple pneumococcal strains are around 25% to 50% depending on setting [21][22][23]. Secondly, when the duration of colonization is relatively short (e.g. of the order of months in S. pneumoniae), factors affecting transmission will play a large role in shaping allele frequencies in the overall bacterial population. The impact of bacteriocins at the epidemiological scale is not well established and may depend on the specific bacteriocin system. For example, the cib locus provides established colonizers an advantage over invaders in a murine model of pneumococcal colonization [24], suggesting these bacteriocins play a defensive role. On the other hand, an observational study of pneumococcal strains co-colonizing humans suggests that the blp locus has a limited impact on which strains can co-colonize hosts [25]. Furthermore, the impact of the cost of bacteriocin production and immunity on epidemiological parameters and the effects of kin recognition are unclear. Exploring bacteriocin dynamics at the epidemiological scale therefore requires a flexible modelling approach that can capture a range of effects.
While the general rock-paper-scissors model (figure 1a) can be interpreted as an epidemiological model (with patches corresponding to hosts), previous work has focused on patch invasion and not explored the possibility that toxin-production and immunity associated costs might impact epidemiological parameters. Bacteriocin-specific within-host models (figure 1b) do not accurately capture interactions at the epidemiological scale: at the within-host scale, toxins lead to the death of the non-producer and freeing up of resources, while at the epidemiological scale, toxins lead to invasion of the host occupied by the non-producer and thus only benefit the producer itself. To some extent, spatial structuring at the within-host scale can limit the benefit of freed-up resources to the producer strain, but this remains different from interactions at the epidemiological scale: the benefit is only private when all neighbouring cells are also producers-which cannot be true across the whole population if multiple strains coexist. Thus, bacteriocin-specific models of non-transitive competition at the epidemiological scale (figure 1c) are necessary for understanding bacteriocin dynamics and diversity in S. pneumoniae and species with similar ecology.
A further motivation for developing epidemiological models of bacteriocin dynamics is their potential role in determining the duration of colonization episodes. We have previously suggested that duration of colonization affects selection pressure for antibiotic resistance in species which are carried asymptomatically most of the time, such as S. pneumoniae [26]. Strains that colonize hosts for longer gain a greater fitness benefit from resistance, leading to an association between long duration of colonization and antibiotic resistance (including multi-drug resistance [27]). As a result, balancing selection maintaining diversity at a locus that affects duration of colonization could also maintain diversity at a resistance locus [26], providing a potential explanation to the long-standing puzzle of why antibiotic resistance has not yet reached fixation [28]. However, this explanation is not complete, because the genetic determinants of duration of colonization have not been fully identified [29]. A more complete understanding of resistance dynamics therefore requires identifying loci that contribute to variation in duration of colonization. The role of toxin production in killing competing bacteria and the role of immunity in preventing this killing makes bacteriocin loci a promising candidate.
This paper is organized in two parts: the first part (section 2) develops an epidemiological model of bacteriocin dynamics in S. pneumoniae (or a species with similar ecology). We explore the circumstances in which this model allows coexistence of strains with different bacteriocin profiles. In the second part (section 3), we investigate differences in duration of colonization between bacteriocin profiles and their impact on resistance dynamics.

Bacteriocin dynamics (a) Epidemiological modelling of bacteriocin dynamics
We begin by considering a bacteriocin system with two components: a toxin gene and an immunity gene. There are therefore three possible strains: a producer strain (P), with both toxin and immunity; an immune strain (M), with immunity but no toxin; and a non-producer strain (N), with neither toxin nor immunity (i.e. the entire bacteriocin locus is absent). Strains with the toxin but no immunity are assumed to be lethal to themselves and therefore not included. We formulate a generic epidemiological model of competition between three strains rather than focusing on what is known about a specific bacteriocin system, and consider how ecological differences between bacteriocin profiles can be reflected in this generic model.
In this model formulation, strain i colonizes uncolonized individuals X at rate β i and is cleared at rate μ i . In addition, strains can invade already colonized individuals and displace the resident strain at rate β i k ij for strain i replacing strain j, where k represents the probability of successful invasion relative to colonization of an uncolonized host. We assume the dynamics of this replacement are fast; the invading strain is therefore modelled as replacing the resident strain instantaneously. Hosts are therefore only ever colonized with a single strain at a time and there is no co-infection. Note that our qualitative results are robust to relaxing this assumption (electronic supplementary material, section 2.1). The dynamics of this general model are described by the following equations: Here, the X and I variables represent proportions of the total host populations: X + I P + I M + I N = 1.
Differences between bacteriocin profiles (i.e. producers P, non-producers N and immune strains M) can be reflected in two ways in this model (figure 1c). Firstly, through asymmetric invasion probabilities k (0 ≤ k ≤ 1) arising from differences in within-host fitness: toxins allow P to invade N more often than the other way around (k PN ≥ k NP ); similarly, the cost of toxin production means M out-invades P (k MP ≥ k PM ); and the cost of immunity means N out-invades M (k NM ≥ k MN ). Secondly, the costs of toxins and immunity could also affect epidemiological parameters-i.e. transmission and clearance rates: These inequalities reflect typical assumptions about bacteriocin dynamics and, unless otherwise indicated, we have assumed they hold in our analysis.
In this structure, the impact of toxins (and other withinhost fitness differences) can be modelled as either offensive (i.e. enabling invasion of already colonized hosts) or defensive (i.e. preventing invasion by another strain) or a combination of both. Figure 1 depicts offensive interactions: invasion is only possible when the invading strain has a fitness advantage. This is also how we parameterize the model in main text results (i.e. k NP = k PM = k NM = 0). It is worth noting that the distinction between offensive and defensive only impacts epidemiological dynamics if the transmission rate differs between strains: when transmission rates are equal, the interaction between two strains-say strain A and strain B- The dynamics therefore depend only on the relative rates of invasion. Under these circumstances, modelling bacteriocins as offensive is mathematically equivalent to modelling them as defensive. This equivalence does not hold when transmission rates are not equal (i.e. when epidemiological costs affect transmission rate). However, we find that in practice, the distinction between offensive and defensive bacteriocins does not have an impact on which strains are observed to coexist (electronic supplementary material, section 2.2).

(b) Epidemiological differences between strains allow a wider range of outcomes
Our aim is to understand the circumstances under which multiple strains coexist at equilibrium in this system. We address this through linear stability analysis (unless otherwise indicated). In general, we parameterize clearance rate as μ = 1 and transmission rate as β = 3. In time units of month −1 , these are plausible values for S. pneumoniae [28]. All analyses and simulations were performed using Wolfram Mathematica [30]; the code is available in the electronic supplementary material. First, it is useful to note that if strains differ only in invasion probabilities (the k parameters) but not epidemiological parameters (i.e. μ P = μ M = μ N and β N = β M = β P ), the model is structurally identical to the classic rock-paper-scissors model (e.g. [16]; see figure 1). The behaviour of this system is well characterized: as long as the non-transitive competitive structure is maintained, all three strains coexist, with oscillatory dynamics around a stable equilibrium point (see also figure 2). When the non-transitive competition structure is not present, a single strain dominates. Equilibria involving just two of the strains are never stable.
Introducing differences in epidemiological parameters (i.e. transmission rate β or clearance rate μ) between bacteriocin profiles has two effects. Firstly, oscillation is not observed when the strains differ in epidemiological parameters (figure 2; electronic supplementary material, figure S7). Thus, like spatial structure in previous models [15][16][17], differences in the epidemiological parameters of strains act to stabilize coexistence. Secondly, the range of possible outcomes increases. Stable coexistence no longer requires all three strains; for some parameter ranges, the producer and non-producer strain can coexist without the immune strain ( figure 2b).
It is worth noting that the properties of the two and three-strain equilibria are different. The general rockpaper-scissors model ( figure 1a) is known to give rise to a 'survival of the weakest' effect: increasing the competitive advantage of a strain (i.e. ability to invade) decreases its equilibrium frequency [16]. This occurs because of the non-transitive competitive structure. For example, a smaller k PN and thus less frequent displacement of N by P increases the frequency of N and therefore the rate at which M is displaced. This, in turn, decreases the frequency of M and displacement of P, leading to a higher frequency of P despite its lower competitive advantage. We find that this effect also applies to the relationship between transmission rate and equilibrium frequency (electronic supplementary material, figure S6 and section 2.3): increasing the transmission rate of a strain decreases its equilibrium frequency when all three strains are present at equilibrium. This survival of the weakest effect is not observed for two-strain equilibria.
Finally, the number of stable strain combinations is even greater if the bacteriocin-producer can be associated with a lower cost than the immune strain: this allows the producer to exclude the other two strains (figure 2b). Such circumstances would arise when bacteriocin-production has benefits beyond inter-strain competition-such as bacteriocidal action against other species, leading to the producer having increased success in colonizing hosts which are not carrying the focal species. This would result in an increased transmission rate for the producer strain, thus offsetting some of the cost of bacteriocin production.

(c) The effects of kin recognition ( pherotype)
To explore the effect of kin recognition on bacteriocin diversity, we expand the original model to include two different signalling molecules and corresponding receptors (i.e. two distinct pherotypes), yielding a six-strain model. The dynamics between strains of the same pherotype remain as described in equation (2.1). Across pherotypes, immunity is not effective: the signal to turn on expression of immunity proteins is not recognized. Producer strains are therefore able to out-compete immune strains with a different pherotype; in other words, the interaction between producer and immune strain is the same as the interaction between producer and non-producer strain. Under the assumption that bacteriocins are required for invasion, this results in the model depicted in figure 3a (see electronic supplementary material, section 1.1 for equations and additional discussion).
The inclusion of pherotype in the model decreases the parameter space in which coexistence is observed: the producer strain excludes the other strains in a large part of the explored parameter space (figure 3b). The inclusion of royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 289: 20221197 pherotype also increases the number of combinations in which strains can coexist, allowing stable equilibria with the producer and immune strains without the non-producer strain. These effects arise because the addition of pherotype expands the range of strains that are susceptible to the toxins, thus allowing the producer strain to exist without the non-producer strain.

Potential role for bacteriocins in resistance dynamics (a) Bacteriocins profiles differ in duration of colonization
We now turn to the potential role of bacteriocins in the dynamics of antibiotic resistance. We have previously  The inclusion of pherotype further increases possible outcomes but decreases the parameter space in which multiple strains coexist. As in figure 2, the x-axis represents how well the producer strain invades (k P ¼ k P where the A and B superscripts indicate pherotype); the y-axis represents the transmission cost for the producer strain; and the colours indicate which strains are present at equilibrium. Other parameters are β = 3, μ = 1 (in time units of month −1 , these are plausible values for S. pneumoniae [28]), c M = 0.1 and all other k parameters 0, similar to figure 2b. Due to the difficulty of computing equilibrium solutions and eigenvalues for this model, these results are derived through simulating rather than stability analysis (simulation until t ¼ 10 5 months).
royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 289: 20221197 suggested that the duration of colonization affects selection pressure for antibiotic resistance: modelling predicts that fitness effect of resistance depends on a strain's duration of colonization, with longer colonization being associated with a greater benefit from resistance [26]. Indeed, a strain's duration of colonization correlates with the frequency of resistance within the strain in multiple pneumococcal datasets [26,31]. Variation in the duration of colonization among coexisting strains could therefore maintain the coexistence of antibiotic sensitive and resistant strains, with long-colonizers being antibiotic resistant and short-colonizers antibiotic sensitive. Bacteriocins, due to their inhibitory effects on others strains, are a plausible source of such stable variation in duration of colonization.
To assess the potential role of bacteriocins in resistance dynamics, we examine the extent of variation in duration of colonization between different bacteriocin profiles. In the model represented by equation (2.1), the mean duration of colonization of strain i (D i )-i.e. the inverse of its overall clearance rate-is given by

ð3:1Þ
Bacteriocin profiles can thus affect duration of colonization through two mechanisms: (i) directly through any effects of the cost of toxin production and immunity on clearance rate and (ii) through strain displacement, the impact of which depends on invasion rates (k), the effect of the cost of immunity and toxin production on transmission rate, and the prevalence of the invading strain. As a result, the extent of variation in duration of colonization and which strain is associated with the longest colonization depends on parameter values (figure 4; electronic supplementary material, section 2.5). When costs affect transmission, the longest duration of colonization is associated with either the producer or immune strain, and the shortest with either the non-producer or immune strain. This is reversed when costs affect clearance, with the longest duration of colonization associated with the non-producer or immune strain, and the shortest with the producer or immune strain. In broad terms, the greatest variation in duration of colonization is observed when the cost associated with bacteriocin production is high. A more detailed discussion of the relationship between parameters and duration of colonization can be found in electronic supplementary material, section 2.5.

(b) The effect of bacteriocins on antibiotic resistance dynamics
The observed variation in the duration of colonization of different bacteriocin profiles suggests that diversity in  We model a species that is carried asymptomatically most of the time (e.g. S. pneumoniae)-the antibiotic exposure of hosts is therefore independent of whether they are colonized and equal to the average antibiotic consumption rate in the population (τ). Interactions between the three bacteriocin profiles are the same as previously described. In addition, each profile can be either antibiotic sensitive (S) or antibiotic resistant (R), giving rise to six possible strains. Antibiotic sensitive strains are subject to an additional clearance rate τ (we assume immediate clearance in response to antibiotics). Resistance carries a fitness cost, which can be associated with clearance and/or transmission rate. The full model structure is given in electronic supplementary material, section 1.2.
Consistent with previous results, the between-strain variation in duration of colonization allows coexistence of antibiotic sensitivity and resistance. Antibiotic resistance is associated with the bacteriocin profile(s) with longer duration of colonization ( figure 5). The range of antibiotic consumption rates over which coexistence is observed and the frequency of antibiotic resistance in this region of coexistence is highly dependent on parameters. The finding that diversity at the bacteriocin locus can maintain intermediate antibiotic resistance frequencies is generally robust, with the exception of some specific cases when the cost of antibiotic resistance affects clearance rate. These are discussed in detail in electronic supplementary material, section 2.6.

Discussion
This paper was motivated by two main questions: firstly, the role of epidemiological processes in bacteriocin dynamics and diversity in S. pneumoniae and similar species; and secondly, whether these epidemiological processes could generate variation in duration of colonization between different bacteriocin profiles and, as a result, contribute to the observed coexistence of antibiotic sensitive and resistant strains.
We have shown that bacteriocin diversity is readily maintained in epidemiological models. This finding is robust to assumptions about how the inhibitory activity of bacteriocins and costs associated with bacteriocin production and immunity translate to the epidemiological scale (e.g. within-host fitness differences leading to one strain displacing another versus differences in epidemiological parameters versus a mixture of the two): coexistence of different bacteriocin profiles is found under a range of assumptions. We find that the specifics of this coexistence (i.e. which bacteriocin profiles coexist and at what frequencies) are highly dependent on parameter values. This sensitivity to the properties of the bacteriocins-such as the cost of toxin productionmay explain the observed variability in characteristics (e.g. prevalence and composition of bacteriocin loci) of bacteriocin systems.  [28]. Due to the difficulty of computing equilibrium solutions and eigenvalues for this model, these results are derived through simulating rather than stability analysis (simulation until t ¼ 10 5 months).
royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 289: 20221197 In addition to the robust diversity in bacteriocin profiles in epidemiological models, we can also draw other general insights from these results.
Firstly, the inclusion of transmission dynamics introduces additional mechanisms through which strain diversity can be maintained, allowing different combinations of strains to coexist. Previous work at the within-host scale had suggested that all three strains ( producer, non-producer and immune) must be present to allow diversity to persist: the maintenance of diversity depended on a non-transitive competitive structure (rock-paper-scissors dynamics) [15,17]. Coexistence of the producer and non-producer has previously been observed in a two-strain meta-population model, through a mechanism related to the dynamics of co-colonization [13]. In our three-strain epidemiological model, coexistence of the producer and non-producer strain without the immune strain is also possible. The mechanism at play here is neither rock-paper-scissors dynamics nor co-colonization but rather a competition-colonization trade-off [32]: the producer strain has an advantage at the within-host scale (i.e. it can displace the non-producer) but a disadvantage at the between-host scale (i.e. it pays a fitness cost which reduces transmission rate and/or increases clearance rate). This type of trade-off gives rise to negative frequency-dependent selection-the benefit of the within-host advantage depends on the prevalence of the non-producer strain-and thus maintains coexistence of the competitors [33].
Secondly, in S. pneumoniae, some bacteriocin systems, such as the pld locus [11] and the circular toxin pneumocyclicin [12], are found on only a subset of strains. Others, such as the blp [4][5][6][7] and cib [8] loci, are ubiquitous. Such ubiquity is not consistent with our epidemiological model under standard assumptions about bacteriocin ecology: all stable equilibria require the presence of the non-producer strain, because toxin production is only beneficial when susceptible strains are present. Our results highlight two mechanisms which allow toxin-producers to exist in the absence of nonproducers (in other words, mechanisms that allow the bacteriocin locus to fully invade the bacterial population). Firstly, the producer can exclude the other two strains if it beats the immune strain in head-to-head competition-such circumstances could arise if toxins are also effective against other bacterial species present in the host, thus providing a benefit beyond their effect on the non-producer. Secondly, in the pherotype model, immunity is not protective against toxins released by a different pherotype, and strains of one pherotype are therefore susceptible to toxins from another pherotype. Toxin production is therefore beneficial even in absence of the non-producer strain. Indeed, the ubiquitous cib and blp loci are both associated with a signalling locus (i.e. pherotype), suggesting this may be the explanation for their ubiquity.
In section 3, we have shown that bacteriocins are a plausible candidate for involvement in resistance dynamics. We predict differences in duration of colonization of different bacteriocin profiles, with the magnitude of these differences highly sensitive to parameter values. The duration of colonization is predicted to modulate the fitness of antibiotic resistance [26]; we therefore expect an association between antibiotic resistance and bacteriocin profile and indeed observe this in a model incorporating both bacteriocin and resistance dynamics.
Testing this predicted association empirically is possible in theory, but would prove challenging given our current understanding of bacteriocin systems. Firstly, the direction of association between bacteriocin profile and antibiotic resistance depends on parameters. It is therefore unclear whether we expect antibiotic resistance to be associated with the producer, non-producer or immune phenotype, and indeed, this may differ between different bacteriocin systems. Secondly, the mapping between these modelled phenotypes and observed genotypes is non-trivial: although toxin and immunity genes can be identified, the effect on phenotype is less clear when multiple toxin and immunity genes are present on a genome (either because of the presence of a multi-toxin bacteriocin system or the presence of multiple different systems on one genome). Furthermore, in some bacteriocin systems, additional factors may affect the association between toxin and immunity genes; for example, a large proportion of blp systems may not be able to secrete pheromones and toxins because of an impaired transporter [34] and there is evidence of regulatory interplay between the blp and cib systems [35,36]. Thus, although we have shown bacteriocins are a plausible candidate for involvement in resistance dynamics, more specific predictions about this involvement will require a more detailed understanding of specific bacteriocin systems.
It is worth highlighting two potential limitations of our modelling approach. Firstly, we do not model co-colonization: within-host strain displacement is assumed to be fast. This assumption may not hold: co-colonization with multiple strains is known to occur (e.g. 25-50% of colonized hosts depending on setting [21][22][23]), although whether such co-colonization is possible between strains with different bacteriocin profiles is not clear. Furthermore, previous theoretical work suggests that assumptions about within-host dynamics and competition can have considerable impact on predictions about strain diversity at the between-host scale [37][38][39]. We tested the impact of allowing slower within-host dynamics and thus co-colonization in all three models (bacteriocin dynamics, pherotype and antibiotic resistance) and find our results are qualitatively robust (electronic supplementary material, section 2.1).
Secondly, we consider toxin-production and immunity (and antibiotic resistance) as binary traits: strains either possesses genes encoding for toxins and immunity or they do not. This modelling approach reflects observed variation in the absence/presence of toxin and immunity genes in pneumococcal genomes [4][5][6][7][8]11,12]. An alternative approach would be to treat these traits as continuous variables; this would correspond to assuming genes are present on all strains and that effects are tunable (e.g. through modulation of gene expression), leading to different levels of toxicity, immunity and fitness cost. This would allow modelling the evolution of the level of toxicity and immunity. Such approach would require assumptions about how tunable bacteriocin effects are; trade-offs between fitness cost and levels of toxicity and immunity; and the relative time scale of evolutionary and epidemiological processes. Such assumptions will be subject to considerable uncertainty. For example, the extent to which bacteriocins are tunable is unclear-indeed, the relationship between inhibition and toxin concentration appears threshold-like rather than gradual [7], suggesting inhibitory effects may not be readily modulated. Nevertheless, the impact of assuming bacteriocin traits are evolvable royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 289: 20221197 on a similar time scale to epidemiological processes is an open and interesting question.
In summary, we have shown that diversity in bacteriocin profiles arises robustly in epidemiological models of bacteriocin dynamics and can be maintained through either rock-paper-scissors dynamics or a colonization-competition trade-off. The specifics of the predicted diversity are sensitive to assumptions about how bacteriocins affect epidemiological processes and to bacteriocin-related parameters (e.g. costs of bacteriocin production and immunity), providing a potential explanation for differences between bacteriocin systems. We have also demonstrated that diversity at bacteriocin loci is a plausible candidate for also maintaining diversity at resistance loci. These insights arise from modelling that approaches bacteriocins dynamics from a high level of abstraction, rather than representing a specific bacteriocin system. Generating more specific insights into particular bacteriocins systems will require models informed by the biology of the specific system. Therefore, a more complete understanding of the role of bacteriocins in bacterial ecology will need a more specific characterization of their effects on transmission, invasion, within-host competition and duration of colonization.
Data accessibility. The code and supplementary information are available as electronic supplementary material [40].