Inhibitors Alter the Stochasticity of Regulatory Proteins to Force Cells to Switch to the Other State in the Bistable System

The cellular behaviors under the control of genetic circuits are subject to stochastic fluctuations, or noise. The stochasticity in gene regulation, far from a nuisance, has been gradually appreciated for its unusual function in cellular activities. In this work, with Chemical Master Equation (CME), we discovered that the addition of inhibitors altered the stochasticity of regulatory proteins. For a bistable system of a mutually inhibitory network, such a change of noise led to the migration of cells in the bimodal distribution. We proposed that the consumption of regulatory protein caused by the addition of inhibitor is not the only reason for pushing cells to the specific state; the change of the intracellular stochasticity is also the main cause for the redistribution. For the level of the inhibitor capable of driving 99% of cells, if there is no consumption of regulatory protein, 88% of cells were guided to the specific state. It implied that cells were pushed, by the inhibitor, to the specific state due to the change of stochasticity.

of the mean first passage time 28,29 . Recently, instead of altering the concentration of the critical intracellular variable, it has been reported that the manipulation of internal stochasticity by tuning the rates of transcription and translation can be utilized to cause the redistribution and to lead cells to the desired state in bimodal distribution 30 . This finding provided a new strategy to control the phenotypic variation through engineered gene circuits. In this work, we aim to examine how the change of the stochasticity, caused by the inhibitor, affects the bimodal distribution.

Models
The gene regulatory network and the deterministic model. The genetic toggle switches are composed of two genes encoding regulatory proteins P1 and P2 which inhibit the expression of gene 2 and gene 1, respectively (Fig. 1). More precisely, the trimer of protein P2 represses gene 1, and the trimer of P1 represses gene 2. Such a reaction network, including two opposing fate-determining proteins, is one of the canonical motifs of bistable switches. This system also accounts for two inhibitors, the peptide Y acted on trimer of protein 1 and the peptide X on trimer of protein 2. Accordingly, the reactions of the system are listed in Supplementary Table S1. Following the mass action law, the deterministic model (Supplementary Table S2) was formulated in agreement with the reactions. The nomenclature and the values of parameters, which are based on literature 18,22,[30][31][32][33][34][35][36][37] , are listed in Supplementary Tables S3 and S4, respectively. The steady state responses were obtained by functions fsolve or solve of the symbolic toolbox in Matlab.
The Stochastic Model. Based on the reaction network shown in Fig. 1, we formulated the Chemical Master Equation (CME). Following the Stochastic Simulation Algorithm (SSA) 38 , we conducted the numerical analysis with the software SynBioSS 39 . Each simulation is composed of ten thousand trajectories (trials) and each trajectory is obtained by calculating the system to ninety thousand seconds. It was sampled only at the end point. The parameters of nominal values are listed in Supplementary Table S4. The cellular volume is of 10 −15 liter. The initial conditions of proteins are given as the values of the unstable steady-state obtained from the deterministic model and the DNA starts with the active configuration. For other variables, they are zeros at beginning. Note that the initial conditions have no impact on final results as the system reached the stationary distribution. The outcome was processed by separating cells into two modes, the ON and the OFF. Cells with the level of protein monomer higher than that of the unstable steady state were considered as ON state, otherwise OFF state. The number of cells in each model was then counted. The "free energy" was calculated with accordance to literature 40 . Specifically, the "free energy" = −log (number of trajectories of (P1 − P2)).

Results
The mutually inhibitory system, shown in Fig. 1, was examined. The reaction network is composed of two genes; each promoter is repressed by the protein produced from the other gene. The trimer of the regulatory protein P1 is the repressor on the promoter of gene 2. Similarly, P2 trimer represses gene 1. Consequently, the ON state of P1 leads to the OFF state of P2 and vice versa. This mutually inhibitory network is one of the canonical patterns of bistability 30 and it performs bimodal distribution in population level as well.
Peptides Y and X are inhibitors to P1 and P2 trimers, respectively. The repression of P2 trimer on promoter 1 is lifted when the gene encoding peptide X is overexpressed. Consequently, the cell stays at the ON state of P1. The conventional beliefs suggest this outcome resulted from merely one factor, the consumption of P2 trimer by inhibitor X. However, in this work, we aim to unveil the role of the stochasticity in cell fate determination. More precisely, we explored how inhibitor exploits stochasticity to drive the cell to a specific state.
The Bimodal Distribution of the System. The bimodal distribution in population level was explored by the stochastic simulations. The stationary distributions were shown in Fig. 2a and b. The ON is in blue color and the OFF in green color; the number above the mode indicates the cell counts in the mode. As one of the features of the mutually inhibitory system, the number of cells at the ON state of P1 should nearly the same as that of the OFF state of P2. Moreover, in this study, the values of parameters for gene 1 is the same as that of gene 2 so the number of cells at the ON state of P1 should theoretically equal to that at ON state of P2. Therefore, for both P1 and P2, the number of cells at ON state is almost the same as that of cells at OFF state with parameters at nominal values.
The overexpression of gene X led to the high level of peptide X. Specifically, the transcription rate of gene X was adjusted to 8 folds of the nominal value and the distributions of P1 and P2 were shown in Fig. 2c and d, respectively. It is clear that the increase of peptide X, as an inhibitor to P2 trimer, pushing cells to the ON state of P1. More than 99% of cells were at P1 ON state or P2 OFF state.

The Change of the Stochasticity by the Increase of Inhibitor Pushing Cells to the Specific
State. In addition to altering deterministic steady states, the inhibitor also caused the change of intracellular stochasticity. By using the term "change of stochasticity", it means all kinds of changes captured by the stochastic model. To investigate how inhibitor exploited the stochasticity to drive cells to the specific state is of interest. It is, of course, not easy in practice to examine the influence of stochastic change on bimodal distribution alone but we can always take advantage of a mathematical model to achieve a better understanding. To this end, we assume that the only consumption of complexes X and Y comes from dissociation; as long as this assumption holds, from the deterministic steady state analysis, the level of free P2 trimer stays the same regardless the level of peptide X (Text S1). Figure 3 is the plot of bistability from steady state analysis of equations in Table S2. The surface in Fig. 3a and b represents the steady states of free P1 and P2 monomer, respectively. The red color indicates the unstable steady state. Along the axis of the transcription rate, which is the expression level of gene X, the bistable behaviors show no shifts. It implies that the deterministic approaches cannot capture the change of the stochasticity. Nevertheless, while applying stochastic simulations, the cells in bimodal distribution migrated to the ON state of P1. Figure 4a and b are the distributions of P1 and P2 with the transcription rate of gene X at nominal value. Remarkably, when the transcription rate of gene X increased to 8 folds of the nominal value ( Fig. 4c and d), cells migrated to the ON state of P1 though there is no difference from the aspect of deterministic steady state analysis. Remarkably, it is also possible to guide nearly all cells to the specific state by intracellular stochasticity alone. As shown in Fig. 4e and f, about 99% of cells were driven to the ON state of P1 with the transcription rate of gene X as 16 folds of the nominal value.
The Contribution of Intracellular Stochasticity of Guiding Cells to the Specific State. In nature, the inhibitor uses two ways to push cells to the specific state. One is through the consumption of P2 trimer and the other is by changing the intracellular stochasticity ( Fig. 4c and d). The former altered the deterministic results but the later was only observed by stochastic simulations. Intriguingly, the latter appears to be influential in redistributing bimodal distribution.
The number of cells at the ON state of P1 for different levels of peptide X was plotted in Fig. 5. The yellow bar indicates the circumstance that the only consumption of complexes X and Y comes from dissociation; the blue bar represents the common situation as described by equations in Table S2. In other words, the yellow bar provides us

Discussion
The inhibitor plays a crucial role in gene regulatory network but the mechanisms of how it affects the bimodal distribution were not well studied. We examined the influence of altering stochasticity on the bimodal distribution and found that it is also the main cause of driving cells to the desired state, as shown in Fig. 5. Such a change of stochasticity is moving along the axis of the transcription rate of gene X in Fig. 3, which shows no difference from the aspect of deterministic steady state analysis. In the other hand, the consumption of regulatory protein is moving along the axis of the degradation rate of P2 monomer, which alters the deterministic bistable behaviors. In a word, inhibitors successfully instructed cells to the specific state at least through these two movements.
The discovery of this study is not limited to the inhibitor. It unveiled the possibility of exploiting the change of stochasticity to control the state of cells in a bimodal distribution. The bimodal distribution raised from bistable switches exists in numerous biological systems 3-5, 14, 17, 18 and how to effectively control the portion of cells Figure 5. The influence of inhibitor on bimodal distribution. The y-axis is the number of cells at the ON state of P1 monomer. We examined the influence of inhibitor on bimodal distribution for two scenarios. For both scenarios, the transcription rate constant of Y is of nominal value and that of X is specified in the x-axis where 0.001 is the nominal value. The yellow bar indicates the circumstance with the assumption that the only consumption of complexes X and Y comes from dissociation. Namely, the yellow bar tells us how the change of stochasticity alone affects the bimodal distribution. In comparison to the blue bar which accounts for both the change of stochasticity and the consumption of P2 trimer, it concludes that the change of the stochasticity is also the main cause of pushing cells to the specific state. in each mode is of high value. Unlike the conventional ways which deal with the bimodal distribution only by the deterministic approaches, our finding, from a new angle, suggested the role of stochasticity is critical. We demonstrated that the altering of stochasticity without consuming regulatory protein is sufficient to push nearly all cells to the specific state (Fig. 4c). It implied that the change of stochasticity, by the inhibitor, can be utilized as a method to redistribute cells. In our previous work 30 we reported that the change of stochasticity, through manipulating the rate of transcription and translation, drove cells to the specific state. However, such an artificial method caused huge noise and compromised the bimodal distribution. The strategy reported in the presented work, from observing the natural interaction of inhibitors and regulatory proteins, does not suffer aforementioned drawbacks and is much easier to be applied. Note that the interaction between proteins and ligands is ubiquitous in biological systems and the proposed method is ready to be applied to all these. Moreover, the change of stochasticity is from the competition of reactions, which is everywhere in regulatory networks. One famous example of practicing the change of stochasticity on redistributing cells to the specific state in bimodal distribution is the novel approach of detecting the latent HIV 41 . The quiescent state is a big problem of detecting HIV infection, and the state switches by the change of the stochasticity offer a better solution. It is always beneficial to know how stochastic fluctuations influences the biological system.
Instead of the protein distribution, "free energy" 40 provides comprehension from a different angle. With the assumption that the only consumption of complexes X and Y comes from dissociation, we observed how the system was altered merely due to the change of stochasticity. Figure 6a is the "free energy" for the case which the transcription rate of gene X is of nominal value and Fig. 6b the case that the transcription rate of gene X is 8 folds of the nominal value. The left valley indicates the "free energy" of P1 OFF and the right valley is that of P1 ON. The ordinate zero of the y-axis is arbitrary but consistent within two sub-figures. The onset of the change of stochasticity by the addition of inhibitor lowered the "free energy" of the right valley (Fig. 6b) and cells incline to stay at the ON state of P1. Moreover, it also raised the "free energy" of the left valley and push cells away from P1 OFF state.
Data availability statement. The model and parameters leading to the findings of this study are available in the manuscript and the supplementary information.