Mathematical modelling of SigE regulatory network reveals new insights into bistability of mycobacterial stress response

Background The ability to rapidly adapt to adverse environmental conditions represents the key of success of many pathogens and, in particular, of Mycobacterium tuberculosis. Upon exposition to heat shock, antibiotics or other sources of stress, appropriate responses in terms of genes transcription and proteins activity are activated leading part of a genetically identical bacterial population to express a different phenotype, namely to develop persistence. When the stress response network is mathematically described by an ordinary differential equations model, development of persistence in the bacterial population is associated with bistability of the model, since different emerging phenotypes are represented by different stable steady states. Results In this work, we develop a mathematical model of SigE stress response network that incorporates interactions not considered in mathematical models currently available in the literature. We provide, through involved analytical computations, accurate approximations of the system’s nullclines, and exploit the obtained expressions to determine, in a reliable though computationally efficient way, the number of equilibrium points of the system. Conclusions Theoretical analysis and perturbation experiments point out the crucial role played by the degradation pathway involving RseA, the anti-sigma factor of SigE, for coexistence of two stable equilibria and the emergence of bistability. Our results also indicate that a fine control on RseA concentration is a necessary requirement in order for the system to exhibit bistability. Supplementary Information The online version contains supplementary material available at 10.1186/s12859-021-04372-5.

The following chemical reactions represent the starting point for the derivation of a set of Ordinary Differential Equations describing stress response of SigE regulatory network. Chemical reactions not included in the model proposed by Tiwari and collaborators in [1] are marked with a diamond.

Mathematical model
Under the assumptions of quasi-steady state approximation for mRNA dynamics, the functioning of SigE regulatory network is described by the following set of ODEs: dP P dt = k Pk ap P − k Pk ad P P + k 2 [ERP P ] − k 1 [ER]P P − k pdeg P P (6) dP dt = −k Pk ap P + k Pk ad P P + k 5 [ERP P ] + ν P − k pdeg P (7) dP 2 dt = f P 2 (E) − k 9 C 1 P 2 + k 10 C − k pdeg P 2 (14)

Module 1
Module 2 where and ν P is the production rate of P.

Bistability investigation through nullclines analysis
Nullclines analysis is a powerful tool to determine the number of equilibrium points of a system: the number of equilibrium points is given by the number of points where all of the nullclines intersect. A clear advantage of this technique is that it does not require running multiple simulations with initial conditions sampled over the state space of the system. On the other hand, when dealing with nonlinear systems of high (or relatively high) dimension, determining the number of solutions of the algebraic equations system is not obvious at all. To overcome this problem we decompose SigE regulatory network into two interconnected modules (see Figure 1): • In Module 1 SigE regulation takes place via the two-component system MprA/MprB: free sigma factor SigE (the input to Module 1) regulates transcription of the mprAB operon, and hence the total amount of proteins MprA and MprB. The two-component system then controls the ratio between phosphorylated and unphosphorylated portions of MprA and MprB: only phosphorylated MprA upregulates transcription of sigE, and thus controls the total amount of SigE, which represents the output of Module 1.
• In Module 2 the amount of free SigE is controlled by the anti-sigma factor RseA and by proteins ClpC1 Note that the output of Module 1 is the input of Module 2 and, viceversa, the output of Module 2 is the input to Module 1. In the following, we artificially break down the interconnection among the two modules, and separately analyze them to obtain input-output relationships. In particular, for each module, we first compute nullclines and manipulate them so as to obtain a suitable system of algebraic equations (Section 3.1 for the MprA/MprB regulation module, Section 3.3 for RseA/ClpC1P2 regulation module). Secondly, considering E as the independent variable, we rewrite the algebraic system in such a way that every nullcline is (directly or indirectly) a function of E (Section 3.2 for the MprA/MprB regulation module, Section 3.4 for RseA/ClpC1P2 regulation module). This allows to compute, for both modules, the total amount of SigE as a function of free SigE. Then, for a given a set of parameters, it is immediate to check whether the system admits a unique equilibrium or multiple equilibrium points.
3.1. Module 1: nullclines computation and derivation of the algebraic equations system 3.1.1. Two-component system MprA/MprB. Define the total amount of MprA protein as A T := A P + A, and notice that dA T dt = dA P dt + dA and The nullcline for A T is hence given by: Similarly, define B T := B P + B and notice that dB Recalling that A = A T − A P and B = B T − B P , the equilibrium conditions on A P and B P yields: The equilibrium condition dE T dt = 0 leads to the equation where we have set

Plugging equation (16) into equation (18) the following equation is obtained:
Solving the above equation for B P yields: Now, let us express B T − B P as Plugging equations (20) and (21) into the algebraic equation (17), after tedious computations we obtain a quadratic equation in the unknown A T of the form: where: Unfortunately, it is not easy to determine the sign of A(A P ), B(A P ) and C(A P ) and hence to decide which of the two solutions is the admissible one. However, plotting A T as a function of A P it is immediate to realize that the only solution satisfying the feasibility constraint A T ≥ A P is: Plugging equation (22) into the algebraic equation (15) and solving for E we obtain: and hence A P = φ −1 (E) (unfortunately, it is not possible to derive en explicit, closed form expression for φ −1 (E)). Finally, from equation (19) it follows that which represents the static input-output relationship of Module 1.
Define the total amount of PknB as P T := P P + P + [ERP P ], and notice that dP T dt = dP P dt + dP dt + d[ERP P ] dt = ν P − k pdeg P T , hence at equilibrium it holds P T = ν P k pdeg The equilibrium condition d[ERP P ] dt = 0 leads to the algebraic equation: where we exploited the fact that k pdeg k 2 , k 5 .
Similarly, the equilibrium condition dP P dt + d[ERP P ] dt = 0 yields: Again, since k pdeg k 5 , k Pk ad , we will replace the above equation with the following equilibrium condition: k 5 [ERP P ] = k Pk ap P − k Pk ad P P By substituting P = P T − P P − [ERP P ], and recalling that k 1 [ER]P P = (k 2 + k 5 )[ERP P ], the equilibrium condition results where we exploited the fact that k pdeg k 4 , k 5 .
The equilibrium condition for d[ER P ] dt + d[ER P C] dt = 0 yields: Hence, the nullclines result

Subnetwork involving ClpC1.
For the sake of analytical tractability we introduce the simplifying assumption that proteins ClpC1 and ClpP2 are described by identical differential equations (i.e., identical parameters) and exhibit the same initial concentration, so that C 1 (t) = P 2 (t) for every t ≥ 0. Define the total amount of ClpC1 protein as C T := C 1 +C + [ER P C], and note that: Recalling that C 1 = P 2 , the equilibrium condition d[ER P C] dt + dC dt = 0 leads to: where we exploited the fact that k pdeg k 10 . Recalling that C 1 = P 2 , the equilibrium condition for C is given by: where we first made use of equilibrium condition (26), and then exploited the fact that k pdeg k 7 , k 8 .

Module 2: nullclines solution in terms of E
The algebraic equations system (30)-(35) is composed of 6 algebraic equations in 7 unknowns (i.e., E, [ER], [ER P ], [ER P C], [ERP P ], C 1 and C), and hence it is underdetermined. Taking E as the independent variable, we now solve the algebraic system for the remaining unknowns. In other words, we express all of the nullclines (except for E-nullcline) as function of E.

[ER]-nullcline as a function of E.
From equations (30)-(31), the following quadratic equation in the unknown [ER] can be obtained: Note that, by Descartes' rule of signs, only one of the two solutions is positive, and hence the explicit solution for [ER] is given by: [C] + k pdeg Substituting (39) into equation (34) yields: [C] + k pdeg Secondly, plugging (33) and (34) into the definition of C T , we obtain: where we used the fact that k pdeg k 10 . Upon substituting [ER P C] with expression (39), the previous equation can be rewritten as [C] + k pdeg Since h(E,C) := C T − [ER P C] is non-negative, Descartes' rule of signs guarantees that only one solution of the previous quadratic equation is positive. We hence obtain the following expression for C 1 -nullcline as a function of h(E,C): We now compute C, C 1 and [ER P C] nullclines with the following iterative procedure: Then, h (0) (E,C) = lim C→+∞ h(E,C), namely: 1 resorting to equation (41), namely: Initialize C (k) resorting to equation (40) and C → +∞, namely: A1: Set k = k + 1 and update [ER P C] (k) , h (k) , C (k) 1 and C (k) as: The previous procedure provides [ER P C], C 1 and C nullclines as a function of E. An alternative procedure, which doesn't require the implementation of the iterative algorithm, is designed as follows: -Compute h(E) = lim C→+∞ h(E,C), namely: -Compute C 1 -nullcline as a function of E from equation (41), namely: -Compute C-nullcline as a function of E from equation (40) with C → +∞, namely: -Compute [ER P C]-nullcline from equation (39), namely: 3.4.4. ER P -nullcline as a function of E. It follows from equations (32) and (35) that Solving the above equation for [ER P ] we obtain [ER P ]-nullcline as a function of E, namely:

Total amount of SigE as a function of E.
Recalling the definition of E T , the static relationship among E and E T can finally be computed.

Including RseA dynamics within the model
So far we have assumed that RseA concentration is constant and equal to R T . We now remove such an assumption, and consider RseA basal production and degradation and ClpC1P2-mediated degradation of phosphorylated RseA. To this aim, the following chemical reactions need to be taken into account: RseA binding to SigE: complex formation and dissociation: ClpC1P2-mediated degradation of phosphorylated RseA: From the previous set of chemical reactions, the following differential equation describing dynamic evolution of RseA concentration can be readily obtained: The equilibrium condition dR dt = 0 yields By combining the previous equation with Eqn. (24) we get where we have set R max = ν R δ R , and hence Now, rewrite equation (30) as Plugging equation (45) into equation (44), we obtain Rewrite equation (31) (recalling that now R is not constant) as As in the case with constant RseA, Descartes' rule of signs ensures that only one of the two solutions is positive, and hence the explicit solution for [ER] is given by: Notice that when R max = R T and δ R → +∞, the above quadratic equation in [ER] reduces to the analogous quadratic equation of the case with constant RseA (R ≡ R T ).