Modeling sRNA-Regulated Plasmid Maintenance

We study a theoretical model for the toxin-antitoxin (hok/sok) mechanism for plasmid maintenance in bacteria. Toxin-antitoxin systems enforce the maintenance of a plasmid through post-segregational killing of cells that have lost the plasmid. Key to their function is the tight regulation of expression of a protein toxin by an sRNA antitoxin. Here, we focus on the nonlinear nature of the regulatory circuit dynamics of the toxin-antitoxin mechanism. The mechanism relies on a transient increase in protein concentration rather than on the steady state of the genetic circuit. Through a systematic analysis of the parameter dependence of this transient increase, we confirm some known design features of this system and identify new ones: for an efficient toxin-antitoxin mechanism, the synthesis rate of the toxin’s mRNA template should be lower that of the sRNA antitoxin, the mRNA template should be more stable than the sRNA antitoxin, and the mRNA-sRNA complex should be more stable than the sRNA antitoxin. Moreover, a short half-life of the protein toxin is also beneficial to the function of the toxin-antitoxin system. In addition, we study a therapeutic scenario in which a competitor mRNA is introduced to sequester the sRNA antitoxin, causing the toxic protein to be expressed.


Introduction
Small regulatory RNA (sRNA) plays an important role in gene regulation in organisms from bacteria to mammals by controlling for example translation and/or mRNA stability and the list of known RNA regulation systems keeps increasing at a rapid pace [1][2][3]. sRNA regulation possesses characteristics that are distinct from protein regulation, in particular a threshold-linear response, which provides an ultrasensitive mechanism for regulatory switching, and the possibility of hierarchical crosstalk, which allows prioritizing of expression [4].
Some of the best known sRNA regulation systems are related to plasmid replication and plasmid maintenance in bacteria. In the first case, an sRNA controls whether synthesis of a replication primer proceeds to plasmid replication [5][6][7]. In the second, the plasmid encodes a (type I) toxin-antitoxin system such as the hok/sok system ("host-killing/suppression-of-killing"), encoding a protein toxin and an antisense RNA which acts as an antitoxin by binding to a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 the toxin mRNA and blocking ribosome access, thus preventing toxin synthesis [8]. (Other types of toxin-antitoxin systems have different functions [9], in particular related to the formation of persister cells [10].) The toxin-antitoxin system, which is also known as an "addiction module", maintains the plasmid number through post-segregational killing of plasmid-free progeny due to differential stability of the toxin and antitoxin RNAs. The killing is done by a potent protein toxin that irreversibly damages the cell membrane [11]. In a steady state with a stable plasmid concentration, sRNA antitoxin exists in excessive molar amount compared to target mRNA, such that the latter is entirely sequestered in translationally inactive sRNA-mRNA complexes [12]. However, when a progeny cell becomes plasmid-free through cell division, synthesis of both toxin mRNA and antitoxin sRNA are stopped. The sRNA, which has a very short half-life, is rapidly depleted and the more stable target mRNA can be translated into toxic protein, killing the cell (Fig 1).
To understand the design of the gene circuit encoding the toxin-antitoxin mechanism of post-segregational killing, here we analyze a theoretical model for the dynamics of such a circuit. The model is related to and extends previous models for sRNA-based post-transcriptional regulation [4,13,14]. It allows us to address the parameter dependence of the genetic circuit to identify essential features and criteria for its efficient function, for example, whether there are any criteria beyond the difference between toxin and antitoxin RNA lifetime. It also allows us to test whether the system can be designed in such a way that only complete loss of the plasmids triggers killing and not a low but non-zero plasmid copy number. We also use the model to study a proposed antibacterial strategy [15,16] making use of regulatory crosstalk between RNAs. To answer these questions we performed extensive parameter sweeps, scanning all parameters of the model individually as well as extensively testing randomly chosen parameter sets, thus varying all parameters simultaneously. Overall, the model confirms the known principles for the function of the toxin-antitoxin mechanism as the dominant ones, in particular, the difference in RNA lifetime and the threshold condition on the synthesis rates. In addition, it also indicates a few new design features. Specifically, our analysis of the parameter dependence shows that the stability of the mRNA-sRNA complex plays an important role as well and should exceed the stability of the free antitoxin. Moreover, an unstable toxic protein is superior to a stable one for induction of killing upon plasmid loss. Simplified drawing of the sRNA regulated toxin-antitoxin mechanism for plasmid maintenance. Structural transitions and processing of the mRNA have been omitted [8]. Parameters of the system are labeled in blue: α m , α s , α p are the synthesis rates of mRNA, sRNA and protein (transcription and translation rates, respectively); β m , β s , β c , β p are the degradation rates of mRNA, sRNA, the mRNA-sRNA complex and protein, respectively; h + and h − are the binding and unbinding rate of the complex.

Dynamic Equations and Analytical Solution of Steady State Concentrations
The dynamics of a toxin-antitoxin system as described in Fig 1, is modeled with four coupled ordinary differential equations for four dynamical variables, the concentrations of the toxin mRNA (m), the antitoxin sRNA (s), the toxic protein (p) and of the mRNA-sRNA complex (c), in which the mRNA is silenced. All concentrations are expressed in units of number of molecules per volume of a cell. The four equations describe the synthesis and degradation of the mRNA, the sRNA, and the protein, as well as the formation and dissociation of the mRNA-sRNA complex: The 10 parameters of these equations are as follows: α m , α s , α p are the synthesis rates of mRNA, sRNA and protein (transcription and translation rates, respectively); β m , β s , β c , β p are the degradation rates of mRNA, sRNA, the mRNA-sRNA complex and protein, respectively; h + and h − are the binding and unbinding rate of the complex, and g is the plasmid copy number per cell volume.
A few remarks are in order with respect to the model's generality. First, for many toxinantitoxin systems, including hok/sok and the pAD1 par system, there exist intramolecular structures of the toxin mRNA that prohibit the translation of protein by partially restricting ribosome access. Such structural states of the mRNA are not included in our model. However, averaging over the active and inactivated configurations of the mRNA, we obtain an effective synthesis rate of the protein and an effective degradation rate of the mRNA, which do not qualitatively alter the model. Averaging requires however that the mRNA switches between those configurations sufficiently rapidly. Second, saturation of the toxin-antitoxin binding naturally occurs within this model (provided that binding is sufficiently strong) and is accounted for in the nonlinear binding term h + Á m Á s. As a consequence an excess of antitoxin sRNA over toxin mRNA results in all toxin mRNA being sequestered in complexes, as discussed in more detail below. Third, the effect of dilution by cell growth can be seen as included in the degradation rates. Dilution by cell growth sets a lower limit on all degradation rates, so proteins or RNA effectively cannot have half-lives exceeding the doubling time of the cell (even though their true degradation rate, due to proteolysis and RNA degradation could be lower).
The steady state solution is obtained by setting the time derivatives on the left hand side of the equations to zero. The steady state for this system can be explicitly given, because the nonlinear terms cancel each other in Eqs (1a) and (1b), leaving a quadratic equation which can be solved analytically. This leads to the steady state sRNA concentration: The steady state concentrations of the other three components can be in turn given as: This result includes some limiting cases that have been studied in earlier work on sRNAdependent gene regulation. In the limit of large binding and unbinding rates h + and h − (i.e. the limit of "rapid equilibrium"), the steady state concentration of the RNA complex concentration becomes c Ã = (m Ã s Ã h + )/h − , as previously obtained by Legewie et al. [13]. The limit of irreversible binding, h − = 0, in which only the equations for sRNA and mRNA concentrations need to be considered, was previously studied by Levine et al. [4] and Mitarai et al. [14].

Numerical Methods
To study the dynamic behavior, Eqs (1a-1d) are numerically integrated using the 4th order Runge-Kutta method with an integration time step on the order of 1 × 10 −4 min, for a time span of 300 min. We note that the time unit of the dynamics could be made dimensionless by rescaling all rates relative to one rate that determines the time unit. Correspondingly, the results presented below will typically depend on ratios of time scales or rates. To achieve both stability and speed, an adaptive time step for the numerical integration is implemented, such that when the integration becomes numerically unstable, a smaller time step is used. A tell-tale sign for numerical instability is the appearance of negative concentration values of one or more components. The analytical results for the steady state concentrations are in good agreement (up to floating point error) with the results obtained from the numerical integration.

Results and Discussion
Qualitative Description of the Dynamics of the Genetic Circuit For the model described above, the following scenario is considered: the dynamics of the toxin-antitoxin system begin with a cell which has just acquired one or more copies of the plasmid, but has not yet synthesized any of its products, i.e. we start with m = s = c = p = 0. The dynamics are numerically integrated for a sufficiently long time, so that a steady state of the system is reached. At a certain time point (here at t = 150 min of a total t = 300 min), the plasmid copy number is set to 0 to mimic plasmid loss, e.g. due to insufficient replication or unequal partitioning of the plasmid during cell division. To illustrate the dynamics, realistic values of the parameters are used, which are estimated from the literature ( Table 1). Eqs (1a-1d) are numerically integrated to produce Fig 2. In Fig 2, all four concentrations increase initially until the first steady state is reached. At t = 150 min, due to the loss of all plasmid copies, the synthesis of new RNA molecules stops. Nevertheless, a transient increase in the free mRNA concentration is observed, which results in a transient increase of the protein concentration. The transient increase in mRNA is a result of sRNA degrading much faster than mRNA, so when a mRNA-sRNA complex dissociates, there is a surplus amount of free mRNA molecules released, which in turn allows the synthesis of toxic protein. Thus, the main function of the toxin-antitoxin system, plasmid maintenance via the synthesis of a toxin upon plasmid loss, which results in the removal of plasmid-free cells from the population, is dependent on a transient dynamics rather than on a steady state. This behavior is in contrast to other sRNA-based regulation systems, which are based on similar mechanisms, but control the steady state concentration of mRNA [4].
For a quantitative characterization of the dynamics, specifically that of the protein concentration change after the loss of plasmids, two quantities are measured for each simulated scenario: (i) the peak fold-increase of the protein concentration, R, defined as the ratio of the maximum of the protein concentration after plasmid loss and the steady state protein concentration before the plasmid loss, R = p peak /p Ã ; and (ii) the transient peak width T p , defined as the width of the peak at half maximum. Table 1. Table of parameter values based on the literature for toxin-antitoxin systems and related sRNA regulation systems. For some parameters, the values are know to depend on growth conditions, which is not explicitly incorporated into our model. For components with long half life, the degradation rate is dominated by dilution through cell growth and division, which provides a lower limit on the possible (effective) degradation rates.
Range of values 0.  The effectiveness of the toxin-antitoxin mechanism is partially determined by the protein concentration peak fold-increase R, as opposed to by the absolute concentration of the toxin. A high peak fold-increase allows a toxicity threshold to be set, such that the fluctuations of the steady state protein concentration are very unlikely to cross such a threshold and thereby "accidentally" lead to cell death, while simultaneously permitting the transient increase of protein concentration when the plasmid number is not being properly maintained to easily overcome such a threshold. For example, if the cell will be killed by a 10-fold-increase in toxin concentration (10 being the killing threshold), and if there is little chance of an accidental 10-fold protein concentration increase due to other factors and noise, the dynamics which can produce R > 10 are considered effective in exerting a robust control on the host cell.

Dependence of the protein dynamics on individual model parameters
We first investigate the parameter dependence of the model by varying each parameter in isolation, holding the other 9 fixed at constant values. In all cases, the peak fold-increase R and the peak width T p are determined as functions of the modulated parameters.
First, the dependence of R and T p on the plasmid copy number g (Fig 3) is investigated. An increase of g is equivalent to increasing the synthesis rates of the sRNA and mRNA molecules by the same ratio. R and T p are both positively affected when g is increased. For the default parameters (Fig 3(a)) as well as for the rapid equilibrium case (Fig 3(b)), the rates of change of both R and T p (i.e. the slope of the curves) are higher when g is low, and lower when g is high.
In the case of rapid equilibrium, i.e. when both h + and h − are high, the changes in R and T p at the loss of one or two plasmid copies are more substantial than in the cases of lower binding and unbinding rates. This is to say, in the case of rapid equilibrium, the model could be effective in triggering killing even when the plasmids are not completely lost.
Another way to think about plasmid maintenance is shown in the simulation of sequential loss of plasmid copies under default parameters (Fig 3(c)). The change in plasmid number does not result in the sudden release of large amount of toxin until the last plasmid is lost. Comparing to when all plasmids are lost at once (Fig 2), sequential plasmid loss will result in a less prominent transient peak when the last plasmid is lost. This observation is consistent with the function of the toxin-antitoxin system to enforce plasmid maintenance, i.e. the host cells are forced to retain at least one copy of the plasmid, but there is little dependence on whether there are more or fewer copies.
The dependence of R and T p on the other 9 parameters are shown in Fig 4. In general, two conditions for the existence of a transient peak of the toxic protein concentration can be extracted from these parameter dependence: First, the synthesis rate of sRNA (α s ) must exceed the synthesis rate of mRNA (α m ). This threshold condition can be observed by comparing the varying values to the default values for mRNA and sRNA synthesis rates. In Fig 4(a), when the synthesis rate of mRNA α m exceeds the default synthesis rate of sRNA α s at 6/min/gene copy, the peak disappears. Conversely in Fig 4(b), when α s becomes larger than α m at a default value of 1/min/gene copy, the protein concentration peak starts to appear (R > 1).
This threshold condition is well-known for sRNA regulation [4]. It can be explained as follows. Increasing the number of mRNA molecules while keeping the number of sRNA molecules constant increases both the protein concentration peak p peak as well as the first steady state protein concentration p Ã . However, p Ã increases faster than p peak as α m increases. Therefore the relative increase of the protein concentration R diminishes as α m increases. A high level of free mRNA before plasmid loss will also mean that plasmid-containing cell might "poison itself" without any loss of plasmids. When α m > α s , there will be more free mRNA than what can be "neutralized" by complex formation, and the toxic protein will be expressed at a high level before the loss of plasmids. The toxin-antitoxin mechanism relies on the low amount of surplus of free mRNA after the loss of plasmids, so that when α m is higher than α s , this effect is lost.
Second, the degradation rate of sRNA must exceed the degradation rate of mRNA: β s > β m . This condition ensures that after the loss of plasmids, a pool of free mRNA builds up, since the sRNA is degraded more rapidly. When sRNA is more long-lived than mRNA, that is, when β m exceeds β s at 1/min, or when β s is lower than β m at 0.2/min (Fig 4c and 4d), a pronounced protein concentration peak is not observed. In the following section, when we discuss the results of random sampling of the parameter space, we will show that this condition on the degradation rates of the RNAs is not as strict as the threshold condition on the synthesis rates. This can be seen in Fig 5(a), where some parameter combination exhibit a protein peak despite a ratio β m /β s > 1, while there never is a peak for α m /α s > 1. Besides these main two conditions, high binding rate, low unbinding rate and sufficient stability of the complex are also needed as shown in Fig 4e-4g. This is consistent with the previous knowledge that the binding rate strongly influences the effectiveness of the toxin-antitoxin mechanism [24,25]. Finally in Fig 4(i), when the protein degradation rate β p is very low (equivalent to having a half life of longer than 5 min), R decreases drastically and eventually drops to a constant around 4, where no protein is being degraded. Thus, the proteolysis of the toxic protein contributes to the efficiency of the toxin-antitoxin mechanism for plasmid maintenance. The translation (protein synthesis) rate, on the other hand, has no effect on the ratio R as expected (Fig 4(h)).
The parameter dependence study shows two features of the toxin-antitoxin system that were not known or have not been demonstrated via theoretical studies. Firstly, it is known that the antitoxin sRNA in plasmid number maintenance is very short-lived [26] (which is indeed true for many types of sRNA [27,28]), while the toxin mRNA has an unusually long half-life [26]. However, in the numerical study above we find when the degradation rate of the sRNA exceeds that of the template mRNA, in some range, the mechanism becomes actually less effective: R decreases in Fig 4(d) after β s becomes larger than around 0.5. However, we did not see this effect clearly in the random sampling of the parameter space in the following section. It is therefore possible this effect only applies to situations where some other parameters, for example the synthesis rates of the RNA molecules, take on certain values.  Table 1 and as in Fig 2) are marked by red lines in the plots.  The color gradient shows the magnitude of the protein peak fold-increase R. In  Fig (a), the results are projected onto α m /α s vs. β m /β s on a log-log scale. In (b) the simulation results that fall within the subregion α m /α s < 0.8 and β m /β s < 4 (marked by lines in Fig (a)) are projected onto β c vs. β s . In (c,d) the parameter sets are further restricted to the subregion α m /α s > 0.8, β m /β s > 4 and β c /β s < 2/3 (marked by a dark line in Fig (b)) and the results are projected onto h + vs. h − in Fig (c), and β p vs. β m in Fig (d), respectively. The color scale values are taken to be log 1.4 R for ease of viewing. Dark lines mark the rough division between the region with more high peaks and the region with fewer. In (c) and (d) the dark lines mark h + = h − and β p = β m respectively. Secondly, in some well-studied cases, the binding of mRNA with sRNA leads to the rapid degradation of the complex [29]. However, this cannot be the case here. As shown in Fig 4(g), when the RNA complex is degraded at a rate higher than 0.8 there is no transient toxin peak at all (R = 0). This can be explained by the fact that mRNA must be released from the complex upon plasmid loss, which requires a sufficiently high concentration of the complex. If the complex degradation rate β c is too high, too few mRNA templates will be left for the translation of the toxin. Nevertheless, complex stability is an important parameter beyond this obvious feature, as will be shown by the second parameter dependence study conducted in the following section "Random Sampling of the Parameter Space" Fig 5(b). Indeed, a stabilization of RNAs in the complex has been demonstrated in at least one toxin-antitoxin system, the par system in the E. faecalis [30].

Random Sampling of the Parameter Space
We have seen that the dependence of the toxin-antitoxin circuit on the individual parameters exhibit expected but also surprising behaviors. However, due to the fact that in our single parameter scans only one parameter was varied at a time, while the other 9 are fixed, the generality of our conclusions is up for debate. To find the region of the 10-dimensional parameter space where the toxin-antitoxin mechanism works effectively, and to explore the joint conditions on the parameters, we randomly sample the parameter space. Random scans of multidimensional parameter spaces have previously been used in various contexts including cytoskeletal dynamics [31] and signaling cascades [32].
We run a total of 4025 simulations of the toxin-antitoxin dynamics with randomly generated parameters. The only restrictions imposed on the random parameter sampling are α m > β m and α s > β s , to make sure that average molecule numbers exceed 1 and the differential equations approach is appropriate. The algorithm generates random values for each of the 9 parameters (g = 6 remains fixed) within the range given in Table 2 and with the aforementioned restriction, at an appropriate sampling resolution. The resolution for α m , β m , α s , β s are chosen such that the log of the ratios α m /α s and β m /β s are uniformly distributed. The sampling is therefore random and to a large degree uniform, or log-uniform.
The value of the protein concentration peak fold-increase R for each simulation is recorded and projected onto the space spanned by α m /α s and β m /β s (Fig 5(a)). The observations obtained by random sampling are consistent with the results of the single parameter dependence study. The results confirm that α m < α s is indeed a strong criterion for the existence of a high protein concentration peak. Fig 5(a) shows when α m /α s > 0.8 there is almost no protein concentration peak after the loss of plasmids. β m < β s is also a very important criteria, but exceptions are allowed: The protein level decreases on average when β m > β s , nevertheless, there are distinct protein concentration peaks beyond β m /β s > 1. However, the majority of the simulations with β m /β s > 1 exhibit no peaks. These exceptions could be due to some other competing effects from the dynamics of the mRNA-sRNA complex or from the dynamics of the RNA molecules.
Next, we focus on the cases that satisfy α m /α s < 0.8 and β m /β s < 4, and plot those as a function of the degradation rate of the complex β c and that of the sRNA β s . This will eliminate toxin-antitoxin dynamics with parameter combinations that are inefficient in producing a transient peak due to the criteria discussed above, allowing new effects to be discovered more easily. Fig 5(b) shows that no peaks occur for β c ≳ β s , in contrast to the rest of the region. In fact, the peaks are visually prominent only when 1.5 Ã β c ≲ β s . A likely explanation is as follows: when β c is low, there is a higher level of total mRNA molecules in bound form both before and after the plasmid loss. Because the DNA gene copies are removed when the plasmids are lost, the transient increase of protein synthesis is almost entirely due to the mRNA released from the complex. Therefore, the complex degradation (as opposed to mRNA degradation) is the dominant reason for total mRNA loss. At the same time, lower β c does not increase mRNA expression before plasmid loss, but merely increases the amount of bound mRNA molecules, and will therefore not lead the plasmid-containing cell to "poison itself". If this loss of mRNA due to complex degradation is slower than the loss of sRNA (the loss of sRNA is equal to the gain of the free mRNA when the plasmids are lost), the toxic protein can be synthesized from the surplus free mRNA. Therefore, if antitoxin sRNA is less stable than the total mRNA, there will be a transient peak of mRNA, and correspondingly one for the toxic protein, after plasmid loss. However, when a higher level of free mRNA is present before plasmid loss, due to faster degradation of sRNA, i.e. when β s is large, R is in turn negatively affected. Considering this negative impact, β c ≲ β s therefore does not always guarantee high peaks, which is shown from a wide range of variations in R within this region, and there is no increase in R when β c is increasingly small compared to β s . It is surprising that the stability of the complex plays such an important role in the toxinantitoxin mechanism. It is secondary only to the threshold condition on the synthesis rates and the differential stability of the RNA molecules. This highlights the crucial role played by the RNA complex in the dynamics. In connection to the single parameter study conducted before, this demonstrates that the two known effects of unstable sRNA [26] and unstable RNA complexes [29] need to be understood in relation to each other and not in isolation. Unstable RNA complexes are only useful for the functioning of the toxin-antitoxin system if the sRNA is more unstable than the RNA complex.
Next, we restrict our set of parameter combinations further, by imposing β c /β s ≲ 2/3, α m /α s < 0.8 and β m /β s < 4. Fig 5(c) shows that for a prominent protein concentration peak to occur, h − /h + should not be larger than one. Projections onto β p and β m show that when β p is close to 0 there is no major peak (Fig 5(d)). This means that by making extremely stable toxic proteins we cannot increase the relative increase in protein concentration. It also shows generally β m < β p result in less pronounced peaks than β m > β p . This is usually satisfied in real situations, since mRNA molecules are typically less stable than proteins.
The fact that high peaks do not occur in the subregion with h − /h + > 1 tells us that the binding rate needs to be higher than the unbinding rate, i.e. the binding needs to be strong. However, it also cannot be infinitely strong, i.e. irreversible, because in that case there will be no free mRNA available for protein synthesis after plasmid loss.
Consistent with the results from the single parameter dependence study, the synthesis rate α p does not have an obvious effect on the protein concentration peak. This shows that the protein concentration peak is mostly a result of the dynamics between the mRNA, sRNA and the RNA complex, characterized for example by relations between β m and β s , α m and α s , β c and β s , and not of a high synthesis rate of the protein. This also emphasizes the importance of understanding sRNA regulation purely from the point of view of RNA dynamics instead of protein dynamics, which is traditionally viewed as playing the dominant role in cellular regulation.
To summarize, by randomly sampling the parameter space as well as performing a single parameter dependence study, we determine that the conditions for an efficient toxin-antitoxin mechanism are as follows: 1. The synthesis rate of mRNA should be lower than the synthesis rate of sRNA; 2. The degradation rate of mRNA should be relatively low, compared to the degradation rate of sRNA; 3. The stability of the complex should be higher than the stability of the sRNA antitoxin; 4. A high affinity and irreversibility (but not complete irreversibility) in RNA complex formation and some degrees of protein instability are also helping factors, but are of less significance.
5. Protein synthesis rate does not contribute to the forming of a transient protein concentration peak.
Despite the nonlinear nature of our system, by individually controlling the available parameters, a genetic circuit could be engineered to produce specific effects, such as a higher increase in toxin concentration after the loss of plasmids for an effective duration (from a few minutes up to an hour).

Analytical approximation to the transient protein concentration peak
Using a simplified version of the differential equations after plasmid loss, an analytical expression for the transient peak can be obtained, which qualitatively describes the dynamics. For this approximation, we assume that after the plasmid copies are lost, all mRNA molecules are immediately available in free form. That is, the number of free mRNA molecules after plasmid loss equals the sum of the mRNA steady state concentration m Ã and the complex steady state concentration c Ã . The differential equations which describe the dynamics are as follows: The initial conditions and integration constants required to solve the differential Eq (3) are provided by the steady state concentrations before and after the loss of plasmids, respectively. Explicit expressions for mRNA and protein concentration as a function of time can then be obtained: where m Ã , c Ã and p Ã are the steady state concentrations for mRNA, complex and protein before plasmid loss, and m Ã 2 and p Ã 2 are the steady state concentrations after plasmid loss (see Section "Dynamic Equations and Analytical Solution of Steady State Concentrations"). In the case that all plasmid copies are lost, the steady state concentrations m Ã 2 and p Ã 2 are both zero and the above expressions can be simplified to: These expressions qualitatively describe the transient peak of the protein concentration after plasmid loss, however the height of the peak is strongly overestimated (Fig 6). In the full model, mRNA is released slowly from the complex, which reduces the peak height. This behavior can be mimicked within the approximation by the sudden release of only a fraction of the mRNAs (replacing m Ã by an effective mRNA concentration somewhere between m Ã and c Ã ).
Despite these quantitative shortcomings of the approximation, it is useful to understand some properties of the dynamics. For example, one can easily see from the solution (Eq (5b)), that the translation rate does not affect the peak fold-increase R, because both terms are proportional to α p , while the steady state concentration before plasmid loss, p Ã , is also proportional to α p . Thus, R being the ratio of the two is independent of α p . This result should also be true for the full dynamics since the analytical approximation used here is only considering the RNA dynamics and assumptions about the protein dynamics have not been made. The dependence on β m − β p is also important: An increase of this quantity is equivalent to an increase of p Ã , consistent with the tendency observed in the random parameter sampling (Fig 5(d)).

Introduction of a Competitor mRNA to Increase Toxin Levels
Regulation by small RNAs typically displays crosstalk with multiple mRNAs under the control of the same sRNA [33]. In the case of toxin-antitoxin systems, this can be exploited as an antibacterial strategy [15,16]. By inducing a gene encoding a competitor mRNA (on the same or another plasmid or on the chromosome) that can bind to the antitoxin sRNA and hence allows toxin mRNA to exist in free form, a similar increase in toxin concentration can be induced as with plasmid loss, which may also lead to the killing of the cell. To describe the dynamics of this scenario, the equations from above are extended to include a second type of mRNA, the competitor, and the corresponding mRNA-sRNA complex. The full dynamics is then described by the following equations: where m 2 and c 2 denote the concentrations of the competitor mRNA and the complex it forms with an sRNA molecule, respectively. k + , k − and β c 2 are the binding rate, unbinding rate and the degradation rate of this new complex.
We assume that the competitor is induced at the time t = 150 min. Then, from t = 0 to 150 min, the dynamics is the same as Eq (1). At t = 150 min, Eq (6) are being integrated, which corresponds to the competitor gene being turned on. With zero initial concentrations for all components (at t = 150 min, new variables r = v = 0), an example simulation is shown in Fig 7. A transient protein concentration peak can be observed with some specific combinations of parameters, but the majority of the simulations we run show no transient peak, as in Fig 7. However, a non-zero steady state toxin concentration exists after induction of the competitor which can be effective in killing the host cell, since the mRNA encoding the toxic protein is still produced from the plasmid. Thus, in contrast to the case considered before, here the toxin-antitoxin system does not perform its key function in a transient dynamical fashion, but rather in the steady state, with the competitor either induced or not induced. Therefore, we define a new fold-increase parameter e R of the toxin concentration as the ratio between the steady states of the protein concentration before and after the synthesis of the competitor RNA.
Similar to the procedure in Section "Dynamic Equations and Analytical Solution of Steady State Concentrations", the new steady state concentration for the mRNA f m Ã 2 turns out to be a positive root to the cubic equation: After algebraic manipulations, it can be written in canonical form as follows.
where A, B, C and D are functions of the rates. An analytical formula for the solutions of a cubic equation can always be given explicitly, but the expressions are excessively lengthy and hence not shown here. With a root finding algorithm such as the bisection method one can easily find the three roots to the cubic equation numerically. Using values similar to the experimental values, we found that even in cases where there is more than one positive root, we can still pick out the correct solution because usually one of the two positive roots corresponds to an unrealistically large concentration.
The analytical solution to Eq (7) is then compared to the numerically generated steady state concentrations after integration, and very good agreement (up to floating point error) is found. Once we obtain the steady state mRNA concentration f m Ã 2 , then following e p Ã 2 ¼ f m Ã 2 Á a p =b p one can solve for the steady state protein concentration. The ratio f p Ã 2 =p Ã , where p Ã is given by Eq (2), gives an analytical value for the fold-increase e R. This is how we could theoretically obtain an analytical formula for e R in the case of competitive RNA binding, which is just as effective in killing the host cell at some triggering event as the original toxin-antitoxin mechanism.  We observe that as soon as the synthesis of the competitor is turned on, i.e. α 2 > 0, there is a substantial fold-increase in the toxin concentration, i.e. e R > 1. As α 2 becomes large, the increase in e R is diminished because all sRNA molecules are bound. Changing the degradation rate of the competitor mRNA β 2 results in small variation in the toxin foldincrease response e R, meaning that for the default parameter combination, there are very few free competitor mRNA molecules and most are bound in complexes. The degradation rate of the competitor mRNA-sRNA complex β c 2 positively affects e R, because it increases the destruction of the antitoxin in bound state, allowing more target mRNA to be in free form available for protein synthesis.
The parameter variation in Fig 8 shows that the synthesis rate of the new competitor mRNA typically plays the dominant role compared to other new parameters, a general feature of cross-talk in sRNA regulation [33].

Conclusion and Summary
We have studied the toxin-antitoxin mechanism for plasmid maintenance through post-segregational killing using a theoretical model. The model extends previously studied models for sRNA-dependent regulation that have considered the limiting cases of rapid equilibrium or irreversible binding, but have thereby largely ignored the nonlinear nature of the RNA complex binding and the dynamics of the complex itself. Here we have taken the full dynamics into account and given an analytical solution for the steady state concentrations (where the plasmid copy number remains constant). We simulated the dynamics before and after plasmid loss by numerical integration and showed that the toxin-antitoxin system performs its function by inducing a transient peak in the toxin concentration, which must exceed a threshold for host killing. The formation of this peak depends on the release of the mRNA template of the protein toxin from the mRNA-sRNA complex. The accumulation of the mRNA molecules is only possible when sRNA is degraded at a faster rate than the mRNA.  We notice that the scenario of post-segregational killing of host cells that have lost the plasmid is not the only possible function of a toxin-antitoxin system of a plasmid. Our simulations of sequential loss of plasmid copies show a gradual increase of the toxin concentration without exceeding the threshold for killing. It is possible that an these sublethal toxin concentrations also contribute to plasmid maintenance: If sublethal toxin concentrations slow down growth of the host cell, this would provide additional time for plasmid replication, allowing to restore the plasmid number to the level before the loss of one or more copies. This effect is not considered in our present model, which does not included growth modulation by sublethal toxin concentrations or the control mechanism of plasmid replication, but will be studied in future work.
Using two kinds of parameter variations to study the parameter dependence of the system (individual parameter are varied in isolation, and all parameters are varied simultaneously in a random manner within a given range), we can draw a number of conclusions with regards to which ones of the 10 system parameters have significant effects on the system, and why when they take on certain ranges of values, the system functions more optimally than in other cases. In general we have found that, for an efficient toxin-antitoxin mechanism, the synthesis rate of toxin's mRNA template should be lower than that of the sRNA antitoxin, the mRNA template should be more stable compared to the sRNA antitoxin, and the mRNA-sRNA complex should be more stable than that of the sRNA antitoxin. Analytically approximating the protein peak by allowing all mRNA to be released at once gives us an analytical expression for the peak, which despite overpredicting its height nonetheless gives qualitative insights that are consistent with previous numerical observations. Finally, we also studied the possibility of inducing the toxic protein with a competitor mRNA, which sequesters the sRNA antitoxin. Such a mechanism has been proposed as an antibacterial strategy [15,16]. Here the effectiveness in killing the host cell depends on the ratio between the steady states of the toxin concentration before and after introducing the competitor RNA. An analytical solution can be given for this ratio by solving a cubic equation. We performed a single-parameter variation study for this scenario and found that the dominant parameter dependence here is on the synthesis rates. A sufficiently high synthesis rate results in all antitoxin being sequestered and the synthesis of protein toxin. In such a way the cross-talk inherent in sRNA regulation could be utilized towards inducible killing of the cells.