Modeling sRNA-regulated Plasmid Maintenance
MModeling sRNA-regulated Plasmid Maintenance
Chen Chris Gong ∗ Max Planck Institute of Colloids and Interfaces,Science Park Golm, 14424 Potsdam, Germany andInstitute of Physics and Astronomy, University of Potsdam,Karl-Liebknecht-Straße 32, 14476 Potsdam, Germany
Stefan Klumpp † Max Planck Institute of Colloids and Interfaces,Science Park Golm, 14424 Potsdam, Germany andInstitute for Nonlinear Dynamics, Georg August University of G¨ottingen,Friedrich-Hund-Platz 1, 37077 G¨ottingen, Germany (Dated: March 17, 2018)
Abstract
We study a theoretical model for the toxin-antitoxin (hok/sok) mechanism for plasmid main-tenance 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 regulationof expression of a protein toxin by an sRNA antitoxin. Here, we focus on the nonlinear nature of theregulatory circuit dynamics of the toxin-antitoxin mechanism. The mechanism relies on a transientincrease in protein concentration rather than on the steady state of the genetic circuit. Through asystematic analysis of the parameter dependence of this transient increase, we confirm some knowndesign 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, themRNA template should be more stable than the sRNA antitoxin, and the mRNA-sRNA complexshould be more stable than the sRNA antitoxin. Moreover, a short half-life of the protein toxin isalso beneficial to the function of the toxin-antitoxin system. In addition, we study a therapeuticscenario in which a competitor mRNA is introduced to sequester the sRNA antitoxin, causing thetoxic protein to be expressed. ∗ [email protected] † [email protected] a r X i v : . [ phy s i c s . b i o - ph ] S e p . Introduction Small regulatory RNA (sRNA) plays an important role in gene regulation in organismsfrom bacteria to mammals by controlling for example translation and/or mRNA stabilityand the list of known RNA regulation systems keeps increasing at a rapid pace [1–3]. sRNAregulation possesses characteristics that are distinct from protein regulation, in particular athreshold-linear response, which provides an ultrasensitive mechanism for regulatory switch-ing, 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 andplasmid maintenance in bacteria. In the first case, an sRNA controls whether synthesis of areplication primer proceeds to plasmid replication [5–7]. In the second, the plasmid encodesa (type I) toxin-antitoxin system such as the hok/sok system (“host-killing/supression-of-killing”), encoding a protein toxin and an antisense RNA which acts as an antitoxin bybinding to the toxin mRNA and blocking ribosome access, thus preventing toxin synthesis[8]. (Other types of toxin-antitoxin systems have different functions [9], in particular relatedto the formation of persister cells [10].)The toxin-antitoxin system, which is also known as an “addiction module”, maintainsthe plasmid number through post-segregational killing of plasmid-free progeny due to dif-ferential stability of the toxin and antitoxin RNAs. The killing is done by a potent proteintoxin that irreversibly damages the cell membrane [11]. In a steady state with a stable plas-mid concentration, sRNA antitoxin exists in excessive molar amount compared to targetmRNA, such that the latter is entirely sequestered in translationally inactive sRNA-mRNAcomplexes [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 avery short half-life, is rapidly depleted and the more stable target mRNA can be translatedinto toxic protein, killing the cell (Figure 1).To understand the design of the gene circuit encoding the toxin-antitoxin mechanism ofpost-segregational killing, here we analyze a theoretical model for the dynamics of sucha 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 dependenceof the genetic circuit to identify essential features and criteria for its efficient function, forexample, whether there are any criteria beyond the difference between toxin and antitoxin2NA lifetime. It also allows us to test whether the system can be designed in such away that only complete loss of the plasmids triggers killing and not a low but non-zeroplasmid 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 weperformed extensive parameter sweeps, scanning all parameters of the model individuallyas well as extensively testing randomly chosen parameter sets, thus varying all parameterssimultaneously. Overall, the model confirms the known principles for the function of thetoxin-antitoxin mechanism as the dominant ones, in particular, the difference in RNA lifetimeand the threshold condition on the synthesis rates. In addition, it also indicates a fewnew design features. Specifically, our analysis of the parameter dependence shows that thestability of the mRNA-sRNA complex plays an important role as well and should exceedthe stability of the free antitoxin. Moreover, an unstable toxic protein is superior to a stableone for induction of killing upon plasmid loss.
FIG. 1: Simplified drawing of the sRNA regulated toxin-antitoxin mechanism for plasmid maintenance.Structural transitions and processing of the mRNA have been omitted [8]. I. Model and Methods1. Dynamic Equations and Analytical Solution of Steady State Concentrations
The dynamics of a toxin-antitoxin system is modelled with four coupled ordinary differ-ential 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 whichthe mRNA is silenced. All concentrations are expressed in units of number of molecules pervolume 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-sRNAcomplex: ˙ m = α m · g − β m · m − h + · m · s + h − · c (II.1a)˙ s = α s · g − β s · s − h + · m · s + h − · c (II.1b)˙ c = h + · m · s − h − · c − β c · c (II.1c)˙ p = α p · m − β p · p (II.1d)The 10 parameters of these equations are as follows: α m , α s , α p are the synthesis ratesof 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 theplasmid copy number per cell volume.The steady state solution is obtained by setting the time derivatives on the left hand sideof the equations to zero. The steady state for this system can be explicitly given, because thenonlinear terms cancel each other in equations (1a) and (1b), leaving a quadratic equationwhich can be solved analytically. This leads to the steady state sRNA concentration: s ∗ = − A ± (cid:115) A + 4 β s β m α s g h + β c h − + β c β s h + β c h − + β c , where A = ( α m g − α s g ) h + β c h − + β c + β s β m . m ∗ = α m g − α s g + β s s ∗ β m c ∗ = m ∗ s ∗ h + h − + β c p ∗ = α p m ∗ β p (II.2)This result includes some limiting cases that have been studied in earlier work on sRNA-dependent 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 complexconcentration 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 mRNAconcentrations need to be considered, was previously studied by Levine et al. [4] and Mitaraiet al. [14].
2. Numerical Methods
To study the dynamic behavior, equations II.1(a-d) are numerically integrated using the4th order Runge-Kutta method with an integration time step on the order of 1 × − min, for a time span of 300 min. We note that the time unit of the dynamics could bemade 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 scalesor rates. To achieve both stability and speed, an adaptive time step for the numericalintegration 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 ofnegative concentration values of one or more components. The analytical results for thesteady state concentrations are in good agreement (up to floating point error) with theresults obtained from the numerical integration.5 II. Results and Discussion1. Qualitative Description of the Dynamics of the Genetic Circuit
For the model described above, the following scenario is considered: the dynamics ofthe toxin-antitoxin system begin with a cell which has just acquired one or more copies ofthe 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 asteady state of the system is reached. At a certain time point (here at t = 150 min of atotal t = 300 min), the plasmid copy number is set to 0 to mimic plasmid loss, e.g. dueto insufficient replication or unequal partitioning of the plasmid during cell division. Toillustrate the dynamics, realistic values of the parameters are used, which are estimatedfrom the literature (Table I). Equations (II.1) are numerically integrated to produce Figure2 . FIG. 2: Simulated dynamics of the toxin-antitoxin gene circuit: The different colors show the concentrations (inmolecules per cell volume) of mRNA (orange), sRNA (blue), mRNA-sRNA complex (red) and protein (green).At t = 150 min, the plasmid is lost, described by resetting the copy number g to zero, inducing a transient peakin the protein concentration. The protein concentration peak fold-increase R is . and the peak width T p equals min. The parameters are as listed in the last row of table I and g = 6 . , h + = 20 . , h − = 1 . , β c = 0 . . In Figure 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 molecules6 m [genecopy − min − ] β m [min − ] α s [genecopy − min − ] β s [min − ] α p [min − ] β p [min − ]Range of values 0.1-15.8 0.034-0.693 0.1-15.8 0.35-1.4 1-7 0.0346References [4, 17, 18] [19, 20] [4, 17, 18, 21] [8, 19, 22, 23] [8]Default parametersused here 1.0 0.2 6.0 1.0 5.0 0.035 TABLE I: Table of parameter values based on the literature for toxin-antitoxin systems and related sRNAregulation systems. For some parameters, the values are know to depend on growth conditions. For componentswith long half life, the degradation rate is dominated by dilution through cell growth and division. stops. Nevertheless, a transient increase in the free mRNA concentration is observed, whichresults in a transient increase of the protein concentration. The transient increase in mRNAis a result of sRNA degrading much faster than mRNA, so when a mRNA-sRNA complexdissociates, there is a surplus amount of free mRNA molecules released, which in turn allowsthe synthesis of toxic protein. Thus, the main function of the toxin-antitoxin system, plasmidmaintenance via the synthesis of a toxin upon plasmid loss, which results in the removal ofplasmid-free cells from the population, is dependent on a transient dynamics rather than ona steady state. This behavior is in contrast to other sRNA-based regulation systems, whichare 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 con-centration change after the loss of plasmids, two quantities are measured for each simulatedscenario: (i) the peak fold-increase of the protein concentration, R , defined as the ratio ofthe maximum of the protein concentration after plasmid loss and the steady state proteinconcentration 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.The effectiveness of the toxin-antitoxin mechanism is partially determined by the proteinconcentration 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 ofthe 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 ofprotein concentration when the plasmid number is not being properly maintained to easily7vercome such a threshold. For example, if the cell will be killed by a 10-fold-increase in toxinconcentration (10 being the killing threshold), and if there is little chance of an accidental10-fold protein concentration increase due to other factors and noise, the dynamics whichcan produce R >
10 are considered effective in exerting a robust control on the host cell.
2. Dependence of the protein dynamics on individual model parameters
We first investigate the parameter dependence of the model by varying each parameterin 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 (Figure 3) is investigated.An increase of g is equivalent to increasing the synthesis rates of the sRNA and mRNAmolecules by the same ratio. R and T p are both positively affected when g is increased. Forthe default parameters (Figure 3 (a)) as well as for the rapid equilibrium case (Figure 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 − arehigh, the changes in R and T p at the loss of one or two plasmid copies are more substantialthan in the cases of lower binding and unbinding rates. This is to say, in the case of rapidequilibrium, the model could be effective in triggering killing even when the plasmids arenot completely lost.Another way to think about plasmid maintenance is shown in the simulation of sequentialloss of plasmid copies under default parameters (Figure 3 (c)). The change in plasmidnumber does not result in the sudden release of large amount of toxin until the last plasmidis lost. Comparing to when all plamids are lost at once (Figure 2), sequential plasmid losswill result in a less prominent transient peak when the last plasmid is lost. This observationis 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 littledependence on whether there are more or fewer copies.The dependence of R and T p on the other 9 parameters are shown in Figure 4 . In general,two conditions for the existence of a transient peak of the toxic protein concentration canbe extracted from these parameter dependence:First, the synthesis rate of sRNA ( α s ) must exceed the synthesis rate of mRNA ( α m ).8his threshold condition can be observed by comparing the varying values to the defaultvalues for mRNA and sRNA synthesis rates. In Figure 4(a), when the synthesis rate ofmRNA α m exceeds the default synthesis rate of sRNA α s at 6/min/gene copy, the peakdisappears. Conversely in Figure 4(b), when α s becomes larger than α m at a default valueof 1/min/gene copy, the protein concentration peak starts to appear ( R > p peak as well as the firststeady state protein concentration p ∗ . However, p ∗ increases faster than p peak as α m increases. FIG. 3: Dependence on plasmid copy number g : (a,b) protein concentration peak fold-increase R and width T p as a function of the plasmid copy number g for our default binding and unbinding rates h + = 20 , h − = 1 (a),and for rapid equilibrium with h + = 1000 , h − = 10 (b). (c) the simulation of sequential plasmid loss underdefault parameters. The plasmid copy number is reduced from → → → → . A substantial increase intoxin concentration only occurs when the last plasmid is lost. R diminishes as α m increases.A high level of free mRNA before plasmid loss will also mean that plasmid-containing cellmight “poison itself” without any loss of plasmids. When α m > α s , there will be more freemRNA than what can be “neutralized” by complex formation, and the toxic protein will beexpressed at a high level before the loss of plasmids. The toxin-antitoxin mechanism relieson the low amount of surplus of free mRNA after the loss of plasmids, so that when α m ishigher 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 (Figures 4(c), 4 (d)), a pronounced protein concentration peak is not observed. We notice that thiscondition is not as strict as the threshold condition on the synthesis rates. FIG. 4: Dependence on individual model parameters: Each model parameter is varied in isolation while theother 9 parameters are held fixed. The toxic protein concentration peak fold-increase R (purple) and the widthof the peak T p (blue) are plotted as functions of each parameter in the 9 subfigures. The default parametervalues (from Table I and as in Figure 2) are marked by red lines in the plots. β p is very low (equivalent to having a half life of longer than 5 min), R decreases drasticallyand eventually drops to a constant around 4, where no protein is being degraded. Thus, theproteolysis of the toxic protein contributes to the efficiency of the toxin-antitoxin mechanismfor plasmid maintenance. The translation (protein synthesis) rate, on the other hand, hasno effect on the ratio R as expected (Figure 4(h)).A few surprises arise in the parameter dependence studys. First of all, it is known thatthe antitoxin sRNA in plasmid number maintenance is very short-lived [26] (which is indeedtrue 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 sRNAexceeds that of the template mRNA, in some range, the mechanism becomes actually lesseffective: R decreases in Figure 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 thefollowing section. It is therefore possible this effect only applies to situations where someother parameters, for example the synthesis rates of the RNA molecules, take on certainvalues.Another surprise comes from the complex stability. In some well-studied cases, the bind-ing of mRNA with sRNA leads to the rapid degradation of the complex [29]. However, thiscannot be the case here. As shown in 4(g), when the RNA complex is degraded at a ratehigher than 0 . R = 0). This can be explained bythe fact that mRNA must be released from the complex upon plasmid loss, which requiresa sufficiently high concentration of the complex. If the complex degradation rate β c is toohigh, 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 shownby the second parameter dependence study conducted in the following section III 3 Figure5(b). 11 . Random Sampling of the Parameter Space We have seen that the dependence of the toxin-antitoxin circuit on the individual param-eters exhibit expected but also surprising behaviors. However, due to the fact that in oursingle 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-dimensionalparameter space where the toxin-antitoxin mechanism works effectively, and to explore thejoint conditions on the parameters, we randomly sample the parameter space.We run a total of 4025 simulations of the toxin-antitoxin dynamics with randomly gen-erated 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 thedifferential equations approach is appropriate. The algorithm generates random values foreach of the 9 parameters ( g = 6 remains fixed) within the range given in Table II and withthe 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 uniformlydistributed. 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 recordedand projected onto the space spanned by α m /α s and β m /β s (Figure 5(a)). The observationsobtained by random sampling are consistent with the results of the single parameter depen-dence study. The results confirm that α m < α s is indeed a strong criterion for the existenceof a high protein concentration peak. Figure 5(a) shows when α m /α s > . β m < β s is also a very importantcriteria, 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 > α m /α s < . β m /β s <
4, and plot thoseas a function of the degradation rate of the complex β c and that of the sRNA β s . Thiswill eliminate toxin-antitoxin dynamics with parameter combinations that are inefficient inproducing a transient peak due to the criteria discussed above, allowing new effects to bediscovered more easily. Figure 5(b) shows that no peaks occur for β c (cid:38) β s , in contrast to12 IG. 5: Parameter dependence found by random parameter sampling: The color gradient shows the magnitudeof the protein peak fold-increase R . In Figure (a), the results are projected onto α m /α s vs. β m /β s on a log-logscale. In (b) the simulation results that fall within the subregion α m /α s < . and β m /β s < (marked by linesin Figure (a)) are projected onto β c vs. β s . In (c,d) the parameter sets are further restricted to the subregion α m /α s > . , β m /β s > and β c /β s < / (marked by a dark line in Figure (b)) and the results are projectedonto h + vs. h − in Figure (c) , and β p vs. β m in Figure (d), respectively. The color scale values are taken to be log . R for ease of viewing. Dark lines mark the rough division between the region with more high peaks andthe region with fewer. In (c) and (d) the dark lines mark h + = h − and β p = β m respectively. . ∗ β c (cid:46) β s .A likely explanation is as follows: when β c is low, there is a higher level of total mRNAmolecules in bound form both before and after the plasmid loss. Because the DNA genecopies are removed when the plasmids are lost, the transient increase of protein synthesisis almost entirely due to the mRNA released from the complex. Therefore, the complexdegradation (as opposed to mRNA degradation) is the dominant reason for total mRNAloss. 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 theplasmid-containing cell to “poison itself”. If this loss of mRNA due to complex degradationis slower than the loss of sRNA (the loss of sRNA is equal to the gain of the free mRNA whenthe 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 transientpeak 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 degradationof sRNA, i.e. when β s is large, R is in turn negatively affected. Considering this negativeimpact, β c (cid:46) β s therefore does not always guarantee high peaks, which is shown from awide range of variations in R within this region, and there is no increase in R when β c isincreasingly small compared to β s .It is surprising that the stability of the complex plays such an important role in thetoxin-antitoxin mechanism. It is secondary only to the threshold condition on the synthesisrates and the differential stability of the RNA molecules. This highlights the crucial roleplayed by the RNA complex in the dynamics. In connection to the single parameter studyconducted before, this demonstrates that the two known effects of unstable sRNA [26] and Parameter name α m β m α s β s h + h − β c α p β p Maximum Value 20 14 20 14 200 10 2 30 2Miminum Value 0.001 0.001 0.001 0.001 0.1 0.001 0.001 0.01 0.001Sampling Resolution 0.0005 0.0005 0.0005 0.0005 4 0.005 0.001 0.9 0.02
TABLE II: Parameter space being sampled and the rates of sampling. β c /β s (cid:46) / α m /α s < . β m /β s <
4. Figure 5(c) shows that for a prominent protein concentrationpeak to occur, h − /h + should not be larger than one. Projections onto β p and β m showthat when β p is close to 0 there is no major peak (Figure 5(d)). This means that bymaking extremely stable toxic proteins we cannot increase the relative increase in proteinconcentration. 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 stablethan proteins.The fact that high peaks do not occur in the subregion with h − /h + > α p does not have an obvious effect on the protein concentration peak. This showsthat 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 emphasizesthe importance of understanding sRNA regulation purely from the point of view of RNAdynamics instead of protein dynamics, which is traditionally viewed as playing the dominantrole in cellular regulation.To summarize, by randomly sampling the parameter space as well as performing a singleparameter dependence study, we determine that the conditions for an efficient toxin-antitoxinmechanism 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 degradationrate of sRNA;3. The stability of the complex should be higher than the stability of the sRNA antitoxin;15. A high affinity and irreversibility (but not complete irreversibility) in RNA complexformation and some degrees of protein instability are also helping factors, but are ofless significance.5. Protein synthesis rate does not contribute to the forming of a transient protein con-centration peak.Despite the nonlinear nature of our system, by individually controlling the available pa-rameters, a genetic circuit could be engineered to produce specific effects, such as a higherincrease in toxin concentration after the loss of plasmids for an effective duration (from afew minutes up to an hour).
4. Analytical approximation to the transient protein concentration peak
Using a simplified version of the differential equations after plasmid loss, an analyti-cal expression for the transient peak can be obtained, which qualitatively describes thedynamics. For this approximation, we assume that after the plasmid copies are lost, allmRNA molecules are immediately available in free form. That is, the number of free mRNAmolecules after plasmid loss equals the sum of the mRNA steady state concentration m ∗ and the complex steady state concentration c ∗ . The differential equations which describethe dynamics are as follows: ˙ m = β m · m ˙ p = α p · m − β p · p. (III.1)The initial conditions and integration constants required to solve the differential equations(III.1) 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 timecan then be obtained: m = m ∗ + ( m ∗ + c ∗ − m ∗ ) e − β m t (III.2a) p = p ∗ − ( m ∗ + c ∗ − m ∗ ) α p β m − β p e − β m t + ( p ∗ − p ∗ + ( m ∗ + c ∗ − m ∗ ) α p β m − β p ) e − β p t , (III.2b)where m ∗ , c ∗ and p ∗ are the steady state concentrations for mRNA, complex and proteinbefore plasmid loss, and m ∗ and p ∗ are the steady state concentrations after plasmid loss(see Section II 1). In the case that all plasmid copies are lost, the steady state concentrations16 ∗ and p ∗ are both zero and the above expressions can be simplified to: m = m ∗ e − β m t (III.3a) p = − ( m ∗ + c ∗ ) α p β m − β p e − β m t + ( p ∗ + ( m ∗ + c ∗ ) α p β m − β p ) e − β p t . (III.3b)These expressions qualitatively describe the transient peak of the protein concentrationafter plasmid loss, however the height of the peak is strongly overestimated (Fig 6). In thefull model, mRNA is released slowly from the complex, which reduces the peak height. Thisbehavior can be mimicked within the approximation by the sudden release of only a fractionof 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 under-stand some properties of the dynamics. For example, one can easily see from the solution(Eq.III.3b), that the translation rate does not affect the peak fold-increase R , because bothterms 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 . Thisresult should also be true for the full dynamics since the analytical approximation used hereis only considering the RNA dynamics and assumptions about the protein dynamics havenot been made. The dependence on β m − β p is also important: An increase of this quan-tity is equivalent to an increase of p ∗ , consistent with the tendency observed in the randomparameter sampling (Figure 5d).
5. Introduction of a Competitor mRNA to Increase Toxin Levels
Regulation by small RNAs typically displays crosstalk with multiple mRNAs under thecontrol of the same sRNA [30]. In the case of toxin-antitoxin systems, this can be exploitedas an antibacterial strategy [15, 16]. By inducing a gene encoding a competitor mRNA (onthe same or another plasmid or on the chromosome) that can bind to the antitoxin sNRAand hence allows toxin mRNA to exist in free form, a similar increase in toxin concentrationcan be induced as with plasmid loss, which may also lead to the killing of the cell. Todescribe the dynamics of this scenario, the equations from above are extended to include asecond type of mRNA, the competitor, and the corresponding mRNA-sRNA complex. The17
IG. 6: Analytical approximation of the dynamics after plasmid loss: With the same parameter values as inFigure 2, Solution III.3b gives a very large peak in approximated protein concentration (purple) compared to theactual simulated protein concentration change (green). It should be noted that this approximation considers noreaction time delay after plasmid loss. full dynamics is then described by the following equations:˙ m = α m · g − β m · m − h + · m · s + h − · c ˙ s = α s · g − β s · s − h + · m · s + h − · c − k + · m · s + k − · c ˙ c = h + · m · s − h − · c − β c · c ˙ p = α p · m − β p · p ˙ m = α · g − β · m − k + · m · s + k − · c ˙ c = k + · m · s − k − · c − β c · c , (III.4)where m and c denote the concentrations of the competitor mRNA and the complex itforms with an sRNA molecule, respectively. k + , k − and β c are the binding rate, unbindingrate and the degradation rate of this new complex.We assume that the competitor is induced at the time t = 150 min. Then, from t = 0to 150 min, the dynamics is the same as Equation (II.1). At t = 150 min, Equation (III.4)are being integrated, which corresponds to the competitor gene being turned on. With zero18nitial concentrations for all components (at t = 150 min, new variables r = v = 0), anexample simulation is shown in Figure 7 . FIG. 7: Introduction of a competitor mRNA: Concentrations of sRNA, mRNA, mRNA-sRNA complex,competitor mRNA, complex formed by competitor mRNA and sRNA and protein as a function of time. Notethe non-zero steady state concentrations after the triggering event resulting from the new dynamics. Theparameters in equation (III.4) are as follows: the first 10 parameters are the same as in Figure 4 . Parametervalues for competitor mRNA are chosen from within the experimental values of mRNA synthesis rate and lifetime, with synthesis rate α = 4 . and degradation rate β = 0 . . The competitor mRNA-sRNA complexbinding and unbinding rates are k + = 60 and k − = 1 . . The competitor RNA complex degradation rate is β c = 0 . . Toxin mRNA level is found to be low throughout the simulation. A transient protein concentration peak can be observed with some specific combinationsof parameters, but the majority of the simulations we run show no transient peak, as inFigure 7 . However, a non-zero steady state toxin concentration exists after induction of thecompetitor which can be effective in killing the host cell, since the mRNA encoding the toxicprotein 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 dynamicalfashion, but rather in the steady state, with the competitor either induced or not induced.Therefore, we define a new fold-increase parameter (cid:101) R of the toxin concentration as the ratiobetween the steady states of the protein concentration before and after the synthesis of the19ompetitor RNA.Similar to the procedure in Section II 1, the new steady state concentration for the mRNA (cid:102) m ∗ turns out to be a positive root to the cubic equation (III.5): α g (1 + k − β c ) m − β s β h + (cid:26)(cid:20) β (1 + k − β c ) − k + h + β m (1 + h − β c ) (cid:21) m + k + h + α m g (1 + h − β c ) (cid:27) × (cid:26) − β m β s h + m − β m (1 + h − β c ) m − ( α s − α m − α ) gh + β s m + α m g (1 + h − β c ) (cid:27) = 0 (III.5)After algebraic manipulations, it can be written in canonical form as follows. Am + Bm + Cm + D = 0 (III.6)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 explic-itly, but the expressions are excessively lengthy and hence not shown here. With a rootfinding algorithm such as the bisection method one can easily find the three roots to thecubic equation numerically. Using values similar to the experimental values, we found thateven in cases where there is more than one positive root, we can still pick out the correctsolution because usually one of the two positive roots corresponds to an unrealistically largeconcentration.The analytical solution to Equation (III.5) is then compared to the numerically generatedsteady state concentrations after integration, and very good agreement (up to floating pointerror) is found. Once we obtain the steady state mRNA concentration (cid:102) m ∗ , then following (cid:101) p ∗ = (cid:102) m ∗ · α p /β p one can solve for the steady state protein concentration. The ratio (cid:101) p ∗ / p ∗ ,where p ∗ is given by Equation (II.2), gives an analytical value for the fold-increase (cid:101) R . Thisis how we could theoretically obtain an analytical formula for (cid:101) R in the case of competitiveRNA binding, which is just as effective in killing the host cell at some triggering event asthe original toxin-antitoxin mechanism.Figure 8 shows the single parameter variation response plots generated with the samemethod as used in Figure 4 . We observe that as soon as the synthesis of the competitoris turned on, i.e. α >
0, there is a substantial fold-increase in the toxin concentration, i.e. (cid:101)
R >
1. As α becomes large, the increase in (cid:101) R is diminished because all sRNA molecules arebound. Changing the degradation rate of the competitor mRNA β results in small variationin the toxin fold-increase response (cid:101) R , meaning that for the default parameter combination,20here are very few free competitor mRNA molecules and most are bound in complexes. Thedegradation rate of the competitor mRNA-sRNA complex β c positively affects (cid:101) R , becauseit increases the destruction of the antitoxin in bound state, allowing more target mRNA tobe in free form available for protein synthesis.The parameter variation in Figure 8 shows that the synthesis rate of the new competitormRNA typically plays the dominant role compared to other new parameters, a generalfeature of cross-talk in sRNA regulation [30]. FIG. 8: Dependence on the parameters of the competitor: each new parameter of the system is increased alongthe x axis. Toxic protein concentration peak fold-increase (cid:101) R is plotted on the y axis for each run correspondingto the value the new parameter takes on. The default or fixed values of the parameters are the same as inFigure 7. IV. 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 modelsfor sRNA-dependent regulation that have considered the limiting cases of rapid equilibriumor irreversible binding, but have thereby largely ignored the nonlinear nature of the RNAcomplex binding and the dynamics of the complex itself. Here we have taken the fulldynamics into account and given an analytical solution for the steady state concentrations(where the plasmid copy number remains constant). We simulated the dynamics beforeand after plasmid loss by numerical integration and showed that the toxin-antitoxin system21erforms its function by inducing a transient peak in the toxin concentration, which mustexceed a threshold for host killing. The formation of this peak depends on the release of themRNA template of the protein toxin from the mRNA-sRNA complex. The accumulationof the mRNA molecules is only possible when sRNA is degraded at a faster rate than themRNA.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 simultaneouslyin a random manner within a given range), we can draw a number of conclusions with re-gards to which ones of the 10 system parameters have significant effects on the system, andwhy when they take on certain ranges of values, the system functions more optimally thanin 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 anti-toxin, the mRNA template should be more stable compared to the sRNA antitoxin, and themRNA-sRNA complex should be more stable than that of the sRNA antitoxin. Analyticallyapproximating the protein peak by allowing all mRNA to be released at once gives us ananalytical expression for the peak, which despite overpredicting its height nonetheless givesqualitative insights that are consistent with previous numerical observations.Finally, we also studied the possibility of inducing the toxic protein with a competitormRNA, which sequesters the sRNA antitoxin. Such a mechanism has been proposed asan antibacterial strategy [15, 16]. Here the effectiveness in killing the host cell depends onthe ratio between the steady states of the toxin concentration before and after introducingthe competitor RNA. An analytical solution can be given for this ratio by solving a cubicequation. We performed a single-parameter variation study for this scenario and foundthat the dominant parameter dependence here is on the synthesis rates. A sufficiently highsynthesis 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 induciblekilling of the cells. 22 . Bibliography [1] S. Altuvia. Identification of bacterial small non-coding RNAs: experimental approaches.
CurrOpin Microbiol , 10(3):257–61, 2007.[2] J. Livny and M. K. Waldor. Identification of small RNAs in diverse bacterial species.
CurrOpin Microbiol , 10(2):96–101, 2007.[3] Lauren S. Waters and Gisela Storz. Regulatory RNAs in bacteria.
Cell , 136(4):615–28, 2009.[4] E. Levine, Z. Zhang, T. Kuhlman, and T. Hwa. Quantitative characteristics of gene regulationby small RNA.
PLoS Biol , 5(9):e229, 2007.[5] David K. Summers.
The Biology of Plasmids.
Blackwell Science Ltd, Oxford, 1996.[6] Johan Paulsson and M˚ans Ehrenberg. Noise in a minimal regulatory network: plasmid copynumber control.
Q Rev Biophys , 34:1–59, 2 2001.[7] Yutaka Eguchi, Tateo Itoh, and Junichi Tomizawa. Antisense RNA.
Annu Rev Biochem ,60(1):631–52, 1991.[8] K. Gerdes and E. G. Wagner. RNA antitoxins.
Curr Opin Microbiol , 10(2):117–24, 2007.[9] R. D. Magnuson. Hypothetical functions of toxin-antitoxin systems.
J Bacteriol , 189(17):6089–92, 2007.[10] K. Lewis. Persister cells, dormancy and infectious disease.
Nat Rev Micro , 5(1):48–56, 2007.[11] K. Gerdes, F. W. Bech, S. T. Jorgensen, A. Løbner-Olesen, P. B. Rasmussen, T. Atlung,L. Boe, O. Karlstr¨om, S. Molin, and K. von Meyenburg. Mechanism of postsegregationalkilling by the hok gene product of the parB system of plasmid R1 and homology with the relFgene product of the E. coli relB operon.
EMBO J , 5(8):2023–29, 1986.[12] K. Gerdes, A. Nielsen, P. Thorsted, and E. G. Wagner. Mechanism of killer gene activation.Antisense RNA-dependent RNase III cleavage ensures rapid turn-over of the stable hok, srnBand pndA effector messenger RNAs.
J Mol Biol , 226(3):637–49, 1992.[13] S. Legewie, D. Dienst, A. Wilde, H. Herzel, and I. Axmann. Small RNAs establish delays andtemporal thresholds in gene expression.
Biophys J , 95(7):3232–8, 2008.[14] N. Mitarai, J. A. M. Benjamin, S. Krishna, S. Semsey, Z. Csiszovszki, E. Mass´e, and K. Snep-pen. Dynamic features of gene expression control by small regulatory RNAs.
Proc Natl Acad ci USA , 106:10655–59, 2009.[15] Julia J. Williams and Paul J. Hergenrother. Artificial activation of toxin-antitoxin systems asan antibacterial strategy. Trends Microbiol , 20(6):291–98, 2012.[16] O. R. Faridani, A. Nikravesh, D. P. Pandey, K. Gerdes, and L. Good. Competitive inhibitionof natural antisense Sok-RNA interactions activates Hok-mediated cell killing in Escherichiacoli.
Nucleic Acids Res , 34(20):5915–22, 2006.[17] S. T. Liang, M. Bipatnath, Y. C. Xu, S. L. Chen, P. Dennis, M. Ehrenberg, and H. Bremer.Activities of constitutive promoters in Escherichia coli.
J Mol Biol , 292(1):19–37, 1999.[18] S. Lin-Chao and H. Bremer. Activities of the RNAI and RNAII promoters of plasmid pBR322.
J Bacteriol , 169(3):1217–22, 1987.[19] K. Gerdes, T. Thisted, and J. Martinussen. Mechanism of post-segregational killing by thehoklsok system of plasmid R1: sok antisense RNA regulates formation of a hok mRNA speciescorrelated with killing of plasmid-free cells.
Mol Microbiol , 4(11):1807–18, 1990.[20] J. Bernstein, A. Khodursky, P. Lin, S. Lin-Chao, and S. Cohen. Global analysis of mRNAdecay and abundance in Escherichia coli at single-gene resolution using two-color fluorescentDNA microarrays.
Proc Natl Acad Sci USA , 99(15):9697–702, 2002.[21] E. G. Wagner, S. Altuvia, and P. Romby. Antisense RNAs in bacteria and their geneticelements.
Adv Genet , 46:361–98, 2002.[22] J. Jasiecki and G. Wegrzyn. Growth-rate dependent RNA polyadenylation in Escherichia coli.
EMBO Rep , 4(2):172–7, 2003.[23] S. Lin-Chao and S. N. Cohen. The rate of processing and degradation of antisense RNAIregulates the replication of ColE1-type plasmids in vivo.
Cell , 65:1233–42, 1991.[24] T. Hjalt and E. G. Wagner. The effect of loop size in antisense and target RNAs on theefficiency of antisense RNA control.
Nucleic Acids Res , 20(24):6723–32, 1992.[25] S. Nordgren, J. G. Slagter-J¨ager, and G. H. Wagner. Real time kinetic studies of the interactionbetween folded antisense and target RNAs using surface plasmon resonance.
J Mol Biol ,310(5):1125–34, 2001.[26] K Gerdes, K Helin, OW Christensen, and A Løbner-Olesen. Translational control and differ-ential RNA decay are key elements regulating postsegregational expression of the killer proteinencoded by the parB locus of plasmid R1.
J Mol Biol , 203(1):119—29, September 1988.[27] M. Brenner and T. Tomizawa. Quantitation of ColE1-encoded replication elements.
Proc Natl cad Sci USA , 88(2):405–409, 1991.[28] P. Stougaard, S. Molin, and K. Nordstrom. RNAs involved in copy-number control andincompatibility of plasmid R1. Proc Natl Acad Sci USA , 78(10):6008–6012, 1981.[29] D. H. Lenz, K. C. Mok, B. N. Lilley, R. V. Kulkarni, N. S. Wingreen, and B. L. Bassler. Thesmall RNA chaperone Hfq and multiple small RNAs control quorum sensing in Vibrio harveyiand Vibrio cholerae.
Cell , 118(1):69–82, 2004.[30] E. Levine and T. Hwa. Small RNAs establish gene expression thresholds.
Curr Opin Microbiol ,11(6):574–9, 2008.,11(6):574–9, 2008.