A Minimal Model for Gene Expression Dynamics of Bacterial Toxin-antitoxin Systems

Toxin-antitoxin (TA) modules are part of most bacteria’s regulatory machinery for stress responses and general aspects of their physiology. Due to the interplay of a long-lived toxin with a short-lived antitoxin, TA modules have also become systems of interest for mathematical modelling. Here we resort to previous modelling efforts and extract from these a minimal model of TA dynamics on a timescale of hours, which can be used to describe time courses derived from gene expression data of TA pairs. We show that this model provides a good quantitative description of TA dynamics for the 11 TA pairs under investigation here, while simpler models do not. Our study brings together aspects of Biophysics with its focus on mathematical modelling and Computational Systems Biology with its focus on the quantitative interpretation of ’omics’ data. This mechanistic model serves as a generic transformation of time course information into kinetic parameters. The resulting parameter vector can, in turn, be mechanistically interpreted. We expect that TA pairs with similar mechanisms are characterized by similar vectors of kinetic parameters, allowing us to hypothesize on the mode of action for TA pairs still under discussion.


Introduction
The vast majority of free-living bacteria contain a number of toxin-antitoxin (TA) gene pairs [1][2][3][4] . The toxin products target key cellular functions inhibiting cell growth and eventually leading to cell death, while the corresponding antitoxin neutralizes the toxin's effect, thus, forming a toxin-antitoxin system whose accurate expression regulation is vital to the survival of the cell 5 . These TA systems are currently classified in six groups (types I, II, III, IV, V, VI) 2 according to the mechanism used by the antitoxin to neutralize the toxin. Types I-III are considered to be well-established TA systems 3, 6-9 while types IV-VI consist of newly discovered types [10][11][12][13][14] . Type II TA systems are the largest and best studied TA system class. Type II antitoxins are proteins. They typically have two domains, one that binds DNA and a second that binds and inhibits the activity of the cognate protein toxin 2,3,9 . The presence of toxin-antitoxin systems is considered to be associated to persistence, i.e. the multidrug tolerance of bacteria, which obviously compromises the effectiveness of antibiotics on many pathogenic bacteria 15 . It is believed 4,15,16 that when antibiotics are applied, a small sub-population of bacteria, called persisters, enters a dormant, non-dividing state and thus are protected from being killed. Experiments have shown a connection between persister formation and the competition between a toxin and its antitoxin inside an E. coli cell. Toxins inhibit cell growth and most antibiotics target the cell during the growth phase. Cells entering this persistent state seem to be immune to antibiotics but this immunity is different from the one obtained through advantageous mutations that result in antibiotic resistance since it is not permanent or inherited 17 . Knowledge about toxin-antitoxin systems in bacteria is still accumulating 18 . This is true for the discovery of new toxin-antitoxin modules 19 , their classification 5,20 , their functional roles 21,22,25,26 as well as their detailed molecular mechanisms 27 .
In E. coli, there are more than ten well-characterized type II TA systems 1 . These include relE-relB, yafQ-dinJ, yoeB-yefM, hipA-hipB, yafO-yafN, hicA-hicB, higB-higA, ypjF-yfjZ, mqsR-mqsA, ymcE-gnsA and ydaT-ydaS. 10,[28][29][30][31][32][33][34][35][36][37]57 . The genomic location of each of these TA systems is indicated in Figure 1. It is of considerable practical importance to understand the dynamics of toxin-antitoxin systems and several plausible models for toxin-antitoxin dynamics and persister formation have been proposed (see, for example [38][39][40] and references therein). It is also important that the proposed model predictions are compared to, nowadays available, high-throughput data. In this paper, we present a minimal model for the description of toxin-antitoxin type II dynamics in E. coli. The basic characteristics of the minimal model is that it assumes: (a) regulation of toxin and antitoxin production rate by means of a negative feedback through DNA binding of the toxin-antitoxin complex (b) toxin induced growth rate modulation. The model's predictions are compared to the RNA-Seq gene expression data published in 41 (see Results and Discussion).
Toxin-antitoxin dynamics have been of interest to mathematical modelling for a long time. So far, the focus of research has been on the basic dynamical properties of TA modules 39,40,42,43 Figure 1. Genomic locations of the 11 TA modules studied in our investigation. Each gene is depicted as a radial dash on the circular chromosome. The location of each TA module is labelled with a red box. The OriC-Ter axis (from the origin of replication to the terminus of replication) is indicated for reference as a dashed blue line. environmental stimuli (e.g., 44 ), rather than the agreement with high-throughput data. For high-throughput data, in particular gene expression patterns, the dominant avenue of research has been to compare these patterns with large-scale regulatory networks or classes of regulatory mechanisms. In the case of bacterial gene regulation, successes have been understanding and experimentally confirming the role of small regulatory devices like feedforward loops 45,46 , the discovery of an interplay the regulatory network and chromosomal structure [47][48][49][50] and the organization of gene expression along the axis from the origin (OriC) to the terminus (Ter) of replication 50 .
Toxin-antitoxin systems are often embedded in an intricate network of regulatory processes 5 and part of functional regulatory modules 51 . There is evidence of collective behaviors arising from the interplay between toxin-antitoxin systems. Such a model of coupled toxin-antitoxin systems has for example been studied in 21 and in 44 . Simple ODE models of (type-II) toxin-antitoxin systems have for example been formulated in 21 with an emphasis on coupled systems and the spontaneous switching occurring in stochastic dynamics, in 40 , where conditional cooperativity of the RelBE system has been studied and its response to environmental stimuli (e.g., nutritional stress), in 53 , which contains a simplified system capable of excitable dynamics, as well as in 39 and 44 with a focus on bistability. For type-I TA systems, a mathematical model has been developed in 54 , offering insight in time scales involved.
Here we study the long-term dynamics of toxin-antitoxin pairs in time-resolved RNA-Seq data for E. coli. Our question is, whether the dynamics of all TA pairs in the data can be described by the same model, or whether qualitatively different models have to be assumed for the different TA modules. Figure 2 shows a schematic of the basic characteristics of the minimal model of type II toxin-antitoxin gene expression. Toxin T and antitoxin A are expressed by neighbouring genes. It is known 1, 39 that toxins are more stable than the antitoxins, thus, the latter have to be constantly expressed in order to neutralize the toxin effects. The toxin and antitoxin form a complex AT which inhibits toxin and antitoxin production. More complex toxin-antitoxin interaction (such as conditional cooperativity 39,40 or cooperation between multiple toxin-antitoxin systems 17 ) are not included in the minimal model. Moreover, the presence of toxin has an inhibitory effect on the cell growth. This last fact is found to be an essential characteristic of an acceptable minimal model.

Methods
We denote the concentration of the antitoxin A with the variable y 1 , that of the toxin T with y 2 and, finally, the concentration of the toxin-antitoxin complex AT with y 3 . The system of ordinary differential equations (ODEs) that describes the system is: Equation (3) is a standard chemical kinetics equation. We assume that the production rate of the complex y 3 is proportional to the product of the concentrations of y 1 and y 2 , thus the term k 3 y 1 y 2 where k 3 is the respective rate constant. We also assume that the complex degrades to its constituents A and T with a rate constant d 3 .
The inhibitory action of the AT complex is modelled through the inclusion of negative feedback terms such as . The existence of toxin T in the cell reduces all protein production and decreases protein dilution by decreasing cell growth. Thus, the toxin concentration will have an inhibitory impact on the production rates of toxin, antitoxin, and on the cellular growth rate. We introduce an inhibition factor 1/(b ′ m y 2 + 1) in Eqs. (1) -(2). The parameter b ′ m represents the redaction of protein expression due to the presence of toxin molecules. We also assume that growth inhibition will influence the toxin degradation rate, and we introduce a factor (b ′ c y 2 + 1) that modulates the toxin degradation rate in Eq. (2), while we assume that the degradation rate of the free antitoxin remains the same. This is in agreement with a recent finding from 55 that importantly, although free antitoxin is readily degraded in vivo, antitoxin bound to toxin is protected from proteolysis, preventing release of active toxin.
However, Eqs. (1) -(3), if one includes the unknown initial conditions for the quantities y 1 , y 2 , y 3 at t = 0, contain 13 adjustable parameters. Our aim is to estimate the model parameters using experimental RNA-Seq data obtained from 41 . These experimental data (10 data points for each toxin antitoxin pair) would render such an estimation problematic, since such a model is structurally unidentifiable 56 .
In order to reduce the number of adjustable parameters we rescale the unobserved variable y 3 by setting y 3 = (k ′ 2 /d 3 )z 3 and rescale the variables y 1 , y 2 by the same factor β = k ′ 2 , i.e. by setting y 1 = k ′ 2 z 2 and y 2 = k ′ 2 z 2 . Thus, we arrive at a system of ODEs for the rescaled variables z 1 , z 2 , z 3 which is: where the new kinetic constants are related the those in Eqs.
Moreover, we assume that z 1 and z 2 at time t = 0 are equal to zero and allow the unobserved complex concentration z 3 (0) to be equal to a constant c 0 which will be determined from the fitting of the solution of Eqs. (4) - (6) to the data. Henceforth, we will refer to the above model (Eqs. (4) -(6)) as the Z-model. The model is essentially a rescaled version of the model proposed in 39,40 with the additional assumption that the antitoxin bound to toxin is protected from proteolysis.
Our numerical investigations have shown that the Z-model (Eqs. (4) - (6)) is the simplest model able to represent the complete set of the experimental data that we have in our disposal with reasonable accuracy. Omission of any of the above basic ingredients of the model (e.g. setting b m and b c equal to zero) leads to plausible models, which may describe adequately the time evolution of the concentrations of some of the toxin-antitoxin pairs, but fail to describe the expression of the entire set. It is obvious to the reader that the Z-model and its variants that we examine in this manuscript are deterministic models. We will not deal with the important topic of investigating a stochastic variant of the Z-model through a Monte Carlo approach based on the Gillespie algorithm. Our modeling decision is based on the fact that the RNA Seq data that we will use to fit the model parameters are not single cell sequencing data. As one can see in the detailed description of the experimental data used in this study, each RNA seq "read" represents multi-cell averages on a time scale of hours. Of course for single cell RNA seq experiments a stochastic modelling approach would be more appropriate although admittedly much more difficult. There is, however, important progress in the direction of using stochastic models and the inference of parameter values from noisy data, see for example 23 .
For our analysis we used experimental RNA-Seq data obtained from 41 (GEO accession number: GSE65244). The RNA Seq data used here are for the wild-type(wt) strain and obtained after the culture growth in rich medium during the stationary phase. The system of Eqs. (4) -(6) was solved numerically with custom code written in Python using the scipy python module 24 . Fitting of the numerical solutions of the ODE's was performed as part of the code using the Nelder-Mead minimization algorithm as implemented in scipy. Since the task of performing fits for all TA pairs and all model variants is quite demanding the code was parallelized using the dask.distributed python module. All numerical simulations were performed on a workstation equiped with 2 Intel Xeon Gold 6140 Processors (72 cpu cores in total). Figure 3 shows the concentrations of toxin and antitoxin for 11 known TA pairs of E. coli as a function of time. Symbols represent experimental RNA-Seq data obtained from 41 (GEO accession number: GSE65244). The above list is exhaustive meaning that it includes all the TA pairs for which there are experimental measures in the dataset. All data have been rescaled (multiplied by the same constant c = 10 5 in order to avoid numerical errors during the fitting process). Lines are the numerical solutions of the ODE system, Eqs. (1) -(3). The kinetic constants of the system were estimated so that the weighted sum of the squared differences between the experimental data and the model predictions becomes minimum. We calculate weighted least squares since we have to fit two different experimental curves simultaneously whose y-axis values may differ considerably. Thus, we first calculate the mean values for each curve and then the weighted sum of the squared differences. Otherwise, curves with low mean values are practically ignored during the fitting process. Thus, the lines represent the "best" fit of the model to the data. We observe a very good agreement between the model predictions and the experimental data. As mentioned above, we assume that z 1 and z 2 at time t = 0 are equal to zero. This is a rather harsh, and probably unrealistic, condition to impose. If more data points were available the more natural and appropriate choice would be to use the RNA seq measurements of the earliest available timepoint as our initial conditions. This is indeed the approach we took in our analysis in Appendix B (Supplementary Materials).  , to the RNA-Seq data. Each box shows the "dispersion" of eleven values, one per toxin-antitoxin pair. We observe a wide distribution of parameter values across the different toxin-antitoxin pairs. This is rather common in biological systems, where the kinetic constants of various metabolic reactions can differ by several orders of magnitude. Therefore, the same underlying differential equations lead to quite different dynamics precisely due to the broad range of the kinetic constants. In Appendix A we include a detailed discussion of the estimated covariances and standard deviations of the fitting parameters (see also the attached files in supplementary materials). Dashed lines show the corresponding variable y 2 (t) for the toxin. Dotted lines show the corresponding variable y 3 (t) for the TA complex. We observe a variety of different dynamics, but interestingly enough in all cases the complex concentration z 3 seems to be lower than that of both the toxin and the antitoxin. For the majority of cases the antitoxin concentration is higher than that of the toxin. There are, however, exceptions, namely the relB-relE, mqsR-mqsA and the ymcE-gnsA pairs. The ydaT-ydaS pair also exhibits higher toxin expression for the most part of the observation time and only at the final stage the toxin level drop below that of the antitoxin. It is also quite intriguing that the Z-model predicts expression states where the toxin is constantly quite higher than the antitoxin (e.g. ymcE-gnsA) without resorting to the mechanism of conditional cooperativity 2, 39 , although it is quite well-established that certain TA pairs (e.g. the relB-relE pair) exhibit conditional cooperativity and, obviously, such effects are not accounted for in the Z-model.

Results
Next, we are interested in examining simpler versions of the proposed model and assessing their ability to describe the experimental data. We compare the Z-model to 7 simpler (i.e. with less adjustable parameters) variants, which we obtain from Eqs. (4) -(6) by forcing constraints on some of the constants, i.e. by fixing their numerical value or by setting them numerically equal to other constants. We describe these simpler variants below: • Model "s1=s2" is obtained by forcing the constants s 1 and s 2 to have the same numerical value.
• Model "s1=s2 no bm" is obtained by forcing the constants s 1 and s 2 to have the same numerical value and by dropping the b m constant, i.e. setting b m = 0.
• Model "s1=s2 no bc" is obtained forcing the constants s 1 and s 2 to have the same numerical value and by setting b c = 0.
• Model "s1=s2 no bm bc" is obtained by forcing the constants s 1 and s 2 to have the same numerical value and by setting both b m = 0 and b c = 0.
• Model "s1!=s2 no bm" is obtained by setting b m = 0. Note that now constants s 1 and s 2 are allowed to have different numerical values.
• Model "no s1 s2 bm bc" is the simplest variant and is obtained from the Z-model ODEs by setting s 1 = 1, Models, where the parameter b m is identically zero, do not take into account the reduction of protein expression due to the existence of toxin, while variants, where the parameter b c is identically zero, ignore the effect of growth inhibition. Fig 6 shows the minimum values of the objective function (i.e. the sum of weighted squared differences between model predictions and the experimental data) for all toxin-antitoxin pairs and for the 7 model variants described above.
The objective function of the Full Z-model is always lower than that of the variants, as expected. We should also mention that the algorithms (basinhopping in combination with a local Nelder-Mead algorithm) used for the minimization of the objective function are guaranteed to find local, not global, minima. Although we have performed a rather extensive search of the parameter space, there is always the chance that there are sets of parameters that will lead to lower values of the objective function than those reported here. We see that there are toxin-antitoxin pairs for which simpler variants are capable of fitting the data with results comparable to those of the Full Z-model. However, the Full Z-model is the appropriate choice if one wants to describe the expression of the entire set of TA pairs.
Since we want to compare models with different numbers of parameters, it might be plausible to examine two widely used model selection criteria, the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) for the Full Z-model and its seven variants. These are calculated as follows: where χ 2 stands for the sum of the squares of the residuals (i.e., the objective function discussed above), N is the number of data points (common for all model variants) and N v is the number of adjustable parameters for each model. each variant. The full Z-model has the highest value, i.e. N v = 10. The most appropriate model is considered to be the one with the lower AIC or BIC value since both these criteria penalize the a large N v number and reward a low objective function. Generally, the Bayesian information criterion is considered the most conservative of the two statistics. Figure 7 shows the AIC and BIC for the "collective" description of the TA gene expression set, i.e. when we describe the complete set of TA-pair with N = 10 * 11 = 110 data points and χ 2 is the sum of the objective functions of all the TA pairs.
Finally, it is helpful to compare the values of the constants that we obtained from the minimal ODE model for the different TA pairs. To this end we may view them as a "vector" characterizing the TA pair and we use an unsupervised learning method, namely a Principal Component Analysis (PCA), a statistical procedure that uses an orthogonal transformation to convert a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables called principal components 58 . PCA is routinely applied to experimental measurments directly for reasons of dimensionality reduction. Using PCA, however, to interpret the parameters of a deterministic ODE model consists a novel approach which has been recently used to interpret the parameters of a fractal kinetics SI model of Covid-19 spreading 59 . Figure 8 shows a plot of the two largest PCA components.
Typically in a PCA plot we try to identify clusters and perceive them as an indication of similar underlying causal behavior. For cluster identification, to avoid subjectivity, we applied a clustering identification algorithm i.e. DBSCAN with parameter eps = 0.8 60 . For DBSCAN the number of clusters is not predefined but decided by the algorithm. Here, the clustering algorithm has identified one cluster of 7 TA pairs, namelydinJ-yafQ, relB-relE, yafN-yafO, higA-higB, hipB-hipA, hicB-hicA, and mqsA-mqsR, which form a large central cluster, and four outliers i.e. the three pairs yefM-yoeB, ydaT-ydaS, and ymcE-gnsA, which have a negative PC2 component, and yfjZ-ypjF with relatively large PC1 and PC2 values.
In Table 1 we summarize this distinction between a main cluster and several outliers, together with the associated functional classification of the TA pairs. This distinction can serve as a starting point for comparing this statistical result with the wealth Figure 6. The values of the objective function for all toxin-antitoxin pairs and for 7 model variants. The Z-model defined by Eqs. (4) -(6) is marked with the label "Full" in the x-axis.
of biological information available for each of these toxin-antitoxin modules. For the toxin-antitoxin module hipB-hipA for example, the mode of action has been debated over the last years 62, 63 , but is still not clear 1 . The similarity of estimated parameters to higA-higB, hicB-hicA and other members of the main cluster may be seen as evidence of a functional classification of this TA system as RNA interferases and guide further attempts of functional elucidation, in particular a better understanding of superfamilies of type-II TA systems 64 .  Figure 7. Collective Akaike Information Criterion (AIC) and Baeysian Information Criterion (BIC) for the Full Z-model and its seven variants. The Full Z-model has clearly the lowest AIC and BIC among all variants studied.
Appendix B (Supplementary Materials) contains the results for another time-resolved gene expression data set, namely the data from 69 which are available at GEO (accession number: GSE131992).
In Appendix C (Supplementary Materilas), we present in tabular form the biological information relevant to the members of the clusters identified in Fig. 7 as obtained from The Universal Protein Resource (UniProt), a comprehensive resource for protein sequence and annotation data (https://www.uniprot.org).

Conclusions
We have proposed a minimal model that is able to capture the dynamics of toxin-antitoxin systems in E. coli and agrees with experimental high-throughput RNA-Seq data reasonably well. We find that a minimal acceptable model of toxin-antitoxin regulation should at least include a negative feedback loop through a TA pair formation and the effect of toxin induced growth modulation. Despite the obvious over-simplifications of the model, e.g. we study each toxin-antitoxin pair in isolation, and we do not account for the influence on cell growth due to the remaining toxin proteins, the model is able to replicate a variety of experimental curves.
With the availability of more time-resolved high-quality gene expression data, the description of time courses of systemic components with the help of simple mathematical models can provide an important instrument for the interpretation of such high-throughput data and thus bridge the gap between Theoretical Biology, Statistical Physics and Systems Biology.  Figure 8. Plot of the largest (PC1) vs. the second largest (PC2) principal components. A distinction between one main cluster and a set of outliers can be discerned: Central Cluster dinJ-yafQ, relB-relE, yafN-yafO, higA-higB, hipB-hipA, hicB-hicA, and mqsA-mqsR. Outliers yefM-yoeB, ydaT-ydaS, ymcE-gnsA, and yfjZ-ypjF.

Competing interests The authors declare no competing interests
Data availability The datasets analysed and the custom code used during the current study are available from the corresponding author on reasonable request.