Population Dynamics Model and Analysis for Bacteria Transformation and Conjugation

We present a two-species population model in a well-mixed environment where the dynamics involves, in addition to birth and death, changes due to environmental factors and inter-species interactions. The novel dynamical components are motivated by two common mechanisms for developing antibiotic resistance in bacteria: plasmid {\it transformation}, where external genetic material in the form of a plasmid is transferred inside a host cell; and {\it conjugation} by which one cell transfers genetic material to another by direct cell-to-cell contact. Through analytical and numerical methods, we identify the effects of transformation and conjugation individually. With transformation only, the two-species system will evolve towards one species' extinction, or a stable co-existence in the long-time limit. With conjugation only, we discover interesting oscillations for the system. Further, we quantify the combined effects of transformation and conjugation, and chart the regimes of stable co-existence, a result with ecological implications.

We present a two-species population model in a well-mixed environment where the dynamics involves, in addition to birth and death, changes due to environmental factors and inter-species interactions. The novel dynamical components are motivated by two common mechanisms for developing antibiotic resistance in bacteria: plasmid transformation, where external genetic material in the form of a plasmid is transferred inside a host cell; and conjugation by which one cell transfers genetic material to another by direct cell-to-cell contact. Through analytical and numerical methods, we identify the effects of transformation and conjugation individually. With transformation only, the two-species system will evolve towards one species' extinction, or a stable co-existence in the long-time limit. With conjugation only, we discover interesting oscillations for the system. Further, we quantify the combined effects of transformation and conjugation, and chart the regimes of stable co-existence, a result with ecological implications.  [14] 5 R δ −→ 0 + P R lysis, releasing a copy of P estimated in model

I. INTRODUCTION
The synthesis of large varieties of antibiotics is one of the most important medical inventions in human history. However, the emergence of antibiotic-resistant bacterial population, infamously known as the "superbugs," poses an alarming challenge to public health authorities [1,2], whose efforts in combating resistance are outpaced by the rapid development of new resistance in the bacteria cells [3][4][5]. Bacteria become resistant to antimicrobial agents as a result of exchanging genetic materials. Most of the bacterial DNA is contained in the chromosome, which provides the genetic identification for each cell. In addition to the chromosome, there are plasmids -small circularized DNA molecules independent of chromosomes -that have essential genes for plasmid functions and accessory genes [6,7]. It is the accessory genes that can confer antibiotic resistance to the host bacteria. There are two main mechanisms of horizontal gene transfer (HGT) through which an individual bacterium cell acquires antibiotic resistance from plasmids: transformation, where a cell incorporates a plasmid from its surroundings, and conjugation, where a plasmidcarrying cell transfers the plasmid (a single strand of the plasmid DNA) to a plasmid-free one by direct cell contact [8]. A schematic of the two processes is shown in Fig.1. The plasmids are transcribed and translated by the cellular machinery of the host cell and thus propagate in the population through cell division. Due to the short cell doubling cycle of bacteria, a small initial fraction of plasmid-carrying cells will be rapidly amplified in the population through division. There is an obvious advantage for plasmid-carrying cells in the presence of antibiotics. In the absence of antibiotics, however, the plasmid-carrying cells can still quickly outnumber the plasmid-free ones as evidenced by both laboratory and clinical findings [15][16][17]. The domination of either genotype depends on the respective doubling time as well as the rates of transformation and conjugation, which vary according to the plasmid concentration, the medium and the type of bacteria population. It is therefore crucial to quantitatively chart the parameter space and identify how the two populations compete under different conditions, which will bring further insights into the emergence of the resistance-dominant population.
The dynamics of biological systems has been extensively studied for decades using nonlinear differential equations such as the venerable Lotka-Volterra equation [18,19], agent-based modeling [20,21], and numerous tools developed in statistical physics [22,23]. Such interdisciplinary efforts have brought deeper insights into understanding biolog-ical systems of different scales, ranging from neuron networks to global pandemic spread, as well as expanding the understanding in generic complex systems. In this study, we take the biological processes of transformation and conjugation among bacteria cells as motivation, and construct a theoretical model for two growing bacterial populations of different genotypes inhabiting a well-mixed environment (e.g. a shaken flask): the plasmid-carrying cells that are antibiotic-resistant (R), and the plasmid-free ones (S) that are susceptible to antibiotics and can be converted to R through either transformation or conjugation. The competition between the two sub-populations depends on the dynamics of the two HGT mechanisms and their respective doubling time, which reflects the plasmid carriage cost with varying benefit according to environmental conditions such as the level of antibiotics. By using a combination of deterministic calculation and numerical methods of solving coupled nonlinear differential equations, we examine the effects that these parameters have on the overall population in Sections II and III. One of the main questions we attempt to answer in this study is with the interplay of transformation and conjugation in various growth regimes, how the parameter space is configured such that the overall system displays one species' extinction (or "fixation" by the other, which is commonly used when the changes are genotypic) or coexistence with one species dominating. We summarize and discuss the biological relevance in Section IV.

II. POPULATION DYNAMICS MODEL
We formulate a model at the population level that incorporates growth, death and the mechanisms of HGT between the S and R populations in a well-mixed culture. Initially, the system is seeded with S 0 susceptible and R 0 resistant cells, which have growth rates b S and b R . The microscopic growth rates of S and R are accessible through experiments by measuring the population doubling time µ S/R of the two genotypes: We focus on two HGT mechanisms through which cells develop antibiotic resistance: transformation and conjugation. When there is a large number of plasmids (P ) in the environment, the transformation rate depends on the cell's innate ability to take up extracellular DNA. Therefore we choose a constant transformation rate α(P ) = α 0 to model the unlimited plasmid supply or that the plasmid has a very high affinity to the cell. Accounting for the availability of free plasmids P in the environment and the kinetics of plasmid uptake, we also study the scenario where the transformation rate depends linearly on P when the concentration of plasmids is low and transitions to the constant rate as P increases, namely α(P ) = α 0 P/(P + K P ). Here K P reflects the plasmid affinity to the cell and is the concentration of plasmids at which transformation occurs at the half-maximum rate. This form is analogous to the Michaelis-Menten kinetics [24] in enzymatic activities. Other P -dependence scenarios can be studied with a similar approach.
During the process of conjugation, the double stranded plasmid DNA in the donor cell becomes two single strands, one of which is transferred from the donor to the recipient through direct cell-to-cell contact mediated by a tubelike structure called "pilus" that pulls the recipient and perforates it [14,25]. Once the single strand of plasmid DNA is in the recipient, the DNA replication process restores the complementary strand and both the donor and the recipient possess the entire plasmid DNA. We model the conjugation process with constant rate γ in the well-mixed population.
Finally, R cells lyse 1 with rate δ. Unlike the phenotypical switching where a bacterium changes from a normal cell to a persistent cell in the presence of antibiotics with the same genetic content [26], a resistant cell dies when its cell wall breaks down. An R-cell is then removed from the population upon lysis and the plasmid it carries is released back to the environment. In the plasmid-limiting scenario, the recycling of plasmids proves to be significant in maintaining a resistant population. When there are antibiotics of various concentrations in the environment, the lysis rates of S and R populations lead to more interesting dynamics which we will not focus on in this article.
The above reactions are summarized in Table I. By using a combination of theoretical analysis for insights and numerical exploration for a wider range of parameters, we aim to provide a comprehensive picture of the roles transformation and conjugation play. To connect our theoretical model to experiments, we also include typical values of the parameters, such as b S , b R , α and γ, in our model. However, the cell lysis process is usually harder to ascertain and we will estimate the parameter δ in our model. More details and discussions will be included in Section III.
The dynamics of the two populations, S(t) and R(t), can be described by the following set of equations after properly non-dimensionalize the parameters without loss of generality: In addition, the plasmids in the environment vary due to the uptake by S cells, as well as the release by R cells upon their lysis: To analyze the above set of coupled differential equations, we used a combination of theoretical analysis and numerical solvers. For the latter, we coded the equations using Python (v3.8) [27]. We used Euler's method [28] and the ODE solver odtint from Scipy (v1.5) [29] to perform the numerical simulations, which yielded consistent results. The numerical results presented below are from Euler's method.
Before delving into a systematic analysis of the full model, we analyze the results of two limiting scenarios which provide guidance on the choice of parameters in the later simulations.
A. Transformation with constant α = α0 We first consider the transformation-only case, namely γ = 0, and the transformation rate is a constant α 0 . This corresponds to the case where the environment concentration of plasmids is very high, or K P is very small. With the effective growth rateb S ≡ b S − α 0 , the dynamics of the S population is an exponential function Similarly, the growth of the R population involves an exponential term withb R ≡ b R − δ and the influx of the transformed S cells: The ratio between the two cell populations f ≡ R/S is given by: where ε ≡b S −b R , the difference between the two effective growth rates. When ε < 0, namely the effective growth rate of R is greater than that of S, then R population fixates, or f (t) = ∞, in the long-time limit with a characteristic time scale of |1/ε|, shown in Fig.2. For the case of ε = 0, f (t) = R 0 /S 0 + α 0 t, indicating that the dominance of R increases linearly with t, slower than the ε < 0 case above. Still in the long-time limit, 1/f (t → ∞) = 0, leading to R-fixation. Fig.3 shows the linear increase in f (t), with slopes given by α 0 . Note the different scales of f (t) in Figs. 2 and 3.
When ε > 0, the system displays an S-R coexistence as f (t → ∞) = α 0 /ε. Either S or R can be the majority in the entire population depending on the competition between transformation and growth differentials: In the case of R-dominance, α 0 > ε. Together with the S-R coexistence condition ε > 0, this gives And similarly, S-dominance is reached when α 0 ∈ 0, (b S −b R )/2 . In Fig.4, we show two cases where R is the majority in steady state (dotted and dashed lines), and one where S dominates (solid line).

B. Conjugation with constant γ
Having seen that transformation-only leads to either R-fixation or stable coexistence, we now turn to the complimentary case where α(P ) = 0 and conjugation is the only HGT mechanism. Eqns.(1) and (2) contain a nonlinear term, similar to the Lotka-Volterra type interaction [18,19]. We find two fixed points: (S * 1 , R * 1 ) = (0, 0) and (S * 2 , R * 2 ) = (−b R /γ, b S /γ). For there to be a coexistence regime where both S and R are non-negative,b R must be negative. In this case, (S * 1 , R * 1 ) is a saddle point which S evolving away from and R into. Furthermore, from Eqns. (1) and (2) we find the trace and the determinant of the Jacobian matrix at (S * 2 , R * 2 ) to be 0 and −b SbR , respectively. Using the standard stability analysis [30], we find that (S * 2 , R * 2 ) is a neutrally stable center.
In this case, the ratio f * ≡ R * 2 /S * 2 = −b S /b R does not depend on the conjugation γ. What maybe slightly counterintuitive is that for f * < 1, i.e. a center with more S-cells than R, b S is constricted to be in (0, −b R ), rather than a greater value. We will return to this point in the next section when we consider the combined effect of transformation and conjugation to the overall population.
To determine the individual closed orbits, we first eliminate t by dividing dS/dt by dR/dt and then separate the variables: where C is the constant of integration. We now have an expression to describe the trajectory for the conjugation-only system: We show in Fig.5 the phase portraits of the two sub-populations forb R < 0. Setting b S = 1 and varying the ratio between −b R and γ using the range of conjugation rates in Table I, we see the competition between growth and conjugation leads to different trajectories, all around the center (S * 2 , R * 2 ), indicated as a red dot. In Fig.6 we observe the oscillations in S and R on one of the closed orbits. The fact that S-R can coexist wheñ b R < 0 is particularly interesting because the increase in R-cells is through converting more S-cells by conjugation. A closer comparison between conjugation-only and conjugation with transformation is provided in the following section.

III. COMBINED EFFECTS OF TRANSFORMATION AND CONJUGATION
Informed by the above discussions on the individual effects of transformation and conjugation, we study the model of which the full dynamics is described in Eqns.(1)-(3) with a range of parameters to chart the different regimes of the sub-populations of S and R cells, with special attention on the S-R coexistence conditions. From the discussions in Section II B, we learned that a physical coexistence of S and R cells requires b R < δ in the conjugation-only case. With transformation included, first let's consider the case of constant transformation, or K P P . Thus α(P ) → α 0 and the new non-trivial fixed points are: The fixed point at (0, 0) remains a saddle point. It is worth noting that the value f * is the same as in the conjugation-only case. Meanwhile, the overall population of the system approachesb Combining transformation and conjugation brings new behaviors in the overall system. Unlike the conjugationonly case where S and R evolve along a closed orbit centered around (−b R /γ, b S /γ) with temporal oscillations, as shown in Figs.5-6, once transformation is introduced, these orbits turn into spirals. We find an S-R coexistence wheñ b R < 0 andb S > 0. To determine the stability of (S * , R * ), we can again examine the Jacobian at this point: the trace and the determinant isb R α b S < 0 and −b SbR > 0. Thus (S * , R * ) will be a stable node or a stable spiral. However, for it to be a stable node,b R > −4(1 +b S /α) 2b S . This means whenb R is almost one order of magnitude greater thañ b S , (S * , R * ) is a stable node. In other cases whereb R ∼b S , the system has a stable spiral towards (S * , R * ).
In Fig.7, the system starts with S 0 = R 0 = 1 and the initial oscillations of the two sub-populations quickly settles into (S * , R * ) = (4. 5,9), as given in Eq. (6). Their trajectory is shown in Fig.8, with the steady state being the red dot at (S * , R * ). We also show the closed orbit for the same set of parameters except for α 0 = 0 as a reference. As predicted, the S-R coexists with a stable ratio of f * = 2.0 in the long-time limit. Noticeably, the spiral is smaller than the conjugation-only orbit due to the damped oscillations in S and R when transformation is introduced (α 0 = 0.1 in this case). The general shape of the spiral will depend on the initial conditions of S 0 and R 0 , which does not affect the final coexistence ratio f * .  As alluded to in Section II B, for S to be the majority in the overall population in steady state, the growth rate b S needs to be less than −b R . For example, when b S = 0.4 while −b R = 0.5, as in Fig.9, the final S-R coexistence ratio f * = 0.8. When b S < −b R , even though S is doubling at a lower rate, that also leads to fewer S being conjugated or transformed into R. The final coexistence state contains more S than R cells. Next let's turn to the case where transformation rate is no longer a constant. Here two parameters are at play: the initial concentration of plasmids in the environment, P 0 , and the transformation kinetics K P . The former obviously depends on the condition of the growth medium, while the latter is a biochemical parameter where a higher K P means a lower plasmid-bacteria affinity. The dynamics of the plasmids in this case is given in Eq. (3).
If the lysis rate δ is sufficiently large, or R is the dominating species, then the plasmid concentration will eventually build up and it becomes similar to the α(P ) = α 0 case as discussed above. However, the transient behavior in this regime contains some interesting details. For an environment with few low-affinity plasmids, K P P , the transformation rate is: The stability condition with conjugation-only case also requires δ > b R while transformation with the low-affinity plasmids is slow. The buildup of plasmids -consequently α(P ) -therefore sees a "burst"-like behavior, shown in Fig.10, where we look at low-affinity plasmids with K P = 10 3 P 0 . The bursts coincide with the decrease of R due to lysis and increase the fastest when S reaches its minimum so that dP/dt is the greatest. Additionally, the burst cycle tracks the S-R oscillations.

IV. SUMMARY AND DISCUSSION
We studied a two-species population model aimed to explore the interplay between bacterial transformation and conjugation in addition to their natural growth cycles. Using a combination of theoretical and numerical analysis, we analyzed the complete effects of transformation: the overall population will either be fixated by R-cells or in stable co-existence. In the latter case, the majority cell type is determined by the transformation rate α 0 and ε, the difference between the effective growth rates of the two sub-populations. The conjugation-only case gives R-fixation, unless when the effective growth rate of R-cellsb R is less than 0, in which case we can actually observe a stable orbit for S and R, with each sub-population oscillating. This could potentially be interesting in the biologically-relevant context: carrying an additional plasmid confers growth disadvantage for R-cells asb R < 0, they can nevertheless survive steadily only by conjugating more S to R. More interestingly, when transformation and conjugation are both considered, the orbit seen in the conjugation-only case turns into a spiral with a stable center, indicating a steady state where the final ratio between the two sub-populations f * reaches a constant.
In our numerical analysis of the model, the reaction rates tested are not directly comparable to the experimentally obtained ones listed in Table I because the equations of the system were first non-dimensionalized. However, our results point to the important regimes where interesting behaviors, such as stable co-existence and oscillation, emerge when the relative reaction rates are taken into consideration. In particular, the difference between the two effective growth rates ε in the case of transformation, and the ratio between the two effective growth rates characterized by f * = −b S /b R in the case of conjugation are direct indication on the overall population behavior.
From a modeling perspective, our findings are novel in that we quantitatively analyze the effects of two main horizontal gene transfer (HGT) mechanisms at a population level. The results are significant especially with the use of a minimal model with few parameters, mostly experimentally accessible. The lysis rate δ is the only estimated parameter here, due to its elusiveness because cell lysis often involves various other environmental and physiological factors. And yet it is crucial in determining whether the system has a steady state, or which sub-population dominates. We thus continue to seek experimental breakthroughs on further clarification on the lysis process.
The results from this study are potentially helpful to a wider range of explorations. In our study, there is no limit on the population growth aside from the dynamical parameters. It is conceivable to impose an environmental constraint such as a carrying capacity and analyze the difference. For this study, the effects of antibiotic resistance is binary, in the form of carrying the plasmid, and deterministic. It is possible for cells to lose a plasmid without lysing (e.g. [31]), in which case there needs to be an additional reaction to capture R −→ S + P process. And in terms of stochasticity, our preliminary study using kinetic Monte Carlo shows consistent results with what we present here, and is more powerful to generalize to more reactions with other complications. Further extensions to include the population spatial structure beyond the existing well-mixed population are also being considered. The spatial aspect is especially compelling to investigate because both transformation and conjugation processes depend on local availability of free plasmids as well as direct contacts between S and R cells. Thus the diffusion of plasmids as well as the cells are expected to be relevant, an aspect not necessary for the well-mixed case.