Oscillations in continuous culture populations of Streptococcus pneumoniae: population dynamics and the evolution of clonal suicide

Agents that kill or induce suicide in the organisms that produce them or other individuals of the same genotype are intriguing puzzles for ecologists and evolutionary biologists. When those organisms are pathogenic bacteria, these suicidal toxins have the added appeal as candidates for the development of narrow spectrum antibiotics to kill the pathogens that produce them. We show that when clinical as well as laboratory strains of Streptococcus pneumoniae are maintained in continuous culture (chemostats), their densities oscillate by as much as five orders of magnitude with an apparently constant period. This dynamic, which is unanticipated for single clones of bacteria in chemostats, can be attributed to population-wide die-offs and recoveries. Using a combination of mathematical models and experiments with S. pneumoniae, we present evidence that these die-offs can be attributed to the autocatalytic production of a toxin that lyses or induces autolysis in members of the clone that produces it. This toxin, which our evidence indicates is a protein, appears to be novel; S. pneumoniae genetic constructs knocked out for lytA and other genes coding for known candidates for this agent oscillate in chemostat culture. Since this toxin lyses different strains of S. pneumoniae as well as other closely related species of Streptococcus, we propose that its ecological role is as an allelopathic agent. Using a mathematical model, we explore the conditions under which toxins that kill members of the same clone that produces them can prevent established populations from invasion by different strains of the same or other species. We postulate that the production of the toxin observed here as well as other bacteria-produced toxins that kill members of the same genotype, ‘clonal suicide’, evolved and are maintained to prevent colonization of established populations by different strains of the same and closely related species.


Mathematical Model. Analysis
The simple version of the model is described by the following set of equationṡ where, Ψ(R) = νR/(K + R) These equations can be rewritten in dimensionless form by rescaling the resource (R), bacterium (B), and toxin (T) variables as: This gives:Ṙ with ψ(R ) defined as ψ(R ) = R /(R + κ) Here we have also rescaled the time variable as t = ωt. The dimensionless combinations of parameters which determine the system's dynamical behavior are: f = (ω + d)/ω, κ = K /C, a = ef ν/yC, b = ν/ω.
The possible stationary points (equilibrium points*) of this system, R * , B * , and T * , are found by puttingṘ =Ḃ =Ṫ = 0, and solving the ensuing algebraic equations. There are three such stationary points.
From this point we dropped the superscript prime (') on each variable. B → B, R → R, T → T . S-2 The first -"resource only"-such point has: B * = T * = 0, and R * = 1 A linearized stability analysis shows this point is stable if, and only if, Failing this, an arbitrary small introduction of bacterium will grow (Ḃ > 0).
The second -"plant-hervibore"-point has This equilibrium is stable if, and only if, Within this domain, any introduction of T decays to zero, and the system returns to the R * , B * equilibrium with damped oscillations if κ > (b − 1) 2 /(2b − 1), or via exponential damping otherwise.
Third, and the most interesting, is a "3 species" (a "plant-hervibore-predator") equilibrium solution with Here, D is defined for notational convenience as D = (a + κ − 1)/2. This equilibrium solution is stable with all the three components persisting together, if and only if For the experimental system, representative values of the parameters are ω = 0.1 ml hour −1 , x = 5 × 10 −6 ml cells −1 a.u.t. −1 , e = 10 −7 µg, y = 4 × 10 −10 ml cells −1 a.u.t. −1 , K = 0.25 µg ml −1 , ν = 1 cells hour −1 , d = 0.1 hour −1 (the a.u.t are defined as "arbitrary units of toxin"). The parameter characterizing the resource/chemostat concentration, C, take various values, but always C >> K . Given that C is the control variable and κ represents it, the definition of a can be changed to a = ακ with α defined below. For the dimensionless parameters of Eq.(8) we thus have f = 2, b = 10, a = ακ (with α = 2000), and κ = 0.25/C, with C having values from a few tens to a few thousands. That is, κ << 1, whereas a can be less or greater than 1, and the defintion of α is as follows: and it is important to remember that according to our definitions K = κC.
The dynamics of this 3-species system (represented in Figure S1), as κ varies (corresponding to the chemostat resource concentration C varying), are interesting, with a rather sharp transition in dynamical behavior at the point where a = 1 − κ ≈ 1. For the illustrative parameters above, this corresponds to C = 500. For a smaller than this (a < 1), and remembering κ 1, for all C-values of interest, we have D = (a + κ − 1)/2 < 0, resulting in R * ≈ 1 − a, B * = 1, T ≈ (b − 1). In this region the system approaches its stable equilibrium point via very weakly damped oscillations. The oscillations have a period, P , approximatelly given by whilst the damping time, τ , is very long of the order Here, these expressions are given in absolute time (not rescaled), hence the factor ω.

S-4
The dynamical behavior changes abruptly when κ --although still much less than 1 -approaches the value such that a = ακ (with α 1) is essentially one: a = 1+ (a few multiples of)κ. In this tiny window, D ∼ K and R * ≈ κ 1/2 , while still B * = 1 and T * = b − 1.
The system still approaches equilibrium via damped oscillations, but now with the period (P ) and the damping time (τ ) are of comparable magnitude: P is much as above, but now Finally, for a > 1 (by a few multiples of κ), we move into a regime where R * ≈ κ/(a−1) << Perturbations to this equilibrium point decay back to it; the characteristic damping time As a increases, β and thence γ decrease, until we have γ = 0 when β = f − (f 2 − f ) 1/2 . Beyond this point, the system returns to its equilibrium point via damped oscillations, whose period (T ∼ 1/ | γ | ω) and damping time (τ ∼ 1/β) have comparable magnitudes. Eventually, once a becomes sufficiently large that a > b[1 − κ/(b − 1)], the 3-species equilibrium is no longer stable, and the toxin T can no longer be maintained in the system. Figure ?? illustrates this range of behavior, and in particular the very abrupt transition around the point a = 1 (which for the illustrative parameter values outlined above, corresponds to C = 500): for C-values above this, there are very weakly damped oscillations; below it, there is exponential damping if 341 < C < 500, and damped oscillations for , corresponding here to C < 50, the toxin can no longer be maintained, and we have a stable "resource -bacterium" system R * -B * . For implausibly small C-values (C < 1/36), the only stable state is the "resource-only" system, with neither bacterium nor toxin. Again, and at the risk of being repetitive, the composite parameteres a, b, and f allows us to examine the transition between the different dynamic behaviours that the system might experience as the result of changes in the resource concentration C.
It should be emphasized that the stability analysis presented above is, strictly speaking, valid only for small perturbations (i.e. a linearized analysis). But -although we have not found a Lyapunov function for the above system -previous experience suggests the global dynamical properties march with the local ones. Extensive numerical simulations support this suggestion. ). This raises the possibility that hydrogen peroxide could be the toxin responsible for the observed oscillations of S. pneumoniae R6 chemostat culture. The results of experiments with catalase reported in our article are inconsistent with this interpretation. Also inconsistent are the experiments we performed with a strain deficient in the production of pyruvate oxidase, an exzyme required for the production of hydrogen peroxide. In chemostat culture this strain oscillates, Figure S2.

Strains deficient in the induction of competence (comC and ComD deficient) and bacteriocin like protein (blpR deficient) systems oscillate
The competence peptide of S. pneumoniae, which appears to be produced in an autocatalytic manner similar to toxin postulated here, is a reasonable candidate for the agent (toxin or regulatory factor involved in the production of the toxin) driving the oscillations.

S-7
To test this hypothesis, we inoculated chemostats with mutants defective in the production of the competence peptide (∆-ComC) or in the receptor of the competence peptide (∆-ComD). These strains oscillate ( Figure S3). A similar behaviour is found in strains in which the regulatory gene of the blp system, also known to operate autocatalitically, has been knocked out. That is, it shows oscillatory dynamics similar to the parental strain ( Figure S4) 4 Strains deficient in the production of murein hydrolase, CbpD, oscillate As estated in the previous section the competence peptide of S. pneumoniae, which appears to be produced in an autocatalytic manner, does not oscillate. Additional evidence is found in a strain that lacks the production o the general response regulator ComE of the competence system (see Figure S5). Moreover, the lysing protein murein hydrolase, produced by CbpD, and seemingly regulated by the competence system, does not seem to prevent the strains from oscillating (see Figure S5), and we consider it can be ruled out as a potential oscillation driving agent in this system.

Assisted suicide protects established populations for invasion of higher fitness clones or species. Model and analysis.
In our report, we postulate that despite the killing of members of the same clone, assisted suicide could have evolved and be maintained as a mechanism to prevent established populations from invasion by bacteria of different clones or species with greater intrinsic fitness. The necessary condition for this is the invading clone is sufficiently sensitive to the toxin at its concentration in an established population of Streptococcus pneumoniae to have a net disadvantage relative to members of the established clone.
To begin exploring these invasion prevention conditions quantitative way, we expanded our model (equations 1-3) to include a second population of bacteria, B 2 and add a subscript 1 to the variable designating the toxin-producing species, B 1 . We assume both species are limited by the same resource which they consume in a Monod-like fashion with for simplicity the same values of k and e. These two populations can have different maximum growth rates, V 1 and V 2 , and different sensitivities to the toxin. With these definitions and assumptions, the rates of change in the densities of the different populations are given by,Ṙ where, Ψ(R) = νR/(K + R) as in the previous model. V 1 and V 2 are the maximum growth rates of these bacteria, and x 1 and x 2 are the killing rates of the producing and invading bacteria respectively. For convenience we assume that the efficiency of conversion of resource into bacteria biomass, e, and the Monod constant, k, are identical for B 1 and B 2 . As in the previous model y is the rate constant of production of the toxin and d is the inverse of the half-life of the toxin.
For analytical tractability in the above model we are assuming that B 2 does not produce the toxin. This assumption can be justified by our concern with the conditions for invasion of established populations of the toxin-producing strain B 1 by B 2 when it is rare. Even if the per-capita rate of production of the B 2 toxin was great and the established population S-10 highly sensitive to it, as long as the density of the invading population is low, the concentration of the toxin produced by it would be too low to cause a perceptible effect on the much more numerous established population.
Suppose B 2 tries to invade the R − B 1 − T system. That is, for entry/growth of a tiny amount of B 2 , we have:Ḃ where, R * , T * are given by the (equilibrium) solutions of the R − B 1 − T system being invaded.
Clearly, invasion is possible if, and only if, But ψ * is given by ν 1 ψ * = ω + x 1 T , so Eq.(23) becomes i.e with T * from our earlier analysis of R − B 1 − T system.
Clearly, if ν 2 ≈ ν 1 and x 2 > x 1 , this Eq.(25) can not be satisfied, and cannot invade. But, if ν 2 is sufficiently larger than ν 1 , then the right hand side (RHS) of Eq.(25) gets bigger AND To be more explicit, we need to put in the value of T * (which, remember, is related to our earlier, RESCALED, T * by T * = ω y T * ) for the R − B 1 − T system with "a > 1" (smaller C values), we have the system stable at T * ≈ b a − 1, whence T ≈ (ν 1 /y − ω)/y.
So, the condition for B 2 to be UNABLE to invade is (for a > 1) S-11 For R − B 1 − T with "a < 1" let us first look at the corresponding criterion for B 2 unable to invade the stable equilibrium point R * − B * 1 − T * , and second consider the case when, in practice, we have cycles of large amplitude, weakly damped. First, for "a < 1", This gives the UNABLE to invade criterion as: Notice that Eq.(27) has the larger denominator, (ν 1 − ω), than Eq.(26) when a > 1, where the corresponding factor is ν 1 a − ω , so it is easier to keep B 2 out when a < 1.
It is also important to notice also that the long time in which the system is cycling will create conditions in which T will be often (roughly half the time) less than T * . Even, so, it is likely that Eq.(26) remains a good approximation to the "exclusion criterion" even when we are in the very-weakly-damped cycle phase of the system. That is, θ > ψ to approximate neglect 10 −5 ! For Eq.(26): a < 1 (big C ↔ most of the simulations) (although complicated by oscillations) cannot invade iff θ > ψ + 8 9 × 10 −5 (ψ − 1) That is, we can assume the term 10 −5 as neglectable and the condition again reduces to θ > ψ

S-12
Even when more analysis can be done with this system our interest is to show that it is possible to find a condition under which a second species can not invade an established population that produces a self-killing toxin. The existance of such condition opens the gates to additional interesting questions on the way such system works and the toxin is maintained.

Toxin sensitivity of other species of Streptococcus
As noted in our report, the naturally occurring clinical isolates of S. pneumoniae we tested oscillate when put into chemostat culture. Based on the spot assay on lawns, other species of Streptococcus are also susceptible to this toxin. To ascertain the relationship between this sensitivity and the genetic relationship of these assayed species we used a maximum likelihood reconstruction protocol for the phylogenetic reconstruction and the parsimony reconstruction of the sensitivity to the toxin overimposed to it.
To ascertain the relationship between toxin sensitivity and the genetic relationship of the assayed Streptococci species, we performed a phylogenetic reconstruction of the relationships among species using the small ribosomal sub-unit (16S) genes, available for these species in GenBank. A Maximum likelihood (ML) reconstruction was performed as implemented in PAML 4.0 [1] from a starting Neighbor Joining (NJ) tree obtained in Mega 4.0 [2]. Support for the nodes was inferred by two procedures: i) bootstrap support for the nodes was estimated in multiple ML searches of the trees in PAML 4.0 and applied a 50% majority consensus rule to assign support to the nodes; ii) Bayesian posterior probabilities as implemented in Mr. Bayes 3.4b [3,4]. Both produced similar results and only the boostrap support is shown in the results. We assayed toxin activity as explained in the materials and methods section of the manuscript, and performed a parsimony reconstruction of toxin sensitivity with Mesquite [5].
The relationships inferred show that S. mitis and S. oralis are the closest related species to S. pneumoniae within the group of Streptococci species considered for the assays (see Figure S6). The colors of the branches correspond to the activity of the supernatant on each one of the terminal species as well as the Maximum parsimony reconstruction of the activity in the deeper nodes (performed in Mesquite, Madison and Maddison 2007). In red it is shown the branches on which activity was detected or inferred, and in black the branches were no activivity was detected or inferred.