Understanding key features of bacterial restriction-modification systems through quantitative modeling

Background Restriction-modification (R-M) systems are rudimentary bacterial immune systems. The main components include restriction enzyme (R), which cuts specific unmethylated DNA sequences, and the methyltransferase (M), which protects the same DNA sequences. The expression of R-M system components is considered to be tightly regulated, to ensure successful establishment in a naïve bacterial host. R-M systems are organized in different architectures (convergent or divergent) and are characterized by different features, i.e. binding cooperativities, dissociation constants of dimerization, translation rates, which ensure this tight regulation. It has been proposed that R-M systems should exhibit certain dynamical properties during the system establishment, such as: i) a delayed expression of R with respect to M, ii) fast transition of R from “OFF” to “ON” state, iii) increased stability of the toxic molecule (R) steady-state levels. It is however unclear how different R-M system features and architectures ensure these dynamical properties, particularly since it is hard to address this question experimentally. Results To understand design of different R-M systems, we computationally analyze two R-M systems, representative of the subset controlled by small regulators called ‘C proteins’, and differing in having convergent or divergent promoter architecture. We show that, in the convergent system, abolishing any of the characteristic system features adversely affects the dynamical properties outlined above. Moreover, an extreme binding cooperativity, accompanied by a very high dissociation constant of dimerization, observed in the convergent system, but absent from other R-M systems, can be explained in terms of the same properties. Furthermore, we develop the first theoretical model for dynamics of a divergent R-M system, which does not share any of the convergent system features, but has overlapping promoters. We show that i) the system dynamics exhibits the same three dynamical properties, ii) introducing any of the convergent system features to the divergent system actually diminishes these properties. Conclusions Our results suggest that different R-M architectures and features may be understood in terms of constraints imposed by few simple dynamical properties of the system, providing a unifying framework for understanding these seemingly diverse systems. We also provided predictions for the perturbed R-M systems dynamics, which may in future be tested through increasingly available experimental techniques, such as re-engineering R-M systems and single-cell experiments. Electronic supplementary material The online version of this article (doi:10.1186/s12918-016-0377-x) contains supplementary material, which is available to authorized users.

(1.2) -RNAP binding to the CR promoter forming a RNAP DNA  complex; (1.3) -D binding to the distal (high affinity) binding site, forming a D DNA  complex; Note that this configuration was neglected in [1], due to a very low probability of detecting a single dimer bound in the wt AhdI. We here include this configuration, since as we perturb (in-silico mutate) the system this configuration can become important. (1.4)a second D recruited to the proximal (low affinity) binding site, forming a C tetramer on the promoter (T DNA  complex); (1.5) -RNAP recruited to its binding site forming a D DNA RNAP  complex.
In equilibrium, the above reactions lead to the following relationships between the relevant concentrations and the dissociation constants (note that the square brackets refers to the reactant concentration): 2  to the free energy of interaction between the C dimer bound to the distal binding site and RNA polymerase (RNAP), with a proportionality constant k (in units of concentration [2]).
The model of transcription regulation is based on the assumption by Shea-Ackers [3], by which the promoter transcription activity is proportional to its equilibrium occupancy by RNAP: Here is a proportionality constant, while  . We will subsequently explore the system dynamics when 1 K is perturbed (decreased), which leads to both monomers and dimers present in the solution (see below).

Thermodynamical model of M promoter transcription regulation
Methyltransferase (M) methylates specific sites in M promoter, which overlap RNAP binding site, thereby repressing transcription of its own genethis mechanism is also found in other R-M systems (reviewed in [4]). We model the control of M expression as an equilibrium negative autoregulation, as described below.
Transcription regulation of M promoter is then described by the following set of reactions (see Figure 2B): The reactions above lead to the following relation between the dissociation constants and the reactant concentrations: We again assume that transcription activity is proportional to the promoter occupancy by RNAP: where  is a proportionality constant, and i) 1 is the statistical weight for the empty promoter configuration, ii) In both equations, the first term on the right-hand side represents transcript synthesis by transcription of the appropriate genes from their promoters, while the second term represents transcript decay by degradation. Protein dynamics is modeled by a similar set of differential equations shown below, in which the first term on the right-hand side describes protein synthesis by transcript translation, while the second term describes protein decay by degradation (for the notation and the constant values, see the Table 1  C and R are produced by translation of the same CR operon transcript, and as their degradation rates are assumed to be the same, ( / ) ,

RC R k k C
 which reduces the number of the differential 5 equations to be solved. Note that C transcript is leaderless, so that following [1] 5 R M C k k k  (see Table 1 below).
The system of differential equations (1.24)-(1.26) and (1.28) is numerically solved in MATLAB, using Runge-Kutta method with the initial conditions set to zero. The model of AhdI regulation is next used to assess the changes of the system dynamics, upon gradually abolishing the characteristic features of the system regulation.  Table 1: Values of the constants are the same as in [1], which were shown to be able to explain the equilibrium measurements for AhdI wild type system. The transcription activities and the transcript decay rates are given in units of transcripts per minute, while the translation initiation rates and the protein decay rates are given in proteins per minute.

Decreasing the dissociation constant of dimerization
We vary the equilibrium dissociation constant of dimerization 1 K , which changes the balance of monomers and dimers in the solution. The total concentration of the synthesized C protein is:  K , which we estimate as the lowest 1 K in which the monomer dimer balance (Eq. (1.30)) leads to the R dynamics that does not notably depart from the one in wt system (estimated as 1max~4 000 nM K  1.12)).

Changing cooperativity in the dimer binding
Cooperativity in the dimer binding is decreased by varying  (see Eq. (1.16)), which translates to varying the parameter c that is related to the establishment of the repressing tetramer complex. As it is estimated that 3000  in the wt AhdI system [5], we varied c in the range from 3000 where positive 2  indicates the stable equilibrium. To quantify the changes in the stability of R steady-state levels, we assess how 2  is changed upon introducing a perturbation in the system (e.g. as 1 K is decreased). Note that the decrease of 2  indicates the decrease in the stability of the R steady-state levels, and vice versa.

The effect on the system switch-like behavior
We here assess how fast the system makes a transition from "OFF" to "ON" state, i.e. from the initial small R values to the values approaching the system equilibrium. For quantifying this effect, we determine the maximum rate of   Rt change during the system evolution. As   Rt curves have sigmoidal shape, this maximal rate corresponds to the time interval in which R makes an approximately linear increase from small to high concentrations.

The effect on the delay of R with respect to M
To quantify the effect of changing the system features on the delay in expression of R with respect to M in the initial interval of the protein synthesis (taken as the first 10 minutes), we define the following parameter (relative delay): Rt correspond to the perturbed and the wt system, respectively. Note that, for the wt system, the relative delay is 1, while values larger than 1 indicate an increase in R delay, and vice-verse.