Regulatory interactions maintaining self-renewal of human embryonic stem cells as revealed through a systems analysis of PI3K/AKT pathway

Bioinformatics, Aug 2014

Motivation: Maintenance of the self-renewal state in human embryonic stem cells (hESCs) is the foremost critical step for regenerative therapy applications. The insulin-mediated PI3K/AKT pathway is well appreciated as being the central pathway supporting hESC self-renewal; however, the regulatory interactions in the pathway that maintain cell state are not yet known. Identification of these regulatory pathway components will be critical for designing targeted interventions to facilitate a completely defined platform for hESC propagation and differentiation. Here, we have developed a systems analysis approach to identify regulatory components that control PI3K/AKT pathway in self-renewing hESCs. Results: A detailed mathematical model was adopted to explain the complex regulatory interactions in the PI3K/AKT pathway. We evaluated globally sensitive processes of the pathway in a computationally efficient manner by replacing the detailed model by a surrogate meta-model. Our mathematical analysis, supported by experimental validation, reveals that negative regulators of the molecules IRS1 and PIP3 primarily govern the steady state of the pathway in hESCs. Among the regulators, negative feedback via IRS1 reduces the sensitivity of various reactions associated with direct trunk of the pathway and also constraints the propagation of parameter uncertainty to the levels of post receptor signaling molecules. Furthermore, our results suggest that inhibition of negative feedback can significantly increase p-AKT levels and thereby, better support hESC self-renewal. Our integrated mathematical modeling and experimental workflow demonstrates the significant advantage of computationally efficient meta-model approaches to detect sensitive targets from signaling pathways. Availability and implementation: FORTRAN codes for the PI3K/AKT pathway and the RS-HDMR implementation are available from the authors upon request. Contact: ipb1{at}pitt.edu Supplementary information: Supplementary data are available at Bioinformatics online.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://bioinformatics.oxfordjournals.org/content/30/16/2334.full.pdf

Regulatory interactions maintaining self-renewal of human embryonic stem cells as revealed through a systems analysis of PI3K/AKT pathway

Regulatory interactions maintaining self-renewal of human embryonic stem cells as revealed through a systems analysis of PI3K/AKT pathway Shibin Mathew 2 Sankaramanivel Sundararaj 2 Hikaru Mamiya 1 Ipsita Banerjee 0 1 2 Associate Editor: Janet Kelso 0 McGowan Institute for Regenerative Medicine, University of Pittsburgh , Pittsburgh, PA 15219 , USA 1 Department of Bioengineering 2 Department of Chemical and Petroleum Engineering, University of Pittsburgh , PA 15261 Motivation: Maintenance of the self-renewal state in human embryonic stem cells (hESCs) is the foremost critical step for regenerative therapy applications. The insulin-mediated PI3K/AKT pathway is well appreciated as being the central pathway supporting hESC selfrenewal; however, the regulatory interactions in the pathway that maintain cell state are not yet known. Identification of these regulatory pathway components will be critical for designing targeted interventions to facilitate a completely defined platform for hESC propagation and differentiation. Here, we have developed a systems analysis approach to identify regulatory components that control PI3K/AKT pathway in self-renewing hESCs. Results: A detailed mathematical model was adopted to explain the complex regulatory interactions in the PI3K/AKT pathway. We evaluated globally sensitive processes of the pathway in a computationally efficient manner by replacing the detailed model by a surrogate metamodel. Our mathematical analysis, supported by experimental validation, reveals that negative regulators of the molecules IRS1 and PIP3 primarily govern the steady state of the pathway in hESCs. Among the regulators, negative feedback via IRS1 reduces the sensitivity of various reactions associated with direct trunk of the pathway and also constraints the propagation of parameter uncertainty to the levels of post receptor signaling molecules. Furthermore, our results suggest that inhibition of negative feedback can significantly increase p-AKT levels and thereby, better support hESC self-renewal. Our integrated mathematical modeling and experimental workflow demonstrates the significant advantage of computationally efficient meta-model approaches to detect sensitive targets from signaling pathways. Availability and implementation: FORTRAN codes for the PI3K/AKT pathway and the RS-HDMR implementation are available from the authors upon request. Contact: Supplementary information: Supplementary data are available at Bioinformatics online. 1 INTRODUCTION Long-term maintenance of human embryonic stem cells (hESCs) in the self-renewal state requires a fine balance of many signaling *To whom correspondence should be addressed. pathways, including PI3K, TGF , WNT and ERK (Singh et al., 2012). Several earlier studies have reported that the PI3K/AKT pathway plays a central role in balancing self-renewal and differentiation but with limited mechanistic details (McLean et al., 2007). The pioneering work by Singh et al. first recognized the presence of molecular switches controlled by the PI3K/AKT pathway that promotes self-renewal in its active state and strengthens the differentiation signals in its inactive state (Singh et al., 2012). Kinase p-AKT is a key effector of the PI3K/AKT pathway and participates in critical functions like survival, metabolism, protein synthesis, cell cycle etc. (Taniguchi et al., 2006). In addition, in hESCs, p-AKT regulates the activity of pluripotency factors like c-MYC and controls the levels of differentiation molecules like pSMAD2/3, p-ERK, p-GSK3 in the self-renewal state (Singh et al., 2012). Therefore, maintaining high levels of p-AKT is necessary for long-term self-renewal of hESCs. In spite of the recognized role of p-AKT, there is limited understanding on the maintenance of p-AKT levels in hESCs by the network of regulatory interactions. The PI3K/AKT pathway includes several positive and negative feedback loops and negative regulators like the molecules PTP, PTEN, SHIP that together influence its state (Taniguchi et al., 2006). Analysis of such regulatory interactions will be helpful in the design of targeted molecules to support self-renewal. This, however, requires a quantitative systems-level approach rather than a restricted study of few interactions. Therefore, in the current work, we have used an integrated experimental and computational approach to identify regulatory interactions maintaining p-AKT levels during hESC self-renewal. This being the first effort toward modeling the PI3K/AKT pathway dynamics in hESCs, our workflow includes the following critical steps: (i) determine a mathematical structure adequately describing the hESC system, (ii) determine the parameter range over which the model captures hESC behavior, (iii) identify the relative importance of components of the pathway in hESCs, and (iv) validate the model-predicted sensitive processes in hESCs. These steps result in a validated mathematical representation of the pathway for hESCs. However, any modeling effort of a biological system is incomplete without understanding how parameter variability affects its predictions. We treat this as an important requirement because of the notorious cell-to-cell variability in hESC systems. We, therefore, (v) analyze the propagation of uncertainty in the pathway under various states Fig. 1. PI3K/AKT pathway (A) Insulin receptor level processes. (B) Intracellular signaling in PI3K/AKT pathway. The reactions marked by donut (negative feedback) and star (PTEN and PTP) are perturbed in experiments mentioned in Figure 7 of the identified sensitive regulators. This will, in turn, allow identification of processes that promote a robust behavior under experimental variability. To accomplish these objectives, we started with a wellestablished, ordinary differential equation (ODE) model of insulin-mediated PI3K/AKT pathway developed for adipocytes by Sedaghat et al. (2002). The model is a compendium of accepted knowledge of the pathway, and has been successfully tested in many mammalian systems. We developed a systematic procedure to adopt the Sedaghat model to a system of self-renewing hESCs. We first performed extensive parameter sampling to identify the mechanisms relevant for hESCs. We next evaluated the most significant contributors to the active levels of key molecules using global sensitivity analysis (GSA). We adopted random sampling high-dimensional model representation (RS-HDMR)-based meta-modeling approach to overcome the large Monte Carlo (MC) sampling requirements of the traditional GSA. The model-predicted sensitive processes were successfully validated by a series of perturbation experiments. Our workflow, thus, demonstrates the application of computationally efficient techniques for mechanism detection in uncertain systems like hESCs. SYSTEM AND METHODS Mathematical model of PI3K/AKT signaling The insulin-mediated activation of the PI3K/AKT pathway can be divided into two modules: Module 1: insulin receptor activation, internalization and recycling, and Module 2: post-receptor signaling cascade involving PI3K/AKT (Fig. 1A and B). On stimulation with insulin, insulin receptors are auto-phosphorylated and become available for further signaling. These active receptors then undergo intracellular trafficking as shown in Module 1. The active receptors on the surface propagate the signal to components of PI3K/AKT pathway as shown in Module 2. This includes phosphorylated IRS1 (tyrosine), kinase PI3K and phosphoinositol lipids like PI(3,4,5)P3 (or PIP3 henceforth). The signal then propagates to important kinases like AKT and PKC- . Negative regulators that catalyze de-phosphorylation reactions include the following: PTP1B or PTP (de-phosphoryate active receptors and active IRS1), PTEN (dephosphorylates PIP3 to PI(4,5)P2) and SHIP (de-phosphorylates PIP3 to PI(3,4)P2). The pathway also activates negative feedback loops by serine phosphorylation of IRS1 via kinases like PKC and a double negative loop from AKT resulting in deactivation of PTP. In this article, we relaxed model assumptions by Sedaghat et al. for a more generalized analysis. The details of the ODEs and relaxed assumptions are given in Sections 1 and 2 of the Supplementary Information. The current version of the model comprises 27 reactions, 20 output species and 31 rate parameters. From the rate parameters, 21 were selected as free inputs for GSA (Supplementary Table S1), and the remaining were functions of these selected inputs. Other input parameters included the concentrations of the molecules PTP, PTEN and SHIP. The output molecules of interest were p-IR, p-IRS1 (Y), p-IRS1 (S) and p-AKT. RS-HDMR meta-model for analyzing global sensitivity of high dimensional models Traditional sensitivity analysis techniques are local in nature, and these evaluate the influence of each free parameter in isolation while the remaining parameters are kept constant at their nominal values. This being the first attempt to model hESCs, it was necessary to estimate global sensitivity measures that are applicable in a wide region of the parameter space and capture parameter interactions. The advantages of traditional GSA based on MC methods are, however, challenged by the large number of parameters and the large number of samples required for accurate estimates of the sensitivity indices. To reduce computational cost, we adopted a meta-modeling technique called RS-HDMR developed by Li and Rabitz (2012). 2.2.1 Meta-model development RS-HDMR is a multivariate regression technique that replaces the ODE model by complex surrogate functions (hence meta-model) that describe the input–output (IO) behavior. RS-HDMR decomposes the selected output of the ODE, Y=fðxÞ (say, p-AKT), into component functions represented by the following hierarchical expansion in the k input parameters, x1; x2; x3; :::; xk as: Y=fðxÞ=f0+ fij xi; xj +:::+f123:::kðx1; x2; :::; xkÞ The component functions represent the influence of each type of input combination on the output and are higher-order integrals evaluated as expected values, Eð Þ, at given input instances; Mean : f0=EðYÞ First order : fiðxiÞ=EðYjxiÞ Second order : fijðxi; xjÞ=EðYjxi; xjÞ First-order functions describe the contribution of a single parameter independently of the others, whereas higher-order functions describe the joint contributions of the parameters. Most importantly, the method rests on the fact that the expansion given in Equation (1) often converges rapidly so that higher-order interactions ( 3) are often negligible (Li and Rabitz, 2012). In this article, we used second-order approximation, as it explained most of the variance in the data (480%). The component functions are approximated as mixtures of orthogonal polynomials whose Fig. 2. Workflow adopted for GSA. 105 MC samples of the parameter vector are generated within a specified range of the nominal values using a uniform distribution. The model output generated for each of the samples is subjected to RS-HDMR analysis as described in the main text and Supplementary Information orders and coefficients are evaluated from the IO data using MC integration and least squares regression (See Supplementary Information, Sections 3, 4 and 5 for details). 2.2.2 Computationally efficient GSA The orthogonality of the RS-HDMR component functions allows the estimation of individual variances contributed independently (first order) and jointly (higher order) by the input free parameters. If 2 is the total variance of the model output, YðxÞ, in ½0; 1 k, the decomposition of total variance after assuming that contributions from higher orders ( 3) are negligible is given by the following equation: Z Z = X The Sobol’ indices describing the strength of each combination of parameter can then be estimated using the relations: Details of the evaluation are given in the Supplementary Information, Section 4. The workflow for the entire meta-model development and GSA is represented in Figure 2. Experimental methods Self-renewing H1 hESCs were treated with insulin (100 nM) and the dynamics of components of PI3K/AKT pathway measured for 120 min. The sensitivity of the processes identified by GSA was validated in these cells using three types of perturbations. The molecules p-PKC , PTEN and PTP were perturbed from their basal levels in self-renewal media, and resulting effect on key molecules like p-AKT was measured. See Supplementary Information, Section 8, for details of the experimental protocol, reagents and normalization method. 3.1 Selection of parameter ranges for describing PI3K/AKT pathway in hESCs The first step in this exercise is to determine whether the structure of Sedaghat model is sufficient to address hESC dynamics. In the Fig. 3. Experimental dynamics of pathway molecules in hESCs. (A) p-IGF1R and p-AKT (B) p-IR, p-IRS1 (Y) and (S). The phosphorylated protein measurements are normalized to corresponding control (see Supplementary Information, Section 8) and further by time 0 original Sedaghat model (Sedaghat et al., 2002), several phosphorylation and de-phosphorylation rate constants of Module 2 were fixed by equilibrium ratios. For extension to hESCs, we decoupled these rates by varying them independently in a realistic range. These rate constants include—k7/IRp, k-7, k8, k-8, k9stim and k-9. To determine their realistic ranges, we first measured the dynamics of key molecules of the pathway under insulin stimulation in H1 hESCs. Representative molecules from different positions of the signaling pathway were selected for analysis: early (p-IR, p-IGF1R), mid (p-IRS1 (Y), p-AKT) and late (p-IRS1 (S)). Figure 3 presents the dynamics of the measured phospho-proteins. p-AKT showed an overshoot behavior and it settled at intermediate values of 3-folds by 120 min (Fig. 3A). The levels of p-IR and p-IGF1R increased rapidly to 2- and 8-fold within 15 minutes and remained at these levels till the end of stimulation (Fig. 3A and B). Among the IRS1 molecules, tyrosine phosphorylation (p-IRS1 (Y)) showed rapid increase with maximum levels reached by 15 min and then a downregulation to intermediate values, which remained constant upto 120 min (Fig. 3B). Serine phosphorylation of IRS1 (p-IRS1 (S)) showed an initial dip followed by upregulation at 30 min and then stabilization to basal levels. A distinct negative correlation was observed between p-IRS1 (Y) and p-IRS1 (S) dynamics. To select parameter ranges, we first explored the behavior of the relaxed Sedaghat model over a broad parameter range (100fold around nominal values). The basal levels of PTP, SHIP and PTEN and initial concentration of insulin were also included as additional parameters, resulting in 25 free parameters. 105 random samples were drawn from a uniform distribution of 25 parameters, and the ODE model was integrated for each sample till 120 min. This being the first modeling effort of the pathway for hESCs, the actual parameter distribution is not known a priori. Hence, a uniform distribution was used here. The profiles were normalized to the maximum level and clustered using k-means algorithm on p-AKT dynamics. There were three major clusters (Fig. 4A–C and Supplementary Table S2 and Fig. 4. K-means clusters for p-AKT. (A–C) Three key dynamic clusters observed in the parameter space. Cluster 3 showed overshoot behavior and intermediate steady state levels. Cluster centroid is shown by red curve and the shaded region shows the cluster extent. Blue lines indicate the overlaid experimental data. Hierarchical clustering analysis gave same type of major clusters (See Supplementary Figure S3 and Supplementary Information Section 7). (D) Parameter k-7 contained in clusters, C1, C2, C3 from (A-C). The black bar indicates the final selected range for sensitivity analysis and the red arrow shows the location of the nominal value Supplementary Information Section 6). Among these, clusters 2 and 3 showed the characteristic overshoot behavior, but only cluster 3 showed intermediate steady state values as seen in the experiments in Figure 3A. We also saw good correlation between experimental and clustered profiles for other molecules, p-IR and p-IRS1 (Y) (Supplementary Fig. S1). Cluster 3 being our primary cluster of interest; we wanted to select parameter ranges defining this cluster. Hence, we first analyzed parameters primarily segregating the three clusters. It was observed that of the 25 parameters, only six were varying between the three clusters: k7/IRp, k-7, k8, k-8, k9stim and k-9. Among these, the activation parameters were high and deactivation parameters were comparatively low in cluster 3. The distribution for de-activation parameter, k-7 for each cluster is shown in Figure 4D. The distributions for remaining parameters are shown in Supplementary Figure S2. Based on these distributions, we narrowed the range of the six parameters around the peaks for cluster 3 (Fig. 4D and Supplementary Fig. S2). The chosen restricted ranges of the parameters are presented in Supplementary Table S1. We next performed GSA in the restricted parameter ranges to identify the key parameters regulating the self-renewal state of hESCs. Meta-model representation and efficient GSA To reduce the computational cost associated with GSA, we adopted RS-HDMR (Li and Rabitz, 2012), to explore the model IO behavior and also to rank the parameters using Sobol’ sensitivity indices. The performance characteristics of the method has been validated for diverse physical systems (Li and Rabitz, 2012). Here, we have applied the technique for a high-dimensional, non-linear signal transduction model. To check the validity of RS-HDMR and its computational efficiency, we compared the Sobol’ indices estimated using RS-HDMR with direct MC evaluations for one of the output molecule, p-IR. First, we analyzed the effect of sample size on the sensitivity index evaluated both by direct MC and by RS-HDMR. Direct MC identified the basal rate of receptor recycling, k-4, as the most sensitive parameter with a Sobol’ index (Si) of 0.38. Figure 5A compares the convergence of Sobol’ indices for this parameter by the two methods, for sample size (N) ranging from 102 to 105. We see that the two methods converge to the same value by a sample size of 104. Relative ranking of parameters are often more informative than the actual value of the sensitivity index. Hence, Figure 5B represents the ranking of the parameters obtained by direct MC analysis for 105 samples and compares it with RS-HDMR predictions at different sample sizes. Overall, it is observed that the parameters with higher sensitivity were predicted with great accuracy even at low sampling of 103. Beyond the fourth-ranked parameter for 103 samples and eight-ranked parameter for 104, there is considerable deviation from MC analysis. The higher sampling size of 105, however, closely predicts the direct MC estimates for most of the parameters throughout the range. Hence, a sample size of 105 was chosen for the remaining work. The primary motivation for adopting RS-HDMR is the computational efficiency of the algorithm. The computational cost of traditional GSA is primarily associated with the need for repeated sampling to evaluate the integrals. In contrast, using polynomial approximations in RS-HDMR, both low- and high-order Sobol’ indices can be estimated simultaneously from one set of MC samples. As an example, we have tabulated the typical time requirements for first-order indices in Figure 5C. Obtaining all the first-order Sobol’ indices by RS-HDMR requires 9 min, although it takes 300 min for direct MC. In RS-HDMR, a critical source of error is the MC integration approximation for the high dimensional integrals. The error of this approximation is inversely proportional to the sample size as N1/2 and favorably independent of the dimension (Li and Rabitz, 2012). Sensitive parameters for key molecules of the pathway On confirming the accuracy of second-order RS-HDMR for the current model, we next determined the globally sensitive parameters for the steady state of four molecules of the pathway: p-IR, p-IRS1 (Y), p-IRS1 (S) and p-AKT (Fig. 6A). The resulting output distributions are presented in Supplementary Figure S4 and overall performance in Supplementary Table S3. We see that the sensitivity contribution of each parameter to the different outputs is different. In Module 1, rates of recycling of the nonphosphorylated receptors, k-4 and active receptor internalization, k4’ were important. These parameters primarily changed the steady state levels of the active receptors p-IR and also affected p-IRS1 (Y) and p-AKT levels to a small extent. The initial insulin levels and the binding rate of insulin to the receptors, k1, were the other parameters affecting p-IR, but they did not affect the other molecules. The important parameters from Module 2 were primarily associated with negative regulators of the pathway. These included many of the deactivation rates: de-phosphorylation of p-IRS1 (Y), k-7, deactivation of PI3K, k-8, and de-phosphorylation of PIP3 to PI(4,5)P2, k-9. All these parameters were upstream of PIP3 and affected the p-AKT levels considerably. p-IRS1 (Y) was affected mostly by k-7 whereas p-IRS1 (S) was affected mostly by k-8 and k-9. Another important set of parameters was associated with negative feedback by p-PKC and subsequent serine phosphorylation of IRS1. These included the Hill equation parameters, Vmax, Kd, n and the IRS1 serine phosphorylation rates, k7’ and k-7’. These parameters significantly affect the p-IRS1 (S) levels followed by p-IRS1 (Y) and p-AKT. Thus, it is seen that for the intracellular molecules, most of the sensitive parameters are ‘negative regulators’ of molecules upstream of PIP3 or they are associated with ‘negative feedback’. Additionally, for p-AKT, there were also important contributions from the second-order indices. For p-AKT, the second-order indices contributed to 28% of the variance. Figure 6B presents the heat-map of the scaled second-order Sobol’ indices. The sensitive parameters were found to cluster in two groups: Group A: interaction between negative regulators upstream of PIP3 ( 7%) and Group B: interaction between negative feedback parameters with negative regulators from Group A ( 6%) (Supplementary Table S4). It is important to note that these parameters also had important first-order contributions. Second-order interactions further increase the sensitivity of these parameters (Supplementary Fig. S5). Experimental validation of sensitive processes in hESCs As seen from the previous section, the model-predicted sensitive perturbations on the system include (i) negative feedback via p-PKC and (ii) negative regulators that affect the de-phosphorylation of PIP3 (Summarized in Fig. 1). To validate the sensitivity of these processes in hESCs, we performed targeted perturbations (See Supplementary Information Section 8), and we carefully chose perturbations that would result in increase in p-AKT levels to avoid differentiation. For negative feedback, the ideal candidate for perturbation is p-PKC and for de-phosphorylation rates upstream of PIP3, two such candidates exist, namely PTEN and PTP. 3.4.1 Influence of negative feedback on self-renewing hESCs Self-renewing hESCs were subjected to 20 min treatment of 10 M GO€ 6983 (PKC- inhibitor) based on Zheng et al. (2000). In our experiments, we see large change in p-AKT (6-fold) and moderate change in p-IRS1 (Y) (1.5-fold) levels for small decrease in p-PKC level (Fig. 7A). Similar trends were predicted by model simulations (Fig 7A, red lines). Thus, small changes in the strength of negative feedback propagated to large changes in the levels of p-AKT. While sensitivity indices indicate the importance of a parameter on a specific output, the directionality of the effect, positive or negative, cannot be directly deduced from Sobol’ indices. The meta-model representation of RS-HDMR is particularly suited for such deduction as this information is contained in the hierarchical component functions of the decomposition. Increasing the strength of negative feedback (increasing k7’) decreases p-AKT and p-IRS1 (Y) (Supplementary Fig. S6). Therefore, we see a positive correlation between p-AKT and p-IRS1 (Y). The sensitivity based on firstorder indices show smaller increase in p-AKT, but the experiments clearly show a large increase in p-AKT levels. We envisage this to be the effect of non-linear influence of the other sensitive processes. For example, model simulations show that under complete inhibition of negative feedback, the levels of p-AKT will rise considerably if the influence of negative regulators like PTEN are low to begin with (low k-9) or the levels of positive regulators like PIP3 are high (Supplementary Fig. S7). 3.4.2 Influence of PIP3 de-phosphorylation on self-renewing hESCs under basal PTP levels Using a direct inhibitor of PTEN, bpV(HOpic), we studied its effects on p-AKT and p-IRS1 (Y) (Schmid et al., 2004). At a low concentration of 100 nM for 24 h, the compound is known to suppress active PTEN, thereby increasing inactive p-PTEN. Our experimental data shows that small change in PTEN resulted in 2-fold increase in p-AKT (Fig. 7A). The levels of p–IRS1 (Y) showed a 1.5-fold decrease. Similar trends were predicted by model simulations (Fig 7A, red lines) and by RS-HDMR (Fig. 7B). This effect is primarily because of the strengthening of negative feedback leading to indirect inhibition of p-IRS1 (Y). This also points to the fact that p-AKT is more sensitive to PIP3 levels as compared with p-IRS1 (Y). 3.4.3 Influence of PIP3 de-phosphorylation on self-renewing hESCs under PTP inhibition Next, we validated the effect of combined PTEN and PTP inhibition to check whether PTP inhibition increases p-IRS1 (Y) when PTEN is still inhibited. At higher concentrations, the same complex bpV(HOpic) can inhibit both PTEN and PTP. We treated hESCs with 1 M inhibitor for 20 min following Schmid et al. (2004). Our experimental results show a proportional increase in all the three molecules, p-AKT, p-IRS1 (Y) and p-PTEN (Fig. 7A). Similar trends were predicted by model simulations (Fig. 7A, red lines) and RS-HDMR (Fig. 7B). The increase in p-IRS1 (Y) is primarily because of PTP inhibition (decrease in k-7), as PTEN inhibition alone resulted in a decrease in p-IRS1 (Y). Thus, under PTEN + PTP inhibition, p-IRS1 (Y) overcomes the effect of increase in downstream negative feedback. In our experimental studies, we observed that the largest increase in p-AKT levels was brought about by the inhibition of negative feedback. This was also predicted by the second-order RS-HDMR component functions for p-AKT. The second-order component functions for one such combination, k-9 and k7’ is presented as a heat-map in Figure 7C for p-AKT and p-IRS1 (Y). For low k-9 and low k7’, there was significant positive contribution to the p-AKT levels but not for p-IRS1 (Y). Combining this with the first-order contribution from k-9 and k7’ in this regime and together with the mean, f0, we get the total p-AKT of 47% [f0 + f (k-9) + f (k7’) + f (k-9, k7’)]. Alternatively, this effect is reduced when the negative feedback is strengthened. For example, in the low k-9 and high k7’ regime, the same contribution to p-AKT amounts to 23%. Therefore, strong negative feedback can considerably decrease the sensitivity of other reactions involving PTEN and PTP. Robustness of system behavior under parameter uncertainty While mathematically representing hESC systems, it is important to consider the effect of variability as observed in different experimental repeats. To test how parameter uncertainty influences the variability in the model output, we chose a biologically realistic log-normal distribution of the parameters centered around the nominal values. From the negative feedback parameters, the value of the most sensitive parameter, k7’ was varied as follows: (i) no negative feedback, k7’ = 0, (ii) nominal level of feedback, k7’ = 0.347, and (iii) 10 times stronger feedback, k7’ = 3.47. Figure 8A presents the probability distribution of steady state levels of p-AKT and p-IRS1 (S), respectively, when exposed to parametric uncertainty, under these different levels of negative feedback. The model shows that strengthening the feedback parameter reduces variance in the distribution for each molecule. Hence, the robustness of the system to input perturbations is enhanced in the presence of a strong negative feedback. In addition, it was observed that the distribution for p-AKT was narrower as compared with p-IRS1 (S). To verify this observation, we plotted steady state levels of p-IR, p-IRS1 (S) and p-AKT from five different experimental repeats in Figure 8B. We see a comparatively high variability in p-IR and p-IRS1 (S) levels but interestingly, p-AKT shows a narrow range of variability (see cell-to-cell distribution in Supplementary Fig. S8). Negative feedback, thus, plays an important role in maintaining robust levels of the important molecules under experimental variability. Additionally, Pearson pairwise correlation between the molecules showed a good agreement between the experiments and model predictions (Fig. 8C and D). Meta-model approximation for GSA of complex biological signal transduction models In this article, we, for the first time, present a detailed analysis of the regulatory interactions in PI3K/AKT pathway actively maintaining the self-renewal state of hESCs. A key step in our workflow is the analysis of the uncertainty associated with the strength of interactions in the PI3K/AKT pathway of hESCs using a meta-model-based GSA. GSA captures the complete non-linear associations between the model parameters in a sufficiently wide region of the parameter space and is suited for nonlinear systems. The traditional methods of GSA are variance decomposition schemes involving exhaustive exploration of the parameter space. This renders the use of a detailed parametric analysis of large-scale ODEs expensive restricting the modeler to fairly simple local analysis. To explore IO relationships efficiently, we adopted meta-model approach called RS-HDMR to obtain accurate information on the sensitive model parameters. RS-HDMR constructs a complex surrogate function to replace the ODE model and evaluates MC integrals of the Sobol’ indices efficiently. The method has been proven to perform well in a variety of engineering systems where parameter uncertainty is a norm. In the current work, the method was explored for a signal transduction model in a 25 dimensional parameter space. We demonstrated that the method is especially accurate in identifying the most sensitive parameters and their functional relationships to the output. Also, the technique, being independent of the total number of parameters can be applied to larger models common in systems biology. Hence, it is a promising alternative to evaluate global sensitivities with computational efficiency, instead of settling for locally based approximations as commonly done in signal transduction studies. Implications of GSA of molecules belonging to PI3K/AKT pathway Through our systems level analysis, we found that parameters associated with post receptor processes can affect the levels of intracellular molecules of the PI3K/AKT pathway more than the receptor level processes for high insulin concentrations. The high sensitivity of the post receptor processes is in support with previous experimental and modeling analyses that show that the functionality of insulin signaling can be severely affected by mutations associated with post receptor signaling molecules (Nyman et al., 2012). Additionally, on removal of equilibrium relationship between the forward and backward reactions, it was observed that the de-phosphorylation reactions of the direct cascade were highly sensitive although the phosphorylation reactions were comparatively insensitive as seen in other systems like the MAPK/ERK pathway (Yoon and Deisboeck, 2009). This is an important relation, as many of the de-phosphorylation reactions are dependent on the concentrations or functionality of phosphatases that can vary widely with the cell type and state and are also implicated widely in diseased states (Yoon et al., 2010). Our analysis also shows that because of the competing nature of the reactions in this module of the pathway, there is a considerable non-linear outcome to simultaneous changes in the sensitive parameters, for example, reduction of sensitivity to PTEN and PTP inhibition by strengthening negative feedback. Such non-linear behavior was seen in an experimental study in hESCs, where inhibition of a direct phosphatase to p-AKT rendered the p-AKT levels insensitive to PI3K inhibition (Yoon et al., 2010). Our current analysis thus highlights the importance of non-linear interactions in determining the effect of perturbations on the pathway components. This is an important outcome of mechanistic modeling and will prove useful in the design of targeted interventions for many systems. Modeling self-renewal in hESCs 4.3.1 Processes affecting pathway dynamics in hESCs Insulin stimulation experiments in self-renewing hESCs showed an overshoot behavior in the dynamics of post receptor molecules. Two possible candidates have been identified to explain such overshoot behavior in other cell types: (i) receptor internalization and (ii) downstream negative feedback from still unknown regulators of receptor de-phosphorylation (Nyman et al., 2012). Usually, there is combined contribution from both the processes. In our hESC system, however, we do not see substantial overshoot behavior in p-IR dynamics, but surely there is a clear decrease in downstream p-IRS1 (Y) levels and an accompanying increase in p-IRS1 (S) levels. This indicates that the negative feedback acting at the level of IRS1 is responsible for decrease in p-IRS1 (Y). This was also seen from the negative correlation between p-IRS1 (Y) and p-IRS1 (S). Our clustering analysis showed that many of the de-phosphorylation reactions above PIP3 had to be maintained at low levels, and it was necessary to couple this with an existing negative feedback to see an overshoot behavior. 4.3.2 Processes affecting p-AKT levels in self-renewing hESCs The central molecule like p-AKT can counterbalance the mechanisms that may lead to differentiation and support mechanisms that can lead to self-renewal (Singh et al., 2012). Any increase in p-AKT levels has been shown to result in increased stability and self-renewal capacity of hESC cultures. For example, the levels of the active form of self-renewal molecule c-MYC can increase with increase in p-AKT levels (Yoon et al., 2010). Yet, there are limited efforts to understand how regulatory mechanisms affect long-term maintenance of selfrenewal in hESCs, which was the focus of the current study. Under the current culture conditions, the receptor level processes were found to be less sensitive. Therefore, a promising strategy to increase p-AKT levels is inhibition of internal signals that suppress p-AKT. Our results suggest that inhibition of negative feedback via PKC- is one such mechanism. A parallel experimental study has recently demonstrated the positive attribute of PKC inhibition in hESC self-renewal, but did not offer any mechanistic insight (Gafni et al., 2013). Additionally, model analysis shows that any perturbation in the phosphorylation and de-phosphorylation reactions of this pathway (for example, PTEN and PTP inhibition) would still need the removal of negative feedback mechanisms to increase sensitivity to these interventions. From uncertainty propagation analysis, we show that negative feedback also increases the robustness of p-AKT levels to variations in the levels of upstream molecules. Mechanisms like negative feedback are known to impart robustness in many biological systems. Interestingly, the steady state correlation between molecules of the pathway held under experimental variability. In conclusion, the strength of negative feedback needs to be maintained in a fine balance. Weakening the negative feedback is favorable for self-renewal but is associated with increased variability. The current work has developed a mathematical structure of the PI3K/AKT pathway, validated by experiments, to describe self-renewing hESCs. Adoption of RS-HDMR, a powerful metamodeling technique, allowed feasible evaluation of GSA of the complex non-linear pathway. An important conclusion from our study is that the maintenance of p-AKT levels, and hence the self-renewal state of hESCs, is controlled by many of the negative processes of the pathway. Additionally, non-linear interactions identified by RS-HDMR show that the existing negative feedback plays an important role of desensitizing the pathway to input perturbations and thus, regulates the steady state distribution of molecules in self-renewing hESCs. Inhibition of negative feedback can significantly increase p-AKT levels and support self-renewal, but with a trade-off associated with increased variability. Such mechanistic analysis of new systems like hESCs is a critical step toward identification of new targets for optimizing cell culture conditions. The authors thank the Stem Cell Core Facility at the University of Pittsburgh for supply of H1 hESCs. Funding: This work was supported by the National Institute of Health New Innovator Award [1DP2OD006491-01 to I.B.]. Conflict of Interest: none declared. Gafni , O. et al. ( 2013 ) Derivation of novel human ground state na€ıve pluripotent stem cells . Nature , 504 , 282 - 286 . Li , G.Y. and Rabitz , H. ( 2012 ) General formulation of HDMR component functions with independent and correlated variables . J. Math. Chem., 50 , 99 - 130 . McLean , A.B. et al. ( 2007 ) Activin A efficiently specifies definitive endoderm from human embryonic stem cells only when phosphatidylinositol 3 kinase signaling is suppressed . Stem Cells , 25 , 29 - 38 . Nyman , E. et al. ( 2012 ) Insulin signaling - mathematical modeling comes of age . Trends Endocrinol. Metab. , 23 , 107 - 115 . Schmid , A.C. et al. ( 2004 ) Bisperoxovanadium compounds are potent PTEN inhibitors . FEBS Lett. , 566 , 35 - 38 . Sedaghat , A.R . et al. ( 2002 ) A mathematical model of metabolic insulin signaling pathways . Am. J. Physiol. Endocrinol. Metab. , 283 , E1084 - E1101 . Singh , A.M. et al. ( 2012 ) Signaling network crosstalk in human pluripotent cells: a Smad2/3-regulated switch that controls the balance between self-renewal and differentiation . Cell Stem Cell , 10 , 312 - 326 . Taniguchi , C.M. et al. ( 2006 ) Critical nodes in signalling pathways: insights into insulin action . Nat. Rev. Mol. Cell Biol ., 7 , 85 - 96 . Yoon , J. and Deisboeck , T.S. ( 2009 ) Investigating differential dynamics of the MAPK signaling cascade using a multi-parametric global sensitivity analysis . PLoS One , 4 , e4560 . Yoon , B.S. et al. ( 2010 ) Optimal suppression of protein phosphatase 2A activity is critical for maintenance of human embryonic stem cell self-renewal . Stem Cells , 28 , 874 - 884 . Zheng , W.H. et al. ( 2000 ) Stimulation of protein kinase C modulates insulin-like growth factor-1-induced akt activation in PC12 cells . J. Biol. Chem. , 275 , 13377 - 13385 .


This is a preview of a remote PDF: https://bioinformatics.oxfordjournals.org/content/30/16/2334.full.pdf

Shibin Mathew, Sankaramanivel Sundararaj, Hikaru Mamiya, Ipsita Banerjee. Regulatory interactions maintaining self-renewal of human embryonic stem cells as revealed through a systems analysis of PI3K/AKT pathway, Bioinformatics, 2014, 2334-2342, DOI: 10.1093/bioinformatics/btu209