1.
Introduction
Pulsatile dynamics, which are characterized by recurrent patterns in molecular abundance or activity, have been demonstrated in diverse biological signaling pathways [1,2,3]. Dynamics per se can encode information from the stimulus and elicit complex downstream response kinetics, as demonstrated in p53 signaling pathway [4]. Rich dynamics, such as damped or sustained oscillations heavily depend on interconnected feedback loops in signaling network [5]. Therefore, identifying the role of feedback loops in p53 signaling may help decode the dynamic responses and provide clues to cell fate decision.
The transcription factor p53 represents a central node for majority of tumor associated pathways and serves as a critical target for pharmacological intervention [6]. The extraordinary diversity by which p53 responds to cellular stress leads to the tumor suppressor phenotype of p53 protein [7]. Under normal conditions, p53 protein is rapidly destabilized by an E3 ubiquitin ligase MDM2 [8]. MDM2 gene is also transcriptionally activated by p53, thereby forming an autoregulatory circuit [9]. The dynamics of p53 may display stereotypic pulses with concise amplitude and duration, extended large pulses or monotonic induction depending on the type and strength of stimuli [10,11]. Recently, Chen et al. have identified a novel negative feedback in p53 signaling pathway during hepatocellular carcinoma (HCC) progression [12]. The glycogen synthase 2 (Gys2) can restrict tumor development in HBV-related HCC by competitively disrupting MDM2-p53 binding and preventing p53 from MDM2 mediated degradation [12]. Furthermore, p53 binds to the Gys2 promoter to inhibit the transcription of Gys2 and the transcriptional repression is dramatically facilitated by Hepatitis B virus (HBV) X protein (HBx) [12]. Gys2 functions as a key enzyme in glycogen synthesis and glycogen accumulation contributes largely to tumor growth suppression in p53-dependent manner [13]. Therefore, unraveling the role of novel p53-Gys2 negative feedback in p53 signaling pathway from theoretical perspective may enrich the clinical significance of Gys2 during HCC progression.
In current work, we constructed a minimal model incorporating p53-MDM2 and p53-Gys2 negative feedback loops. The results showed that p53 displayed complex dynamics ranging from sustained or damped pulses, maintenance at low levels to monotonic induction. Using multi-scale modeling, we recapitulated the key experimental findings both ‘in vitro’ and ‘in vivo’. Overexpressing Gys2 could inhibit tumor progression whereas reducing Gys2 expression facilitated tumor growth in silico. HBx introduction dramatically increased the number of tumor cells and diminished the pool of Gys2. We further showed that optional p53-Gys2 negative feedback provided robustness against oscillatory and impulse disturbances. Instead, p53-MDM2 negative feedback can more effectively reshape amplitude and duration of p53 pulses. Therefore, dual negative feedback loops in p53 signaling can provide both robustness and tunability, whereas strong perturbations to p53-Gys2 feedback may destroy cellular homeostasis leading to tumor progression. Our model has uncovered the important role of p53-Gys2 negative feedback loop during tumor progression and may provide a theoretical framework for homeostatic control in p53 signaling pathway.
2.
Materials and method
2.1. Model construction
We constructed a minimal three-component p53 model which incorporated dual negative feedback loops (p53-MDM2 and p53-Gys2). For details in model construction, equations and parameters, please refer to supplementary materials.
2.2. Multi-scale modeling of p53 signaling network
The agent-based model is composed of four scales to simulate the temporal variation of cells (i.e. agents): (1) molecular; (2) cellular; (3) microenvironment and (4) tissue scale. The molecular scale describes the p53 signaling pathway. The cellular scale characterizes the phenotypic switching of cells. The tumor microenvironment establishes a link between angiogenesis and phenotypic switching. The p53 protein not only affects the migration of cells but also regulates cell division [14,15]. Meanwhile, cells may secret vascular endothelial growth factor (VEGF) to induce angiogenesis (tissue scale). Formulation of multi-scale model was described in supplementary materials.
2.3. Local sensitivity coefficient
A dynamic system can be defined by x’ = F(x, θ), where x and θ denote state and parameter vector, respectively. We used the integrated p53 dynamics or maximal p53 response as the metrics to calculate local sensitivities. Local sensitivity coefficient Sω is defined as
2.4. Disturbance
The disturbance to the minimal system was described by adding a term Ip53 to Eq.1 of the ordinary differential equations (ODEs) as previously described [16]. We explored two types of disturbance: an oscillating disturbance and an impulse disturbance. See supplementary methods for details.
2.5. Bifurcation
The bifurcation is obtained using Oscill8 (https://sourceforge.net/projects/oscill8/). Raw data were exported and bifurcation diagram was plotted using MATLAB.
2.6. Model simulation
The ordinary differential equations (ODEs) were integrated with the ode23s solver embedded in MATLAB. Codes to simulate in vivo and in vitro cell growth were drafted using MATLAB (MathWork, Natick, MA, R2018b).
3.
Results and discussion
3.1. A minimal p53 model shows complex dynamics
To evaluate the role of p53-Gys2 negative feedback, we constructed a minimal model with interconnected feedback loops (Figure 1A and supplemental methods). Notably, the p53-Gys2 and p53-MDM2 negative feedback converged in MDM2 regulation (Figure 1A). The deterministic p53 system revealed sustained pulses with a roughly 5-hour period (Figure 1B). We noted that MDM2 mediated p53 degradation (Vp) and p53 induced MDM2 expression (Vy and K2) dramatically affected p53 dynamics whereas p53 dynamics was less affected by p53-Gys2 loop (rGys, Kp53, dGys and KGys, Figure 1C). Since sensitivity coefficient might be only valid locally, we also adopted a global sensitivity measure termed partial rank correlation coefficient (PRCC). We found that local and global sensitivity showed significantly correlation implying that the local sensitivity coefficients were also applicable (Figure S1). We then introduced random fluctuations in kinetic parameters to investigate cell-to-cell variations in p53 dynamics (see supplementary methods) [17]. Totally, 5000 stochastic runs were performed, and the results showed that p53 exhibited complex dynamics from damped pulses, sustained pulses, maintenance at low levels to monotonic elevation (Figure 1D and Figure S2). Nearly half stochastic runs reserved pulsatile behaviors with characteristic periods (Figures 1E and 1F). Meanwhile, p53 could also be maintained at low levels or become highly activated (Figure 1F). These results suggested that p53 signaling may display complex dynamics when confronted with noises.
3.2. HBx promotes HCC progression
HCC primarily originates from chronic hepatitis B virus (HBV) infection and HBx serves the key oncogenic protein encoded by HBV [18]. HBx may facilitate p53-mediated Gys2 suppression [12]. Therefore, we used HBx expression to model HCC for simplicity by decreasing Kp53 to 0.01 (normally, Kp53 was set to be 0.1, Table S1). In deterministic simulation, we found that p53 and Gys2 expressions were both suppressed (Figures 2A and S3A). We then performed stochastic simulations by sampling all parameters from lognormal distributions (5000 sets). Results showed that p53 also displayed complex dynamics with pulsing or non-pulsing behaviors (Figure 2B). HBx expression significantly decreased Gys2 expression (Figures 1D, 2B, blue curves and Figure S3B), which is qualitatively consistent with experimental findings [12]. To more faithfully evaluate the role of HBx, we used multi-scale modeling. In addition to the molecular scale, a multi-scale model also evaluates cell contact, spatial restriction, microenvironmental effect and angiogenesis, simultaneously [19]. We mimicked the ‘colony formation assay in vitro’ by randomly assigning 30 ‘seeds’ in the lattice, in which each ‘seed’ contained random numbers of cells (27~46). The results showed that HBx expression significantly accelerated HCC cell proliferation (Figures 2C and 2D, two representative simulations were shown). HBx expression markedly increased the number of active cells (i.e. migrating and dividing cells, Figure 2D) and decreased the levels of Gys2 (Figure 2E). Not surprisingly, p53 levels were also substantially downregulated with HBx expression at 30 days (Figure S3C). The results suggested that HBx may promote HCC development by inhibiting p53 and Gys2 expression, which are qualitatively consistent with experimental data [12].
3.3. Regulating Gys2 expression affects cell proliferation
We next investigated whether alterations in Gys2 abundance might influence cell proliferation. Gys2 in cells (i.e agents in silico) was either knockdown (KD, rGys = 0.0002) or overexpressed (OV, rGys = 0.0008) to mimic the experimental conditions. Compared with cell proliferation in wild type (WT) cells, Gys2 KD resulted in explosive increase of active cells whereas Gys2 OV markedly inhibited cell growth (Figures 3A and 3B). As expected, cells with Gys2 KD showed significantly reduced Gys2 expression whereas Gys2 OV elevated Gys2 expression (Figure 3C). The pseudocolor plot additionally revealed altered Gys2 protein levels in WT, Gys2 KD and Gys2 OV cells (Figure 3D). These results suggested that Gys2 played an inhibitory role for cell proliferation and qualitatively matched the experimental findings [12].
3.4. Modeling the ‘in vivo’ effect of HBx and Gys2
We then simulated ‘in vivo’ vascular tumor growth. Angiogenesis, diffusion of microenvironmental factors such as glucose and VEGF, as well as chemotaxis were all incorporated. The results showed that HBx expression indeed facilitated vascular tumor growth compared with those under WT condition (Figure 4A). The numbers of active cells and endothelial cells were significantly higher when HBx was expressed (Figures 4B and 4C). Tumors tended to growth towards vessels where glucose concentrations were higher (Figures 4A and S4). Meanwhile, the branching microvasculature was denser and more complex near the solid tumors where the VEGF was highly accumulated (Figures 4A and S4). When Gys2 was overexpressed, however, the vascular tumor growth was markedly suppressed (Figures 4A-4C). Consistently, Gys2 silence profoundly accelerated tumor volume expansion (Figures 4A-4C). HBx indeed significantly abolished Gys2 expression (Figure 4D). Gys2 OV or Gys2 KD markedly increased or decreased Gys2 levels, respectively (Figure 4D). Glucose concentration was actively consumed at tumor locations where VEGF secretion was also relatively enhanced (Figures 4E and 4F). These data suggested that HBx and Gys2 regulation could alter Gys2 abundance and remodel tumor growth.
3.5. Optimal p53-Gys2 negative feedback provides robustness against external disturbance
We next investigated the role of p53-Gys2 feedback loop. For simplicity of exposition, we introduced a ‘disturbance’ as previous described [11,16]. Oscillatory and impulse disturbance were both considered to explore the resilient properties of p53 dynamics. Combinatorial variations in disturbance period (or 1/frequency) and Kp53 (indicative of negative feedback strength [20]) were applied together with different levels of disturbance strength (10%, 20% and 50%, see supplemental methods), the deviations were then evaluated (Figure 5A). We found that at relatively low disturbance strength (10%), the deviation to initial states was minimal unless the disturbance frequency approached the natural p53 frequency if Kp53 was located within an ‘optimal region’ (i.e. two blue sink regions in the pseudocolor plots, Figures 5A, left and 5B). We further noted that the significant deviation from the initial p53 pulsatile dynamics was largely ascribed to a phase difference without remarkable change in the pulsing period (Figure 5B, middle, the light purple dot in 5A). As the disturbance strength was increased (20%), the deviations became more significant especially for longer disturbance periods (e.g. > 5 h, Figure 5A, middle). When the disturbance strength is relatively large (50%), the sink regions below the natural p53 period has shrunk whereas the region above has vanished (Figures 5A, right and 5C, middle and right). Interestingly, we found that p53 system could effectively filter out high-frequency (low period) disturbances irrespective of the disturbance strength provided that Kp53
was within a certain range (Figures 5A, S5A and S5B). Notably, p53-Wip1-ATM and p53-MDM2 negative feedback loops are both required to produce sustained p53 pulses theoretically and experimentally in DNA damage response [1] although recently there is still an excellent p53 model with only p53-MDM2 circuit [21]. Our model does not incorporate p53-Wip1-ATM negative feedback, which is primarily functional in DNA damage response [22]. However, Wip1 can dephosphorylate Ser395-phosphorylated MDM2 to enhance MDM2 mediated p53 degradation [23]. As a result, MDM2 dephosphorylation by Wip1 may form an additional negative feedback loop in p53 network (Figure S5C). The p53-Wip1 and p53-Gys2 negative feedback loops showed similarities and converge in MDM2 regulation (Figure S5C). We anticipated that the p53-Wip1/p53-Gys2 negative feedback loops may complement or reinforce each other to provide heterogeneous redundancy in robustness [24]. These analyses suggested that p53 signaling was robust against moderate levels of disturbances with optimal p53-Gys2 negative feedback loop.
3.6. Resilient properties of p53-Gys2 negative feedback against impulse disturbance
We then focused on the transient behaviors of p53 dynamics with impulse disturbances. The impulse disturbance was reshaped with different amplitudes and durations (Figure 6A). In default settings, p53 performed recurrent pulses, which could be characterized as closed circles similar to limit cycles in phase plane (Figure 6B). After an impulse disturbance, cellular dynamics may deviate from the initial phase portraits and then reverted to the closed circles. We measured the time needed to return to the limit cycle (‘resilient time’) after the impulse disturbance and quantified its deviation to the resilient time without disturbance (WT, the red dots represented the start at which the system re-entered the limit cycles, Figure 6B). We found that at Kp53 = 0.03, at which the p53-Gys2 feedback strength was outside the ‘optimal region’ (Figure 5A), the deviations were relatively large (roughly 2~8 h, Figure 6C). However, for Kp53 = 0.3 or Kp53 = 0.05, which was located well within the optimal region or close to the boundary respectively, the p53 system after impulse disturbances reverted to the limit cycles similar to that under WT condition (Figures 6B, 6D and S6B). We then explored the bifurcation properties of p53 system and found that the amplitudes and periods of p53 pulses were relatively stable within a certain range of Kp53
values (Figures S6C and S6D). These results suggested that optimal strength of p53-Gys2 negative feedback may improve rejection of impulse disturbance and provide robustness.
3.7. The p53-MDM2 provides tunability to pulsatile dynamics
Coupled negative feedback circuits may conduct diverse regulatory functions [25]. We then explored whether p53-MDM2 loop showed different features. We focused on the pulsatile characteristics of p53 such as duration and amplitude with varying strength of dual feedback loops. Deterministic results showed that the existence of sustained pulse, duration and amplitude were sensitive to changes in p53-MDM2 feedback loop (Figures 7A and 7B). P53 trajectories were less heterogenous across different stochastic runs (Figure 7C, left). Further randomization in Kp53 values (see supplementary materials) did not lead to substantial variations in p53 dynamic features (Figures 7C, middle). However, remarkable stochasticity could be observed if Vy was randomized (Figure 7C, right). The distribution of pulse duration and amplitude indicated that stochasticity in Vy resulted in significantly more heterogeneous p53 dynamics (Figures 7D and 7E). Bifurcation diagram further confirmed that changing Vy could strongly affect pulse amplitude and duration (Figure S7A and S7B). These data suggested that p53-MDM2 negative feedback loop could effectively regulate pulsatile p53 dynamics compared with p53-Gys2 circuit.
3.8. Balance in p53-Gys2 negative feedback maintains cellular homeostasis
We then evaluated the effect of perturbations in p53-Gys2 circuit and associated effects on cellular homeostasis. Previous study has shown that chronic hepatitis B virus infection results in cycles of inflammation in liver and serves as a risk factor for HCC [26]. Expression of HBx significantly lowered the glycogen content to promote HCC development and HBx facilitates p53 mediated transcriptional Gys2 suppression [12]. Meanwhile, cells have evolved multiple mechanisms against spontaneous apoptosis to maintain hepatocyte homeostasis [27]. Consistently, we found that depleting p53-Gys2 negative feedback loop (e.g. Kp53
= +∞) in our model may lead to significant apoptosis (Figures S8A-S8C) and therefore impede normal cell regeneration (Figure S8D). On the contrary, HBx may utilize a different strategy by strengthening p53-mediated Gys2 suppression leading to tumorigenesis (Figure S8D). Consequently, either depletion (Δp53-Gys2) or dramatic enhancement (HBx expression) in p53-Gys2 negative feedback loop might both result in breakdown of hepatocellular homeostasis and unfavorable outcomes. In addition, hepatocytes may replenish themselves gradually in homeostasis and actively against liver injury [28]. Remarkable liver regeneration is a prominent feature against acute injury and is important for drug detoxification [29]. However, regeneration is highly related to tumorigenesis, in which hyperactivated regeneration will possibly lead to tumorigenic phenotypes [30]. P53 can act as a coordinator to block the transition from regeneration to tumors, and p53 levels are fluctuated during liver regeneration [29]. In priming phase, p53 is inhibited to allow cell cycle re-entry whereas p53 becomes elevated to protect genome stability at late stage [29]. We argued that p53 signaling may transiently exit the homeostatic domain to prime cell division and re-enter the optimal region to assure fidelity of regeneration (Figure S8D). The latter assumption should be verified in future. Taken together, mild perturbation in p53-Gys2 circuit might facilitate liver regeneration and homeostasis whereas strong perturbations may lead to tumorigenesis or unfavorable apoptosis.
4.
Conclusion
In current work, we investigated the role of a novel p53-Gys2 negative feedback circuit using multi-scale modeling. The interconnected dual negative feedback loops in p53 signaling pathway may provide both robustness and tunability for cellular homeostasis against tumor progression. The performance of feedback controls in p53 network might be extended to other signaling pathways with similar architectures.
Acknowledgments
This work is supported by National Natural Science Foundation of China (31971185, 31800316) and Quality Engineering Project of Anhui College Education (2018jyssf086).
Conflict of interest
All authors declare no conflicts of interest in this paper.
DownLoad:
