Excitable dynamics through toxin-induced mRNA cleavage in bacteria

Toxin-antitoxin (TA) systems in bacteria and archaea are small genetic elements consisting of the genes coding for an intracellular toxin and an antitoxin that can neutralize this toxin. In various cases, the toxins cleave the mRNA. In this theoretical work we use deterministic and stochastic modeling to explain how toxin-induced cleavage of mRNA in TA systems can lead to excitability, allowing large transient spikes in toxin levels to be triggered. By using a simplified network where secondary complex formation and transcriptional regulation are not included, we show that a two-dimensional, deterministic model captures the origin of such toxin excitations. Moreover, it allows to increase our understanding by examining the dynamics in the phase plane. By systematically comparing the deterministic results with Gillespie simulations we demonstrate that even though the real TA system is intrinsically stochastic, toxin excitations can be accurately described deterministically. A bifurcation analysis of the system shows that the excitable behavior is due to a nearby Hopf bifurcation in the parameter space, where the system becomes oscillatory. The influence of stress is modeled by varying the degradation rate of the antitoxin and the translation rate of the toxin. We find that stress increases the frequency of toxin excitations. The inclusion of secondary complex formation and transcriptional regulation does not fundamentally change the mechanism of toxin excitations. Finally, we show that including growth rate suppression and translational inhibition can lead to longer excitations, and even cause excitations in cases when the system would otherwise be non-excitable. To conclude, the deterministic model used in this work provides a simple and intuitive explanation of toxin excitations in TA systems.


Introduction
Toxin-antitoxin modules are small genetic elements, omnipresent on the genomes of bacteria and archaea, that code for a small intracellular toxin and its counteracting antitoxin [1][2][3]. The antitoxin typically has a higher in vivo turnover rate than the toxin [4]. In type II toxin- antitoxin modules, both the toxin and the antitoxin are proteins and the toxin neutralization occurs through the formation of non-toxic complexes [5]. In several toxin-antitoxin modules one antitoxin can neutralize up to two toxins, forming either the complex AT or the complex TAT. Toxin-antitoxin modules further have an intricate transcriptional regulation: the antitoxin has a DNA-binding domain with which it can bind to the promoter/operator region of the toxin-antitoxin module, and functions as a weak repressor. The toxin can function as a corepressor or a derepressor for the antitoxin, depending on the toxin:antitoxin ratio [2]. Different toxins have different targets in the cell, for example, CcdB poisons DNA gyrase [6], while MazF and RelE cleave mRNA [7][8][9]. Such endoribonuclease toxins will be the focus of this paper. Although toxin-antitoxin modules are widespread in prokaryotes, their biological role is currently still unclear. Toxin-antitoxin modules have been implicated in plasmid maintenance, abortive phage infections, the response of bacterial cells to nutritional stress and the formation of persister cells [3,10]. These are cells that are tolerant to multiple antibiotics because they are in a temporary state of dormancy [11]. Although previously all known type II mRNA endoribonuclease toxins in E. coli K-12 were proposed to be involved in persistence, the role of these toxin-antitoxin modules in persister generation in the absence of stress is currently uncertain [3,12]. Computational studies can be useful to gain insight into the possible dynamics caused by the architecture of the genetic network and the protein-protein, protein-DNA and protein-RNA interactions in a toxin-antitoxin module. Several groups have studied toxin-antitoxin modules computationally, using either deterministic [13,14] or stochastic [15,16] approaches. From these modeling efforts, two possible deterministic explanations have emerged for the elevated free toxin levels that might be linked to the generation of persisters. First, it is plausible that there is bistability between a growing, antitoxin-dominated state and a toxin-dominated state [13,14,16,17]. A critical component to allow the existence of a toxin-dominated state is that higher free toxin levels decrease the cellular growth rates, which in its turn has an effect on the accumulation rate of the toxin itself. Increased noise levels in the presence of stress could lead to stochastic switching between these two states. A second possibility is that the toxindominated state only exists as a transient excursion in the free toxin level [15]. Such deterministic excursions could theoretically be generated through a process called excitability, where noise could act to trigger them. Furthermore, if toxins induce growth rate reduction, the duration of such toxin excursions could be significantly lengthened. So far, theoretical studies have only observed such transient toxin excitations using stochastic simulations [15]. As for low molecule numbers the deterministic limit of stochastic models does not always give an accurate description of the real dynamics [18], a potential link to deterministic excitability remains to be shown. Finally, it is important to note that these different types of deterministic dynamics aim at describing the behavior of single cells. Both bistability and excitability can give rise to bimodal distributions on a population level.
In this article we focus on the effect of the cleavage of mRNA in the presence of elevated free toxin levels, which has recently been shown to cause toxin excitations [19]. We use a simplified system, where we leave out the formation of the complex TAT and the transcriptional regulation, as this is the simplest toxin-antitoxin model system that still displays the spikes in the free toxin concentration. We use a deterministic set of differential equations to describe this system and show that it yields similar results as simulations with a Gillespie algorithm. Combining a deterministic approach with a separation of time scales in the system allows to further reduce the problem to a two-dimensional system, which can be visualized and interpreted in the phase plane. Moreover, it facilitates bifurcation analyses that show how changes in system parameters affect the TA dynamics. For example, we verify how nutritional stress, which causes an increase in the degradation rate of antitoxin, influences free toxin spikes.
Finally, we examine how additional feedback mechanisms like transcriptional regulation by binding of the toxin-antitoxin complex to DNA, or growth rate modulation and inhibition of translation at high free toxin levels can affect the dynamics of the system.

Materials and methods
Ordinary differential equations are simulated with the integrate function of the package scipy [20] of python (Python Software Foundation. Python Language Reference, version 2.7. Available at http://www.python.org). In order to perform bifurcation analysis we use a Newton-Raphson method to solve a set of equations using the jacobian of the system [21].
The stochastic simulations were performed using a Gillespie algorithm, based on treating the biochemical reactions as discrete stochastic events [22], implemented in MATLAB. The concentrations are converted to number of molecules using a volume factor of 3:612 � 10 8 molecules=cell M . This factor is based on an estimated volume of an E. coli cell of 0.6 (μm) 3 [23] and calculated as 6:02 � 10 23 molecules mol � 1000 l m 3 � 6 � 10 À 19 m 3 cell . A detailed description of the different deterministic and Gillespie models used can be found in the supplementary material (S1 File).

Toxin-induced mRNA cleavage leads to toxin excitations
Recently, it has been demonstrated experimentally that the mRNA cleavage by MazF toxins can lead to cell growth heterogeneity [19]. Using stochastic Gillespie modeling, the authors have shown that the growth rate inhibition during stressful conditions may be linked to excitations in the free toxin level, which are only found in simulations if the cleavage of the mazEF mRNA itself is included in the model. In this paper, we expand the analysis of this model [19]. Briefly, it describes the transcription of the operon and the translation of the formed mRNA to produce the antitoxin MazE (A) and the toxin MazF (T), taking into account that both proteins form dimers in solution and that the antitoxin is produced at a higher rate than the toxin [4]. The model also includes the dilution of all components due to cell division. Only the mRNA and the antitoxin have a separate, higher degradation rate since these species are less stable than the toxin. The toxin and the antitoxin can bind to form the non-toxic complexes AT or TAT, the latter corresponding to the MazF 2 -MazE 2 -MazF 2 complex described by Kamada et al. [24]. Since this complex formation has a stabilizing effect on the antitoxin, the degradation rate of the antitoxin in complex is lower than that of the antitoxin alone. Additionally, there is a negative transcriptional regulation, with complex AT binding to the DNA and inhibiting the transcription of the operon. Finally, once the free toxin level exceeds a certain threshold, the toxin will cleave the mRNA of the mazEF operon. We made only three adaptations to the model used to describe the mazEF toxin-antitoxin system in [19], for the sake of simplicity: when DNA binding is included, only one binding site is considered instead of three; we omit the binding of the antitoxin alone to the DNA, which is known to be much weaker than the binding of toxin-antitoxin complexes [25]; and we exclude the slow unbinding of toxin-antitoxin complexes.
In Fig 1 we show the different systems that will be used in this article in order of increasing complexity. The model described above, corresponding to the model used in [19], is shown in the middle panel (Fig 1C). At the top of Fig 1A we show the reduced model, which is our main focus in this article. Finally, for a biologically more complete picture, we analyze the effect of growth rate and translation inhibition at high free toxin levels, as shown in the bottom panel ( Fig 1E).
The resulting behavior is shown in Fig 1B, 1D and 1F for the situations with and without mRNA cleavage, where the highlighted zones represent an episode with an increase of stress. As described above, for the model used in [19], we find that, for increased stress, toxin excitations are mainly present if the cleavage of mRNA by the toxins is included in the model (Fig 1D), suggesting that mRNA cleavage plays an essential role in triggering large spikes in the free toxin levels. This made us wonder whether other interactions such as transcriptional regulation by binding of AT to the DNA and the formation of the second complex TAT are also required to generate toxin excitations. When removing these interactions altogether, we found that the qualitative behavior of the system was not affected (Fig 1B). When including A: Scheme of the minimal system that we use to explain toxin excitations. DNA is transcribed to mRNA which is translated to toxins T and antitoxins A. These can combine to form the complexes AT. If the toxin level reaches a threshold this will cause cleavage of the mRNA. B: Gillespie simulation of the minimal system. In stress conditions (highlighted region) the system is only excitable when mRNA cleavage is included. C: Scheme of the system when including secondary complex (TAT) formation and DNA binding of complex AT, inhibiting the mRNA transcription. D: Gillespie simulation of the system with TAT formation and DNA binding by AT, with and without mRNA cleavage. E: Scheme of the entire system, including translational inhibition as well as inhibition of the growth rate when the toxin level is increased. F: Gillespie simulation of the entire system, with and without mRNA cleavage. Without mRNA cleavage there are increased toxin levels, but pronounced excitory behavior is observed when including mRNA cleavage. For the details of the Gillespie simulations, see supporting information S1 File. Parameters: see Table 1 mRNA cleavage, similar toxin spikes were observed, while they disappeared when also abolishing the mRNA cleavage. These simulations indicate that toxin-induced mRNA cleavage is one mechanism by which toxin excitations can be triggered. Finally, we wondered whether this behavior persists when adding additional known interactions, such as growth rate modulation and inhibition translation. The growth rate of the cell can be slowed down when toxin levels are increased. This has been included in models before, for example by Cataudella et al. [13]. They show that a combination of growth rate inhibition and conditional cooperativity can lead to a bistable system, with one healthy low-toxin state and one high-toxin state. Additionally, high levels of the MazF toxin can have an influence on the global translation rates in the cell, since MazF also cleaves the transcripts for ribosomal proteins and rRNA precursors [26]. If we include these mechanisms in the Gillespie simulations, we observe toxin excitations, whereby the excitation time is increased due to the growth rate inhibition (Fig 1F), corresponding to the result obtained in [15]. The effect of these mechanisms are analyzed in more detail in the last section.

Toxin excitations are the result of underlying deterministic excitability
In order to better understand the dynamical origin of the observed spikes in toxin levels (Fig 1), we constructed and analyzed a deterministic model of the reduced TA network, shown in Fig 1C, consisting of ordinary differential equations (ODEs). By not considering the formation of the complex TAT and the transcriptional regulation for now, we focus on analyzing the role of toxin-induced mRNA cleavage in generating toxin excitations. The reduced deterministic model is given by the following 4 ODE equations: The parameters and variables can be found in Tables 1 and 2, where most parameters have been experimentally measured or motivated, see [15,19]. The parameters are identical to those in [19] unless specified below. DNA is transcribed to mRNA [M] with a constant transcription rate r F , 0.121 s −1 , based on a transcription rate of 70 nucleotides per second [27] and the length of the transcript. This mRNA is then translated to antitoxins and toxins with rates b 1 , 0.122 s −1 , and b 2 , 0.009 s −1 . The toxin translation is lower than the antitoxin translation, this difference was measured experimentally for several type II toxin-antitoxin modules including mazEF using ribosome profiling [28]. A toxin and an antitoxin can bind to form a complex AT with a rate a T , 1.32 × 10 5 M −1 s −1 , corresponding to 0.000365 s −1 after applying the volume factor. This value is based on binding experiments using Surface Plasmon Resonance (SPR) (personal communication) and lies within the expected theoretical range [29,30]. The mRNA has a degradation rate d m of 0.002 s −1 under normal conditions (in the absence of stress conditions), based on a half-life of 5.7 minutes [31]. If the toxin level reaches a threshold, then it is assumed that mRNA cleavage is activated, which is modeled by an increase of the degradation rate (d large ) and by using a Hill-function with threshold K t and coefficient n. Here, we assume that the maximal degradation of the mRNA in the presence of the toxin is 100 times larger than d m (the basal degradation rate of mRNA), while it was fixed at 10 molecules s −1 in [19]. The Hill coefficient on the other hand is fixed at 2 for all simulations in this paper unless specified differently, while it was 3 in [19]. The threshold K t is fixed at 15 toxins per cell, as in [19]. The toxin and complex AT are considered as stable and are diluted at a rate d c of 0.00028 s −1 as the result of cell division, with an estimated cell cycle length of 40 minutes. The antitoxin is degraded with a higher rate d a , which was estimated as 8 times d c . Although the antitoxin is more stable within toxin-antitoxin complexes, we assume that it can still be degraded in this case with a lower rate d a2 . As this depends on the degradation rate of the antitoxin, we assume that the ratio d a2 d a ¼ 0:1 is fixed. Besides the difference in the degradation rates of the toxin and the antitoxin, the toxin translation rate is smaller as well, introducing a difference in the time Table 1. Redefinition of the parameters for the normalized model and the correspondence to the original parameters. These parameters are valid for non-stress conditions. In case of stress, δ a , d a and/or ε are adjusted as indicated in the figure captions.

Normalized model
Original model scales of A and T. In order to quantify this, we use ε ¼ b 2 b 1 ⪡ 1, and reformulate (1) as follows: This model is now normalized in a way that the parameters are approximately of order O (1), except for β and ε. and will approach their equilibrium condition relatively quickly in comparison to the toxin (x) and the complex AT (y) (O(ε)). As a result, we assume that dm dt ¼ 0 and da dt ¼ 0 in order to describe the long term dynamics. This is often called a quasi steady state approximation. Carrying out this further simplification the system becomes: with mðxÞ ¼ 1 Since this reduced system of Eq (3) is two-dimensional, the time evolution of the toxin x and complex y can be depicted in a simple two-dimensional graph (see Fig 2A). For every set of initial conditions (x = x 0 , y = y 0 ) in the plane, the Eq (3) predict how the system will evolve. In Fig 2, these predictions are illustrated by the gray arrows, which represent the so-called "flow" of the system. The arrows not only show the direction in which the system evolves, but also the speed with which it does so (larger arrows represent faster dynamics). Such a twodimensional flow allows us to immediately see by eye how the system will evolve in each situation. Note that the scale used in this figure is logarithmic. Another useful tool to understand the behavior of two-dimensional systems is to plot the nullclines. Nullclines consist of all points (x, y) where dx dt ¼ 0 (NC(x), black line) or dy dt ¼ 0 (NC(y), gray line), as given by: The intersection of both nullclines represents the equilibrium situation of the system, which is called a fixed point. In our case, there is only one such fixed point and it is stable in the sense that small perturbations will immediately return to this equilibrium. In fact, as it is the only fixed point, all initial conditions (all points in the plane) will eventually return to this equilibrium state. However, when perturbing the system such that the dashed part of NC(x) is crossed, the system will make a large excursion in the plane before returning to the fixed point (see red line in Fig 2A). The excursion in the plane corresponds to a spike in the toxin levels, see Fig 2B. This behavior is called excitability [32].
Strictly speaking this analysis in the plane only applies to the reduced two-dimensional system (3). However, we also evaluated the full four dimensional system (2) and projected its dynamics onto the plane (x, y) (see blue line in Fig 2A). The full model behaves very much alike, exhibiting similar excitable dynamics. This confirms the validity of the quasi steady state approximation used in deriving the reduced model system. Small differences are observed for low toxin levels, in the regime x < κ, because the time scale separation used in our reduction no longer applies. In this regime, a different normalization could be used, see supporting information, S2 File.
What happens biologically is that when crossing the threshold for mRNA cleavage (x > κ), mRNA cleavage quickly reduces the levels of mRNA (m), which in turn suppresses the translation of antitoxin (a) and toxin (x), thereby preventing the formation of complexes (y). Even though toxin translation is suppressed, at first toxin levels keep increasing as the complexes fall apart with a rate δ AT . When there is an insufficient amount of complexes left, the dilution rate δ c will dominate so that the toxin level decreases. When the level of toxins drops below the mRNA cleavage threshold again, the toxin excitation is terminated, and the system relaxes to the fixed point with low levels of toxins. These rates δ c and δ AT determine the shape of the toxin excitation (Fig 2B) as can also be seen analytically by simplifying the Eq (3) in the limit x > > κ as follows (for details see supporting information, S2 File): While our findings show that toxin excitations in a deterministic system are related to excitability, it does not yet prove that the stochastic toxin spikes we observed in Fig 1 and in [19] were related to underlying deterministic excitability. Indeed, various works have shown that excitations can also be a purely stochastic phenomenon without the need for any corresponding deterministic oscillatory behavior [33][34][35]. In order to test this, we performed simulations with a Gillespie model, using the parameter values corresponding to the ODE model in Fig 2, for different values of the antitoxin degradation rate δ a . We then visualized the data obtained from the stochastic simulation as a heatmap in the plane (x, y), thus plotting the probability to find the system at a certain location in this two-dimensional space. The results are shown in Fig 3A, where the heatmap illustrates the most probable trajectory of the toxin excitation. This most probable path in phase space nicely overlaps with our deterministic predictions. The corresponding time traces of the deterministic and the Gillespie simulations are represented in Fig 3B and 3C. Together, these simulations confirm that the toxin pulses observed in Gillespie simulations correspond to stochastically triggered pulses, existing due to underlying deterministic excitability.

The toxin-antitoxin system can be monostable, excitable or oscillatory
Next, we explored the robustness of these toxin excitations to changes in the system parameters. In Fig 3, we show that the system can display qualitatively different behavior when changing the antitoxin degradation rate (δ a ). Note that δ AT varies along with δ a , as the ratio d AT d a is fixed. For δ a = 1, the system has one stable solution and stochastic Gillespie simulations do not show any toxin excitations (Fig 3A). This is explained by the fact that the excitation threshold is too large, such that noisy excursions around the stable state are too small to trigger any excitation. We note that the projected Gillespie data does seem to occasionally cross the threshold (dashed line), even though this does not lead to an excitation. This discrepancy between the Gillespie simulations and the reduced 2D ODE system is explained by the fact that perturbed values of m and a in the Gillespie simulation do not immediately reach their quasi steady state condition. Whereas an excess of toxins immediately leads to mRNA cleavage in the reduced deterministic case, there is a small delay in the Gillespie simulations. As long as the antitoxin level is sufficiently large, complex formation is able to reduce the toxin level, thus having a stabilizing effect on the system. When increasing δ a , the excitability threshold is reduced, such that the probability to stochastically excite toxin spikes increases. Indeed, for δ a = 2, toxin excitations can be observed in the Gillespie simulations (Fig 3B). Although the deterministic 2D model accurately predicts the shape of the toxin excitations, the equilibrium state is often not fully reached in the Gillespie simulations. The reason is that the fixed point is situated close to the top of the nullclines, causing lower complex levels when a perturbation crosses NC(y), or an excitation when NC(x) is crossed. When the initial condition is such that it corresponds with a point that is situated in the main band of the Gillespie data, then the deterministic trajectory of the first excitation does correspond well with the Gillespie simulations. This discrepancy disappears when the fixed point is situated further from the top of the nullclines (see also supporting information, S2 Fig).
Finally, for even larger antitoxin degradation rates (e.g. δ a = 3), the reduced ODE system shows oscillatory behavior (Fig 3C). The fixed point becomes unstable and instead the system converges to a limit cycle. Consistent with these deterministic findings, Gillespie simulations also show more regular excitations.

Cellular stress causes more excitations
When a bacterial cell is experiencing nutritional stress such as amino acid starvation, the degradation rate of several antitoxins (δ a ) increases due to the increased activity of cellular proteases such as Lon [7,36]. As shown before in Fig 3, this increases the probability of toxin Excitable dynamics through toxin-induced mRNA cleavage in bacteria excitations. Another way to influence the toxin level is by directly increasing the translation rate of the toxin (by varying ε). Here, using time evolution simulations and bifurcation analysis, we analyze how changes in the system parameters δ a and ε affect the toxin-antitoxin dynamics. Other important biological parameters that affect the dynamics are the parameters determining the mRNA cleavage (β, n, κ) and the binding parameter α. The influence of these parameters on the dynamics and shape of the nullclines are discussed in S1 Fig. In order to analyze the results of the bifurcation analysis (Fig 4) we define the excitation time as the time that the toxin level is higher than the toxin threshold κ used in the Hill function (see Fig 2B). Fig 4A and 4B show the average excitation time of toxin excitations, obtained with deterministic (A) and stochastic Gillespie (B) modeling, as a function of the antitoxin decay rate δ a and toxin translation rate ε. Our simulations show that excitations occur for moderate stress levels, while higher stress leads to oscillatory behavior.
When stress levels are increased even further, the fixed point becomes stable again, although we are currently in a state where free toxin dominates the behavior of the system, and this situation would not be viable for a cell. The loss of oscillations is less clear in the Gillespie data, as some excitations are still detected. This is explained by the fact that the toxin value x at the fixed point lies close to threshold κ, such that noise can still trigger excursions.
Interestingly, in the Gillespie simulations the excitation time remains constant over a large range of ε and δ a , whereas in the deterministic case the excitation time is larger for lower values of δ a and ε. The reason is that here the fixed point is situated closer to the top of the nullclines, so that the amplitude of an excitation is slightly larger than in the Gillespie simulation. This is analyzed in more detail by keeping ε = 0.074 fixed in Fig 4C and by keeping δ a = 2 fixed in Fig 4D. The bifurcation analysis, using Eq (3), shows that the fixed point loses its stability in a supercritical Hopf bifurcation. We found that the amplitude of the limit cycle is increased almost instantaneously in what is called a Canard explosion [37]. The excitability of this system is related to the vicinity of a limit cycle in parameter space. As the period is nondiverging at the bifurcation point, this is an example of type II excitability [38]. The maximum toxin level during an oscillation, the excitation time and the period of an oscillation are compared with the observation in the Gillespie data. They show a good correspondence, although the period of the Gillespie data is larger, as it takes time for a perturbation to cross the excitation threshold. In Fig 4C and 4D we clearly see that the period decreases for increased stress, so that excitations are more frequent. We see a decrease in the excitation time in the deterministic model, while in the Gillespie data the excitation time is more constant. This discrepancy is due to the fact that the Gillespie data does not fully reach the fixed point, as explained in last section. This behavior is illustrated in more detail in S2 Fig, where we explore the system dynamics for different values of ε.

Toxin excitations persist when including additional complex formation and DNA binding
So far, we showed that toxin-induced mRNA cleavage is the main mechanism leading to excitable behavior. Here, we will incorporate the effect of a second complex TAT and the binding of AT to the DNA and show that this does not change the dynamics in a qualitative manner (Fig 1). Each deterministic and corresponding stochastic model is described in more detail in S1 File.
Inclusion of a secondary complex. The system (3) can be extended to incorporate the second complex TAT by using an additional variable z ¼ d m ½TAT� ε . We assume that y and z have the same creation rate α. We also assume that z is reduced with rate δ AT to two toxins x due to antitoxin degradation within toxin-antitoxin complexes. The difference in time scales still Excitable dynamics through toxin-induced mRNA cleavage in bacteria exists, so that m and a can be assumed to be in steady state, and the resulting normalized equations are the following: To show the similarity of this three-dimensional ODE model with the two-dimensional model (3), we look at the total amount of toxin in the complexes, c = y + 2z rather than the dynamics of y and z separately, resulting in the following system: The only difference between these equations and (3) are the terms of α(a+ y)x. However, during an excitation the terms of order O(ε) become important, which have the same shape as in (3). Hence, the mechanism of the excitation does not change by inclusion of the secondary complex z. This is explored graphically in Fig 5, where we show the phase plane representation and the corresponding time traces for the deterministic and stochastic models, for low (δ a = 1.5) and high (δ a = 3) stress levels. The altered term −α(a + y)x does have a stabilizing effect on the system dynamics, as x will be depleted more quickly. This causes a shift in the nullclines such that the excitation threshold is larger if the secondary complex is included in the model (Fig 5A). Therefore, the excitation frequency is lower than in the reduced model when stress is increased (Fig 5B). The plotted nullclines in Fig 5 are approximations, assuming that a and y are of the same order of magnitude. This way the terms α(a + y)x reduce to 2αax. As a result there is not an exact correspondence with the intersection of the nullclines and the fixed point of the system.
Inclusion of binding to DNA. Transcriptional regulation occurs when the complex AT (y) binds to the DNA and thereby inhibits the transcription of mRNA. As such, the translation of the toxins and antitoxins is also indirectly inhibited. In the deterministic system, we model DNA binding by incorporating an inhibiting Hill function r/(r + y) into the translation rate (for derivation, see supporting information, S1 File). This expression decreases from 1 (no suppression) to 0 (complete suppression) for increasing values of the complex y. When assuming quasi steady state conditions for m and a, the system is still described by (3), but the expression for m(x) is changed as follows: As the binding inhibits the translation of toxin and antitoxin, there is less complex formation. This causes a shift in the nullclines leading to an equilibrium with less complexes (Fig 5). Binding of AT to the DNA did not stabilize the system as the threshold for toxin excitations did not increase. Moreover, it seems that inclusion of DNA binding leads to excitations for parameter values where the reduced system is still stable (Fig 5A). Even though there is no stabillizing effect, the advantage of DNA binding is that the cell uses less energy as there is less translation [15]. If both the secondary complex TAT and DNA binding are included in the system, there are less complexes due to the binding and the system is more stable due to the presence of the secondary complex TAT (Fig 5A). When the amount of stress is increased, the excitation frequency is increased as well (Fig 5B), as is the case in the reduced model.

Inclusion of growth rate inhibition allows for excitability without a sigmoidal response of mRNA cleavage by toxins
For the simulations we kept the Hill-coefficient fixed at n = 2. A Hill-coefficient larger than one leads to the bend in the nullclines, so that an unstable branch of NC(x) exists in the phase plane (see also S1 Fig). As explained before, this determines a threshold for an excitation.
When the Hill-coefficient is equal to one, the toxin dependence of the mRNA cleavage is hyperbolic rather than sigmoidal. As a result, the shape of the nullclines is such that no excitable behavior is observed in the models used above. However, even in this case, the system can still be excitable if we include extra biologically relevant mechanisms in these models.
One such mechanism is the slowing down of the growth rate due to increased toxin levels. Cataudella et al. have shown that growth rate inhibition, in combination with conditional cooperativity potentially leads to bistable behavior [13]. We used similar functions as in [13] to incorporate these effects in our model: the translation rate is scaled with a function f m ¼ S m 1þB m x , whereas the dilution rate is scaled with a function f t ¼ S t 1þB t x (for the detailed model, see supporting information S1 File). In Fig 6 we study the effect of these inhibitory functions for Hill coefficient n = 1 ( Fig 6A) and for Hill coefficient n = 2 ( Fig 6B). In contrast to the result of Cataudella et al. [13], we do not find bistability when including these inhibiton functions, which is probably related to the fact that we use different time scales in our model than Cataudella et al., where it is assumed that complex formation happens instantaneously. The growth rate inhibition leads to increased excitation times, corresponding to the observation in Fig 1F. We find that for n = 1 translational inhibition is needed for excitable behavior. The decrease of the dilution rate does not affect the presence of excitable behavior, but causes the cell to be in the toxic regime for a longer period of time. In all the cases where excitable behavior exists, the nonlinearity of the nullclines is of at least degree two, causing the characteristic bend in the nullcline NC(x). This can be obtained by combining mRNA cleavage with n = 1 and translational inhibition, or by using mRNA cleavage with a Hill-coefficient of at least two.

Conclusions
TA modules are small dynamic systems, coding for a toxin and its corresponding antitoxin [1][2][3]. The toxin level can affect a cell in different ways: post-segregational killing, abortive phage infections and the formation of persister cells [3,10], although the latter is currently unclear [3,12]. There exist different dynamic explanations for elevated toxin levels in cells. The first possibility is that this corresponds to a second equilibrium state, due to bistable behavior [13,14,16,17], while the second possibility is that stochastic perturbations cause pulses in the free toxin level, corresponding to excitable behavior [15]. Recently it was found that mRNA cleavage by the toxins can cause toxin excitations, leading to cell growth heterogeneity [19].
In this article we used a deterministic model to analyze the role of mRNA cleavage as a possible mechanism behind such excitations. As there is a difference in time scales between toxin and antitoxin translation and degradation, the model can be simplified to a set of two ordinary differential equations (ODEs), allowing a description in the phase plane. An excitation occurs when a threshold is crossed, mRNA cleavage is switched on, and the repression of translation prevents an immediate return to the equilibrium state. By systematically comparing with Gillespie simulations we showed that even though the real system is inherently stochastic, a deterministic model is capable to describe the observed dynamics. When molecule numbers are low, a deterministic limit of a stochastic model does not always give an accurate description of the dynamics in the presence of multiple time scales [18,39,40]. Moreover, it has been shown that oscillatory behavior in noisy systems can be created by the noise, when no oscillations are observed in the deterministic model [33][34][35]. In this article, we show that the toxin excitations are caused by the underlying deterministic excitability, due to the vicinity of a Hopf bifurcation where a limit cycle is created. Stress can be modeled by varying the antitoxin degradation rate and the toxin translation rate, which increases the probability of excitations, as the fixed point gets closer to the unstable branch in the phase plane. In conclusion, even though this system is inherently stochastic, we provided a deterministic description of the excitable behavior in TA modules due the presence of toxin-induced cleavage of mRNA.
Similar excitable behavior in bacteria was theoretically and experimentally analyzed in the ComK-ComS gene regulatory circuit in Bacillus subtilis, where excitability led cells to be in a transient state in which they were competent to take up DNA from the environment [41,42]. Although the circuitry of interacting genes and proteins in the ComK-ComS system is significantly different than that of the TA systems we studied here, the type of excitable behavior is similarly caused by a combination of fast positive and slow negative feedback loops. By using quantitative fluorescence time-lapse microscopy to observe circuit components in individual cells, and comparing such measurements with mathematical models, significant new insights were gained into how the ComK-ComS gene regulatory circuit works [41,42]. Excitations in the free toxin level due to a difference in time scales also occur in the context of type I toxinantitoxin modules, and they play a role in plasmid maintenance due to post-segregational killing [43,44].
We hope that our model will trigger new experimental efforts in the field of TA systems, especially to measure the dynamics of circuit components on a single cell level, and that they will help in shedding new light on the temporal dynamics of cellular toxin and antitoxin levels and growth rates. Devising the necessary fluorescent reporters to track these protein levels in vivo is however no trivial task [45]. The reason is that adding a fluorescent tag to the toxin or the antitoxin may significantly alter the working of the operon, since both are small proteins, which have to be able to interact with each other, the DNA (for the antitoxin) and the target (for the toxin). An alternative approach is the use of a transcriptional fluorescent reporter, in which the promoter of the toxin-antitoxin system is fused to a fluorescent reporter gene, giving information about the transcription from the operon instead of the protein levels. However, many stress conditions interfere with translation, which may lower the readout from this type of fluorescent reporter. Still, we believe that such efforts would be very valuable, since they would allow to bridge the internal dynamics in individual cells and the dynamics on the level of whole cell populations, where bimodal distributions of fast and slowly dividing cells have been observed [19].
Supporting information S1 Fig. Influence of the parameter n, β, κ and α on the phase plane. These parameters can change the shapes of the nullclines: an increase in n makes the slope steaper, an increase of β leads to a bigger difference between the nullclines in the toxic and normal state, κ changes the threshold itself and α affects the complex level at steady state. However, the overall behavior as explained in the main text remains the same.