Figure 2 illustrates in a simplified and schematic way how the regulatory network of Fig. 1 has been described by Tiwari and colleagues (left) and in the present work (right). The mathematical model proposed by Tiwari and collaborators provides a dynamic description of the subnetwork within the orange rectangle and describes the remaining part of the network (blue rectangle) through the action of a fixed, constant amount of RseA. The mathematical model we are proposing makes a considerable step further by providing a dynamic description also of the subnetwork within the blue rectangle. In particular, our model: (i) describes stress response initiation by both MprB and PknB; (ii) considers regulation in free SigE concentration by ClpC1P2; (iii) accounts for RseA degradation.
The main findings of our work are, firstly, the development of a new model of SigE regulatory network (more accurate with respect to state of the art models) and, secondly, elucidation of the critical role played by RseA degradation pathway for the emergence of bistability and coexistence of two stable equilibria. To disentangle the multiple feedback loops involved within the network, we cast the whole network of Fig. 1 into a block structure with feedback interconnection and adopt an approach based on nullclines analysis to determine, in a reliable though computationally efficient way, the number of equilibrium points of the system.
Mathematical model
Ordinary differential equations describing the behaviour of MprA/MprB twocomponent system are borrowed from [11] and reported here for convenience. We adopt capital letters A and B for state variables representing the concentration level of proteins MprA and MprB, respectively, and we use the subscript P to denote the corresponding phosphorylated forms. As in [11], dynamic equations take into account exogenous phosphorylation (dephosphorylation) of MprA (phosphorylated MprA), autophosphorylation (autodephosphorylation) of MprB (phosphorylated MprB), phosphotransfer from phosphorylated MprB to MprA, phosphatase activity of MprB, upregulation of MprA and MprB synthesis by phosphorylated MprA and SigE (denoted with the capital letter E), and proteins degradation. The resulting differential equations are given by (see Additional file 1  Supplementary Material for further details):
$$\begin{aligned} \frac{dA_P}{dt}&= \frac{k_t}{K_T}AB_P  \frac{k_p}{K_P}A_PB + k_{exp}A  k_{exd}A_P  k_{pdeg}A_P \end{aligned}$$
(1)
$$\begin{aligned} \frac{dB_P}{dt}&= k_{ap}B  k_{ad}B_P  \frac{k_t}{K_T}AB_P  k_{pdeg}B_P \end{aligned}$$
(2)
$$\begin{aligned} \frac{dA}{dt}&= \beta _1 \frac{\left( 1 + f_1\frac{A_P^2}{K_1}\right) }{\left( 1 + \frac{A_P^2}{K_1}\right) } + \beta _2 \frac{\left( 1 + f_2\frac{E}{K_2}\right) }{\left( 1 + \frac{E}{K_2}\right) } + \frac{k_p}{K_P}A_PB +\nonumber \\&\qquad  \frac{k_t}{K_T}AB_P + k_{exd}A_P  k_{exp}A  k_{pdeg}A \end{aligned}$$
(3)
$$\begin{aligned} \frac{dB}{dt}&= \lambda \beta _1 \frac{\left( 1 + f_1\frac{A_P^2}{K_1}\right) }{\left( 1 + \frac{A_P^2}{K_1}\right) } + \lambda \beta _2 \frac{\left( 1 + f_2\frac{E}{K_2}\right) }{\left( 1 + \frac{E}{K_2}\right) } + \frac{k_t}{K_T}AB_P \nonumber \\&\qquad + k_{ad}B_P  k_{ap}B  k_{pdeg}B \end{aligned}$$
(4)
The amount of free SigE (denoted by E) is described by differential equation (5), which takes into account: positive regulation on the synthesis of SigE by phosphorylated MprA, binding by RseA to form the complex denoted by [ER] and dissociation of the [ER] complex (here, \(R_T\) denotes the total amount of RseA), ClpC1P2mediated dissociation of the complex \([ER^PC]\) formed by SigE and phosphorylated RseA (\(R^P\)), protein’s degradation:
$$\begin{aligned} \frac{dE}{dt} = \beta _3 \frac{\left( 1 + f_3\frac{A_P^2}{K_1}\right) }{\left( 1 + \frac{A_P^2}{K_1}\right) }  k_3 E R_T + k_4 [ER] + k_8[ER^PC]  k_{pdeg}E \end{aligned}$$
(5)
The concentrations of phosphorylated PknB and its unphosphorylated form (denoted by \(P^P\) and P, respectively) are described by the following differential equations:
$$\begin{aligned} \frac{dP^P}{dt}&= k_{ap}^{Pk}P  k_{ad}^{Pk}P^P + k_2 [ERP^P]  k_1 [ER] P^P  k_{pdeg}P^P \end{aligned}$$
(6)
$$\begin{aligned} \frac{dP}{dt}&=  k_{ap}^{Pk}P + k_{ad}^{Pk}P^P + k_5 [ERP^P] + \nu _P  k_{pdeg}P \end{aligned}$$
(7)
The above equations, directly derived from chemical reactions, assume constant PknB synthesis (at rate \(\nu _P\)) and account for autophosphorylation (autodephosphorylation) of PknB (phosphorylated PknB), formation and dissociation of the proteins complex degrading RseA (see Additional file 1  Supplementary Material for further details).
Dynamics of the complex formed by SigE and RseA (denoted by [ER]), and of the other complexes involved in RseA degradation pathway (specifically, \([ERP^P]\), \([ER^P]\), \([ER^PC]\)) are described by the following differential equations (obtained from the corresponding chemical reactions, see Additional file 1  Supplementary
Material for further details):
$$\begin{aligned} \frac{d[ERP^P]}{dt}&= k_1 [ER] P^P  k_2 [ERP^P]  k_5 [ERP^P]  k_{pdeg}[ERP^P] \end{aligned}$$
(8)
$$\begin{aligned} \frac{d[ER]}{dt}&= k_2 [ERP^P]  k_1 [ER] P^P + k_3 E R_T  k_4 [ER]  k_{pdeg}[ER] \end{aligned}$$
(9)
$$\begin{aligned} \frac{d[ER^P]}{dt}&= k_5 [ERP^P]  k_6 [ER^P] C + k_7 [ER^PC]  k_{pdeg}[ER^P] \end{aligned}$$
(10)
$$\begin{aligned} \frac{d[ER^PC]}{dt}&= k_6 [ER^P] C  k_7 [ER^PC]  k_8 [ER^PC]  k_{pdeg} [ER^PC] \end{aligned}$$
(11)
Finally, dynamic description of proteins ClpC1, ClpP2 and of the complex ClpC1P2 (denoted by \(C_1\), \(P_2\) and C, respectively) is obtained by taking into account positive regulation on the synthesis of ClpC1 and ClpP2 by SigE, formation and dissociation of the ClpC1P2 complex and of the other complexes (i.e., \([ER^P]\) and \([ER^PC]\)) involved in RseA proteolitic degradation. The resulting differential equations are given by:
$$\begin{aligned} \frac{dC}{dt}&=  k_6 [ER^P] C + k_7 [ER^PC] + k_8 [ER^PC] + k_9 C_1 P_2  k_{10}C  k_{pdeg}C \end{aligned}$$
(12)
$$\begin{aligned} \frac{dC_1}{dt}&= f_{C_1}(E)  k_9 C_1 P_2 + k_{10} C  k_{pdeg}C_1 \end{aligned}$$
(13)
$$\begin{aligned} \frac{dP_2}{dt}&= f_{P_2}(E)  k_9 C_1 P_2 + k_{10} C  k_{pdeg}P_2 \end{aligned}$$
(14)
where \(f_{C_1}(E)\) and \(f_{P_2}(E)\) denote Hill equations, namely
$$\begin{aligned} f_{C_1}(E)&:= \beta _{C1} \frac{\left( 1 + f_{C1}\frac{E}{K_{C1}}\right) }{\left( 1 + \frac{E}{K_{C1}}\right) }\\ f_{P_2}(E)&:= \beta _{P2} \frac{\left( 1 + f_{P2}\frac{E}{K_{P2}}\right) }{\left( 1 + \frac{E}{K_{P2}}\right) } \end{aligned}$$
Differential equations (1)–(14) describe the functioning of the whole SigE regulatory network under the simplifying assumption, borrowed from Tiwari’s model [11], that RseA concentration is constant. RseA is indeed a fixed parameter of the network, denoted by \(R_T\) and appearing in differential equations (5) and (9).
In contrast with the assumption of constant RseA concentration, experimental data show that, in response to surface stress, rseA mRNA levels remain stable while RseA concentration decreases [14, 18]. To take into account experimental evidence of nonconstant RseA concentration, we further develop an alternative model of SigE regulatory network by considering RseA as a state variable endowed with a proper dynamics. Specifically, we describe RseA concentration by the following differential equation:
$$\begin{aligned} \frac{dR}{dt} = \nu _R  k_3 E R + k_4 [ER]  \delta _R R \end{aligned}$$
(15)
The above equation accounts for RseA production and degradation (at constant rates \(\nu _R\) and \(\delta _R\), respectively), formation and dissociation of the complex formed by SigE and RseA. Most importantly, since RseA is proteoliticallty degraded by ClpC1P2, RseA is not retrieved after dissociation of the complex \([ER^PC]\) (differently from what happens to SigE and ClpC1P2, see the terms \(+k_8[ER^PC]\) appearing in Eqs. (5) and (12)). The alternative mathematical model of SigE regulatory network, which takes into account RseA dynamics, is hence composed of differential equation (15) together with differential equations (1)–(14), upon substitution of parameter \(R_T\) with the state variable R in Eqs. (5) and (9).
Parameters setting
Parameters that define regulatory interactions included in the mathematical model by Tiwari and collaborators have been set in accordance with [11]. Numerical values of these parameters is taken from Table S3 in [11] with two exceptions. First, in our model, to ensure perfect balance between exogenous phosphorylation and dephosphorylation fluxes, the exogenous dephopshorylation rate constant of phosphorylated MprA (i.e., parameter \(k_{exd}\)) equals the exogenous phopshorylation rate constant of MprA (i.e., parameter \(k_{exp}\)). This choice is consistent with the analysis carried out in Section 2.5 of [11] and demonstrating bistability of mycobacterial stressresponse network (bifurcation diagrams of Figure 5 in [11] have indeed been obtained for \(k_{exd} = k_{exp}\)). Secondly, for the model with constant RseA concentration, parameter \(R_T\) is slightly increased with respect to the value reported in Table S3 of [11] (however, it still belongs to the bistability region reported in Figure S2 of [11]). This modification is justified by the fact that, with respect to the model in [11], we are considering additional proteins complexes in which RseA is present.
Unfortunately, for parameters defining the regulatory interactions mediated by PknB, ClpC1 and ClpP2 (not included in [11]), estimates based on experimental data are not available in the literature. When similar interactions (namely, similar chemical reactions) are described in [11], then corresponding parameters take similar values, see, e.g., the basal transcription rate and the amplification gain of proteins ClpC1, ClpP2. When this reasoning is not applicable, namely, when a parameter playing a similar role in the model by Tiwari cannot be found, then the order of magnitude is set either to \(10^{3}\) (like parameter \(\frac{k_p}{K_P}\) in [11]) or to \(10^{2}\) (like parameter \(\frac{k_t}{K_T}\) in [11]).
Numerical values of all model parameters, either borrowed from [11] or set according to the previous reasoning, can be found in the Additional file 2  Table I. In addition, to mitigate the effects due to uncertainties on model parameters for which experimental data are not available, numerical experiments have been performed (see “Parameters perturbation experiments” section) in which parameters whose value is not borrowed from [11] are randomly perturbed, either one at a time or multiple parameters at the same time.
Bistability investigation through nullclines analysis
The whole SigE regulatory network can be seen as the feedback interconnection of two modules controlling the overall SigE concentration and the amount of total SigE that is free and hence functionally active. This interpretation is graphically illustrated in the block diagram of Fig. 3. In Module 1 SigE regulation takes place via the twocomponent 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 twocomponent system then controls the ratio between phosphorylated and unphosphorylated portions of MprA and MprB. Only phosphorylated MprA upregulates transcription of sigE gene, and thus controls the total amount of SigE protein, which represents the output of Module 1. In Module 2 the amount of free SigE is controlled by the antisigma factor RseA and by proteins ClpC1 and ClpP2. SigE, playing the role of input to Module 2, is partly bound by RseA and degraded by the protein complex ClpC1P2. The remaining free, and hence functionally active, amount of SigE represents the output of Module 2. Clearly, the feedback interconnection is such 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 order to underpin the mechanisms leading to coexistence of two stable equilibrium states, we artificially break down the feedback interconnection and derive input–output relationships of each module separately. Via nullclines computation and involved manipulations of their nonlinear algebraic expressions, we provide (i) for Module 1, exact implicit expression of total SigE concentration as a function of free SigE, and (ii) for Module 2, approximated implicit expressions of total SigE concentration as a function of free SigE. We remark that the obtained input–output relationship for Module 2 is not the exact input–output relationship but an approximation of the exact relationship. Indeed, the high nonlinearity of the differential equations makes exact computation of the input–output relationship impracticable. However, an accurate approximation of the true input–output relationship can be derived by exploiting time scale separation between proteins degradation and protein complexes formation and dissociation. Then, the number of intersection points between the input–output relationship from Module 1 and the approximated input–output relationship from Module 2 is precisely the number of equilibrium points of the system. Notice that an advantage of our approach based on nullclines analysis is the fact that it allows, for a given set of parameters, to immediately check whether the system admits a unique equilibrium or multiple equilibrium points, without the need to run hundreds of simulations starting from initial conditions that explore, with sufficiently dense sampling, the state space.
Clearly, the feedback architecture of Fig. 3 is maintained independently of the regulations included within Module 2. In Tiwari’s model [11] a fixed, constant amount of RseA controls the level of free SigE and represents the unique regulatory mechanism accounted for by Module 2. On the contrary, in our model PknB and ClpC1P2 mediate additional regulations on SigE, and the antisigma factor RseA is endowed with a proper dynamics. The choice of relaxing the assumption of constant RseA is justified by the data showing that in response to surface stress, while rseA mRNA levels remain stable [24], RseA is proteolitically degraded by ClpC1P2 after its phosphorylation by PknB clearly indicating that in these conditions RseA concentration decreases [14, 18]. On the other hand, from a mathematical point of view, disregarding RseA dynamics and considering it as a fixed parameter means that the positive feedback loop implemented by ClpC1P2 degrading phosphorylated RseA has no actual effect on RseA. For these reasons, we decided to remove the assumption of constant RseA and instead to consider it as a state variable endowed with a proper dynamics. The ordinary differential equation describing RseA concentration accounts for its basal production and degradation (at constant rates), and additional proteolitic degradation by ClpC1P2. Nullclines reported in Fig. 4 clearly show that, under the assumption of dynamic RseA, the closedloop system exhibits three distinct equilibria, two of which are asymptotically stable (while the third one is necessarily unstable). Remarkably, when the level of RseA is kept constant and the same values (borrowed, when applicable, from [11]) for the remaining model parameters are assumed, the system exhibits a unique equilibrium point^{Footnote 2}, as reported in Fig. 5.
It is worth mentioning that high nonlinearity of the mathematical model makes exact analytical derivation of the system’s nullclines impracticable. As a consequence, the curves of Figs. 4 and 5 have been obtained under biologically reasonable approximations (see “Methods” section and Additional file 1  Supplentary Material for further details on the derivations). Validity and accuracy of the approximated nullclines is however testified by the exact equilibrium points reported with black circles in the figures and obtained by numerical simulations of the model with random initial conditions (indicated with black crosses in the figures). In fact, the equilibrium point obtained via numerical simulation lay at the intersection of the two approximated nullclines.
Figures 4 and 5 together point out the importance of taking into account RseA dynamics and, in particular, its degradation pathway implementing a positive feedback loop on SigE regulation. Looking at the curves^{Footnote 3} we realize that, in order for the equilibrium point corresponding to higher SigE levels to appear, the blue nullcline needs to undergo a slowdown so as to form, in logarithmic scale, a sort of “plateau”. This is exactly the role of RseA degradation pathway. On the other hand, in order for the equilibrium point corresponding to lower SigE levels to appear, the total amount of SigE (on the vertical axis) needs to be highly sensitive to small variations in free SigE concentrations (on the horizontal axis), so as for the blue nullcline to exhibit a high slope (in logarithmic scale) in its initial segment, i.e., for small amount of free SigE. Our results suggest that a fine regulation of RseA concentration levels is critical for the coexistence of two stable equilibrium points.
Parameters’ perturbation experiments
Nullclines analysis and model simulations show the need to include RseA degradation pathway into SigE regulatory network’s model in order to capture coexistence of two steady states as experimentally observed. To further investigate the role played by interactions newly introduced into our mathematical model, we perform some numerical experiments in which we perturb model’s parameters and quantitatively evaluate robustness of bistability. When intersection points among nullclines are spaced and clearly distinct, it is reasonable to expect that bistability will be retained in spite of perturbations on the parameters. We hence quantify robustness of bistability by computing the projection on the \(E_{free}\)axes of the distance between intersection points corresponding to stable equilibria. Referring to the nullclines plot in Fig. 4 and denoting by \(E_1\) and \(E_2\) the abscissa of the stable equilibrium points, we introduce the scalar distance function \({\mathcal {D}}(\mathbf{p }) := E_2  E_1\), where \(\mathbf{p }\) is the vector of parameters. Clearly, when monostability arises, such a distance vanishes, i.e., \({\mathcal {D}}(\mathbf{p }) = 0\), since the two stable equilibria coincide and \(E_1 = E_2\). We denote by \(\bar{{\mathcal {D}}}\) the distance computed with nominal values for the system’s parameters. For an \(\varepsilon\)percentage variation on parameters selected through the vector \(\mathbf{v }\) whose entries are either \(1\), 0, or 1, we define the following measure of bistability robustness:
$$\begin{aligned} {\mathcal {R}}_\varepsilon (\mathbf{v }) := \frac{{\mathcal {D}}(\mathbf{p } + \varepsilon \mathbf{v })  \bar{{\mathcal {D}}}}{\bar{{\mathcal {D}}}} * \frac{100}{\varepsilon } \end{aligned}$$
(16)
Note that efficient computation of the above metric is made possible by the derivation of the (implicit or explicit) solutions to the system’s nullclines.
Since our focus is on the regulations mediated by PknB, ClpC1P2 and RseA degradation pathway, in the following only parameters defining these interactions are perturbed, while model parameters that correspond to analogous parameters in Tiwari’s model [11] and whose value has been borrowed from, are kept constant.
Robustness experiment 1
In the first robustness experiment we aim to investigate the effects of local, i.e., with \(\varepsilon\) small (\(\varepsilon = 2.5\%\)) perturbations on the system’s parameters. To this aim, we vary model parameters one at a time, namely we set \(\mathbf{v } = \mathbf{e }_i\) where \(\mathbf{e }_i\) is the canonical vector with the ith entry equal to 1 and all other entries equal to 0, and we compute the robustness measure \({\mathcal {R}}_{2.5}(\mathbf{e }_i)\). Numerical results obtained by performing such an experiment on the mathematical model with dynamic (i.e., nonconstant) RseA concentration are reported in Fig. 6. It results that robustness of bistability is most sensitive to parameters controlling basal and enhanced production of protein ClpC1 (i.e., parameters \(\beta _{C1}\) and \(f_{C1}\)) and to RseA degradation rate (i.e., parameter \(\delta _R\)). Indeed, a slight increase in ClpC1 production (either basal or enhanced by SigE) substantially increases the robustness metric \({\mathcal R}_{2.5}\). Conversely, an increase in RseA degradation rate causes a decrease in the distance metric \(\mathcal D\).
Robustness experiment 2
The second robustness experiment is designed to test the effects of random perturbations that (i) have large absolute value (so as to explore regions of the parameter space farther from approximately linear effects explored in Robustness experiment 1), and (ii) involve more that one parameter at a time. For fixed large value of the perturbation \(\varepsilon\) the procedure outlined in Algorithm 1 is carried out for the mathematical model including RseA degradation pathway.
Algorithm 1

0.
Fix the number N of perturbations tests and set \(i = 1\).

1.
Generate a random vector \(\mathbf{v }\) whose entries are integers drawn from the discrete uniform distribution on the set \(\{1, 0,1\}\).

2.
Compute the distance metric \({\mathcal {D}}(\mathbf{p } + \varepsilon \mathbf{v })\).

3.
If \(i<N\), repeat from step 1.

4.
Plot a histogram of the normalized distance between stable equilibria^{Footnote 4}, i.e., \({\mathcal {D}}(\mathbf{p } + \varepsilon \mathbf{v }) / \bar{{\mathcal {D}}}\).
The resulting histograms provide an estimate of the relative probability of the distance \({\mathcal {D}}(\mathbf{p } + \varepsilon \mathbf{v })\). Figure 7 illustrates results for \(N=10^4\) and percentage variations \(\varepsilon = 10\%\) (Fig. 7a), \(\varepsilon = 15\%\) (Fig. 7b) and \(\varepsilon = 25\%\) (Fig. 7c). It is worth highlighting that bistability is always retained with random perturbations of \(\pm 10\%\) with respect to nominal values; when very large perturbations are considered (\(\pm 25\%\)), bistability is lost on \(22.54\%\) of iterations. These results indicate that our mathematical model is quite robust to parameter variations.
To identify parameters more closely connected to enhanced bistability, we retrieved from perturbation experiments with \(\varepsilon = 10\%\) and \(\varepsilon = 15\%\) all perturbation vectors \(\mathbf{v }\) corresponding to normalized distance metric larger than 1.1. Bar plots summarizing their sign patterns are reported in Figs. 8 and 9. Interestingly, observations resulting from perturbation experiment 1 are confirmed: increased bistability is associated with increased ClpC1 production (either basal or enhanced by SigE) as well as with decreased RseA degradation rate.
To unravel parameters whose variation leads to a loss of bistability, we retrieved from perturbation experiment with \(\varepsilon = 25\%\) all perturbation vectors \(\mathbf{v }\) corresponding to monostability, whose sign patterns are reported in Fig. 10. Consistently with previous observations, decreases in parameters \(\beta _{C1}\) and \(f_{C1}\) (which regulate ClpC1 production) and increases in parameter \(\delta _R\) are frequently reported when bistability is lost. However, RseA binding rate to SigE (i.e., parameter \(k_3\)) appears to predominantly control the loss of bistability, its decrease being associated with monostability.