Introduction

Biological systems have evolved by improving the efficiency with which complex regulatory networks control multiple mechanisms in the cell through the fine-tuned balancing of enzymatic reactions. Phosphoinositides are important lipids that are interconverted into each other by multiple enzymatic reactions, which together constitute an example of such a complex network regulating critical cellular functions. Phosphoinositides are key signalling messengers, and several play important parts in regulating physiological processes including vesicular trafficking, transmembrane signalling, ion channel regulation, lipid homeostasis, cytokinesis and organelle identity as characteristic identifiers for different membranes in the cell1,2,3,4,5. It is thus not surprising that phosphoinositides play critical roles in a number of pathological conditions including immunological defence, mediating replication of a number of pathogenic RNA viruses, in the development of the parasite responsible for malaria, in tumorigenesis, Alzheimer’s disease, diabetes, and numerous others6,7,8,9.

The inositol head of phosphoinositides can be phosphorylated at its third, fourth and fifth carbon, thus creating different subspecies. The responsible pathway connects eight metabolites through a dense network of 21 chemical reactions, which are catalysed by 19 kinases and 28 phosphatases10 (Fig. 1). The resulting degree of complexity prevents simple interpretations and renders intuitive predictions of pathway behaviour and regulation unreliable. It is especially difficult to pinpoint the roles of less abundant phosphoinositides, such as phosphatidylinositol 5-phosphate (PI(5)P) and phosphatidylinositol 3,4-biphosphate (PI(3,4)P2). PI(4,5)P2 is present throughout the plasma membrane and considered a general marker for the cell membrane. By contrast, phosphatidylinositol 3,4,5-triphosphate (PI(3,4,5)P3), marks the basolateral part of a polarized cell’s membrane but is absent from the apical part1,11.

Figure 1
figure 1

Map of the phosphoinositide pathway. Red arrows represent fluxes of phosphorylation and blue arrows fluxes of hydrolysis. For each flux, the name (vi→j) and the group of enzymes that catalysed the reaction are shown. Black arrows represent influxes and effluxes of material entering and leaving the system. SIOSS is a group of phosphatases, consisting of SYNJ 1/2, INPP5 B/J/E, OCRL1, SAC2 and SKIP. PI4K + PIP5KI + DVL denotes a complex formed by the three proteins. Proteins separated by commas catalyzed the same reaction. SYNJ: Synaptojanins; INPP5: Inositol polyphosphate 5-phosphatases; OCRL1: Lowe Oculocerebrorenal Syndrome Protein; SAC2: Suppressor of actin ;SKIP: Skeletal muscle and kidney enriched inositol polyphosphate phosphatase.

Full size image

Other phosphoinositides characterize intracellular membranes (Fig. 2). Phosphatidylinositol 3,5-biphosphate (PI(3,5)P2) is typical for multivesicular bodies and lysosomes, whereas PI(4)P is found in the Golgi, and phosphatidylinositol (PI) is located in the endoplasmic reticulum (ER)12. To achieve this distinctive variability in phosphoinositide composition among different membrane compartments, the cell must be able to modulate phosphoinositide metabolism in a targeted, localized manner.

Figure 2
figure 2

Functions of phosphoinositides in the cell. Phosphoinositides are signalling lipids that are cell membrane identifiers. PI(4,5)P2 marks the plasma membrane, PI(3)P the early endosomes, PI(4)P the Golgi, PI(3,5)P2 the late endosomes, PI the ER; finally, PI(3,4,5)P3 is present in the basolateral part of the plasma membrane and absent from the apical part. Phosphoinositides are also second messengers, precursors to other signalling molecules and membrane protein docking sites and regulators.

Full size image

Phosphoinositides also serve as precursors for various signalling molecules, as in the case of PI(4,5)P2, which can be transformed into diacylglycerol (DAG) and inositol triphosphate (IP3) through the action of phospholipase C (PLC). They are furthermore docking sites in the plasma membrane, for instance, for AKT (also known as Protein Kinase B) in the case of PI(3,4,5)P3.

Interestingly, phosphoinositides are also key regulators of ion channel activity11. The epithelial sodium channel (ENaC) is of interest, as it plays a critical role in Cystic Fibrosis (CF), a genetic condition caused by mutations in the gene encoding CFTR, a chloride channel that also regulates other ion conductance, namely through ENaC, across epithelia. In order to keep ENaC open, lysine residues present at the N terminus of the β and δ subunits need to be bound to PI(4,5)P213. This channel function is upregulated in the lungs of individuals with CF, and the increased absorption of sodium and water is considered to be the major cause of lung disease due to critical dehydration of airway surface liquid (ASL)14. The dehydrated ASL and consequent impairment of mucociliary clearance, in turn, is a major cause of respiratory problems in CF13. Thus, a better understanding of the phosphoinositide pathway is of paramount importance, as it may contribute to ameliorating the CF phenotype by manipulating the levels of PI(4,5)P2, which moderate the action of ENaC.

To address the challenge of complexity, it is advantageous to resort to mathematical models, which indeed have already been proposed for particular components of the phosphoinositide pathway. Narang12, Xu15, Nishioka16 and Purvis17 proposed models mainly focused on understanding the dynamics of PI(4,5)P2, PLC, IP3 and DAG, since these molecules are directly associated with calcium release and protein kinase C (PKC) activation, which are important signalling events. Other models, such as those developed by Araia18 and MacNamara19, focus on PI(4,5)P2, PI(3,4,5)P3, phosphoinositide 3-kinase (PI3K), phosphatase and tensin homolog (PTEN) and their roles in cancer. None of these models accounts for all phosphoinositide species. However, the inclusion of less abundant species is important for understanding the distinctions between membrane compartments and for rationalizing the observed impact of several enzyme knock-downs on PI(4,5)P2-mediated ENaC modulation13.

Here, we propose a mathematical model of the complete phosphoinositide pathway. Our primary goal is to shed light on the dynamics of this pathway. Moreover, the model will facilitate a deeper understanding of the unique composition of membranes in different compartments and thereby provide an effective tool for exploring various physiological conditions and their potential treatments, including possible therapeutic targets for CF, cancer and other diseases in which the phosphoinositide pathway plays a critical role.

Results

The prime result of this study is a mathematical model of the phosphoinositide pathway that contains all allegedly relevant molecular components and captures pertinent features of the pathway documented in the literature. The model is certainly not all-encompassing, but detailed enough to serve as a launch pad for future extensions. For instance, it is known that the actual pathway is distributed and compartmentalized. Here, we simulate it restricted to a 1 μm2 patch of plasma membrane, which we consider spatially homogeneous. Nonetheless, the model is designed in a manner that is flexible enough to simulate membrane patches in different compartments, and once the necessary data become available to allow such an extension, it will be easy to block a given reaction a priori if it is known to be absent in that compartment. Alternately, one may perform the same type of parameter optimization as we have done here, but fit experimental observations in different compartments, and this strategy would lead to very low enzyme activities for the corresponding reactions.

The pathway map underlying the model is exhibited in Fig. 1. To facilitate the presentation and discussion of results, each flux is represented by vi→j, and the group of enzymes catalysing it by Ei→j, where the subscripts i and j identify the phosphorylated positions of the substrate and product phosphoinositide species, respectively. The modelled reaction network is based on a review by Balla11, but expanded with information from other sources2,3,10. In particular, we added four fluxes: v0→45, v45→0, v4→34 and v34→4. The first, v0→45, transforms PI into PI(4,5)P2 through a ternary complex of proteins PI4K, PIP5KI and DVL20. This complex is included as one possible molecular complex facilitating the direct channelling of PI into PI(4,5)P2, and it is possible that other protein assemblies could perform this function as well21. v45→0 represents the opposite reaction, which is catalysed by synaptojanins, which are phosphatases that have both a 5-phosphatase domain and a suppressor of actin 1 (SAC1) domain. Balla11 and Hsu22 speculate that the 5-phosphatase domain can transform PI(4,5)P2 into PI(4)P by feeding the SAC1 domain, which dephosphorylates PI(4)P, into PI. Although v0→45 and v45→0 are based on molecular mechanisms that are not generally accepted, their inclusion in the model turned out to be necessary for the maintenance of the PI(4,5)P2 pool when the level of PI(4)P is low. The inclusion of v4→34 has been suggested by Sasaki10, Shewan and Mostov3. Di Paolo and De Camilli2 reported the existence of both v4→34 and v34→4.

Consistency of the Model with Data

As described in the Methods section, model equations were formulated according to Biochemical Systems Theory (BST)23,24. Initial parameter estimates were derived from the literature and from the BRENDA database25. The parameter values were subsequently optimized with a genetic algorithm such that the model matched reported phosphoinositide steady-state levels (Supplementary Table ST1) and dynamic phenomena reported in the literature (Fig. 3a,c and Supplementary Table ST2). This model successfully mimics steady-state levels and 11 out of 13 observed phenomena. The observations that were not replicated are: 1) when the PI levels are reduced, the drop in PI(4,5)P2 levels is not as evident as reported in the literature (Figs 3a and 2) the knockout of myotubularin MTMR2 effects are only partially replicated (Fig. 3c).

Figure 3
figure 3

Perturbations to the phosphoinositide pathway. Blue lines represent experimental observations and bars represent model predictions. (a) Perturbation of PI levels, PI4K and PI5KI activities and resulting effects in PI(4,5)P2 and PI(4)P. γ→0 is decreased to 50% to trigger a decrease of 50% in PI. (b) Perturbation of input fluxes to the levels of PI(4)P and PI(4,5)P2. After stopping all inputs into PI(4)P and PI(4,5)P2, the inputs are re-activated, one at a time, to test if they are sufficient to restore PI(4,5)P2 levels. Enzyme knockouts were simulated by setting the rate constant of the corresponding flux to zero, except for γ0→4, which was decreased to 20% of its original value, in order to avoid numerical errors in the simulation due to very small levels of PI(4)P. (c) Perturbations to MTMR, SYNJ_TMEM55 and PIKfyve that were used to fit the model to the behaviour of phosphoinositides with small pools: PI5P, PI(3,5)P2 and PI(3)P. (d) Consequences of Golgi PI(4)P input (γ→4) for the levels of PI(4)P and PI(4,5)P2 pools. Golgi PI(4)P has a significant impact on the PI(4)P pool but barely affects the PI(4,5)P2 pool. The graphs were created in R42 and the x axis labels were added with PowerPoint.

Full size image

Model sensitivities

The profile of model sensitivities is a double-edged sword. On the one hand, high sensitivities make the model susceptible to unreasonable responses from small perturbations or noise. On the other hand, if the system has a signalling function, small signals must be amplified to have appropriate effects. The model presented here has a stable steady state that is mostly insensitive to parameter changes (Supplementary Table ST8). In fact, the system is robust even to large changes in parameter values (Supplementary Fig. S1). At the same time, the model does exhibit clusters of high sensitivities that are associated with signalling compounds, which one should expect (Fig. 4).

Figure 4
figure 4

High-sensitivity network. Arrows represent amplifying sensitivities with absolute magnitude greater than 1. Red and blue arrows represent positive and negative sensitivities, respectively. The thickness of each arrow is proportional to the magnitude of the corresponding sensitivity.

Full size image

Analysis of low sensitivities and parameter identifiability

Even though most sensitivities are low, one must question how many of the parameters are actually identifiable. To address this question, we performed a Monte Carlo search of the parameter space, which revealed that only 166 out of 79,993 parameter sets tested yield correct steady-state levels (less than 0.15%), given an acceptable material influx into the pathway. All of these 166 solutions have a worst adjustment score than the manually fitted set and the set found with a genetic algorithm (Fig. 5). These results suggest that the model parameterization is sufficiently specific, given the available experimental information.

Figure 5
figure 5

Sum of squared errors for parameter sets detected through Monte-Carlo exploration of the parameter space. All parameter sets shown comply with the following conditions: 1. phosphoinositide steady-state levels are within the intervals retrieved from the literature; 2. the relative amounts between the phosphoinositide pools match the data; 3. influxes are less than 25% of the corresponding phosphoinositide pools; 4. effluxes are less than 7% of the corresponding phosphoinositide pools. The black bar at zero represents the score of a perfect model, the best set found by the genetic algorithm is shown as a blue triangle and the manually found set is the cyan bar. The boxplot concerns the 116 admissible alternative parameter sets. Figure created in R42.

Full size image

High-sensitivity sub-networks

The pairs of model variables and parameters with high sensitivities (Fig. 4) form a network that clusters into four groups around: 1) PI, which is the source of the phosphoinositides; 2) PI(4)P and PI(4,5)P2, which are responsible for plasma membrane identification and PI(4,5)P2 maintenance; 3) the small lipids pools (PI(3)P, PI(5)P and PI(3,5)P2); and 4) PI(3,4,5)P3 and its derivate PI(3,4)P2.

This high-sensitivity network is reflected in a map of parameters that are best poised to serve as “master regulators” for controlling the variables in the different groups. For example, an increase in the levels of PI(3,4,5)P3 and PI(3,4)P2 is most easily accomplished by altering the kinetic order in the flux V45→345. An increase in V345→45 elicits a reduction of PI(3,4,5)P3, which highlights the importance of PI3KI and PTEN for this part of the pathway. If simultaneous increases in the levels of the three phospholipids PI(3)P, PI(5)P and PI(3,5)P2 are required, a researcher should boost V0→3. As an alternative, he could decrease each phospholipid independently manipulating the respective consumption fluxes.

New Insights into the Phosphoinositide System

The model can be used to shed light on the control of the phosphoinositide pathway. Particularly pertinent insights are described in the following subsections.

PI(4,5)P2 is sensitive to PI, PI4K and PIP5KI

Model simulations replicating reported experimental results demonstrate that PI(4,5)P2 is sensitive to the level of PI and to the activities of phosphoinositide 4-kinase (PI4K) and phosphoinositide 4-phosphate 5-kinase (PIP5KI) (Fig. 3a).

PI4K controls PI(4,5)P2 levels: According to the literature, a knockout of phosphoinositide 4-kinase (PI4K) leads to a decrease in PI(4)P and PI(4,5)P2 to 50% of their basal level26. Decreasing PI4K will cause not only the decrease of v0→4 but also V0→45 because this kinase is part of the protein complex that catalyzes V0→45. The model mimics this phenomenon for PI(4,5)P2 although it predicts a more severe drop in the levels of PI(4)P.

PIP5KI controls PI(4,5)P2 levels: One strategy for reducing PI(4,5)P2 levels is to decrease the amount of PIP5KI. This mechanism is probably viable in vivo because a single allele of the PIP5KIγ gene is sufficient to sustain life in mice embryos, whereas knock-out PIP5KIγ mice die shortly after birth27. The same study also showed that α and β genes are not necessary to maintain viability, and their roles are still unclear. Volpicelli-Daley et al.27 furthermore reported that PI(4,5)P2 levels drop around 50% in PIP5KIγ KO mice. Decreasing the activities of PIP5KI (E4→45) and PI4K/PIP5KI (E0→45) to 50% in the model, reduces PI(4,5)P2 to roughly 50% of its basal level.

PI controls PI(4,5)P2 levels: Kim28 reported that a 50% drop in the PI pool causes a similar decrease in PI(4,5)P2 levels. PI(4,5)P2 in the model is sensitive to a reduction in PI but does not drop as much as reported in the literature. Specifically, a 50% drop in PI will only lead to a reduction of 11% in PI(4,5)P2. A 50% drop of PI in the whole cell would also affect other membrane compartments responsible for the production of PI(3)P and PI(4)P. To include this effect, we closed v→4 and v→3. However, this intervention decreases PI(4,5)P2 only to 88% of its basal level. Interestingly, PI(4)P drops to 47%. To achieve a 50% drop in the PI(4,5)P2 pool we would have to shut down v→4 and v→3 completely and reduce the influx of PI,v→0, to 2% of its original value.

PI(3,4,5)P3 levels are sensitive to the concentrations of PTEN and PI3KI

PTEN has been known to be a tumour suppressor for almost twenty years11. This phosphatase hydrolyzes the third position of the phosphoinositide inositol ring in PI(3,4,5)P3 into PI(4,5)P2 and, to a lesser degree, in PI(3,4)P2 into PI(4)P1,10,11. PI3KI phosphorylates the third position of the inositol ring of PI(4,5)P2 into PI(3,4,5)P3, thereby catalysing the inverse reaction of PTEN. This kinase is known to control the cell energetic state and metabolism and thus playing a key role in tumorigenesis11. Bryant and Mostov1 reported that PI(3,4,5)P3 is present at the basolateral membrane, but absent in the apical part, of polarized epithelial cells. PTEN and PI3K are believed to be responsible for this difference. PTEN is present in the apical part and at the tight junctions, where it transforms PI(3,4,5)P3 into PI(4,5)P2. By contrast, PI3K is located in the basolateral part of the membrane and catalyses the opposite reaction from PI(4,5)P2 to PI(3,4,5)P3.

Regulation of PTEN and PI3KI: Cell polarization is highly regulated through mechanisms involving PTEN, PI3K, PI(4,5)P2 and PI(3,4,5)P318,29,30. We investigated to what degree high activity of PTEN (2.3e-15 mg/µm2) and low activity of PI3KI (6.1e-16 mg/µm2) are sufficient to deplete PI(3,4,5)P3 to about 2 molecules/µm2 and thereby mimic the apical membrane configuration. Conversely, we asked if low PTEN (3.9e-17 mg/µm2) and high PI3KI (1.5e-14 mg/µm2) could replicate the basolateral membrane configuration, which is rich in PI(3,4,5)P3 (760 molecules/µm2). Interestingly, model simulations readily mimicked both membrane configurations, which suggests that the model is a satisfactory approximation of the observed phenomena characterizing epithelial and basolateral membrane states (Fig. 6).

Figure 6
figure 6

PI(3,4,5)P3 is sensitive to PTEN when v345→34 is slow. A decrease in PTEN is sufficient to increase the levels of PI(3,4,5)P3 and change the membrane configuration form apical (low PI(3,4,5)P3) to basolateral (high PI(3,4,5)P3). A fast v345→34 will decrease PI(3,4,5)P3 and make the membrane much less sensitive to a PTEN change. A PTEN knockdown of 98.33% does not alter the amount of PI(4,5)P2 in either fast or slow v345→34 conditions. Fast v345→34 increases the levels of PI(3,4)P2 and makes the levels of this lipid dependent on the PI(3,4,5)P3 pool. Slow v345→34 is modelled as γ345→34 = 1e11 and f345→34 = 0.9982. Fast v345→34 is modelled as γ345→34 = 6e13 and f345→34 = 0.9998. The graph was created in R42 and the x axis labels were added with PowerPoint.

Full size image

Flux v345→34 modulates the effects of PTEN: If the flux v345→34 is accelerated to values close to those ones described in the literature for SH2 domain-containing phosphatidylinositol 5-phosphatase (SHIP1), the model predicts a decrease in PI(3,4,5)P3. The surprising consequence of this prediction is that this decrease will lock the membrane in a basal-like configuration and that a knockdown of PTEN will no longer increase PI(3,4,5)P3 (Fig. 6).

Control of PI(4,5)P2 levels

The proposed model is a powerful tool for exploring how the cell controls the phosphoinositide levels in its cell membrane. Due to the multiple functions of PI(4,5)P2, including ion channel activity regulation, cell polarization, and signalling, the control of this phosphoinositide is of particular relevance.

PI(4,5)P2 can be synthesized from three phosphoinositide species in addition to PI (through v0→45), namely PI(4)P, PI(5)P and PI(3,4,5)P3 (Fig. 1). PI(3,4,5)P3 is present in low concentrations and transformed into PI(4,5)P2 mainly by the phosphatase PTEN. The cellular location of this enzyme is tightly regulated, as it is located in non-polarized cells in the cytosol and nucleus most of the time31. PI(5)P also exists as a small pool and its role is not clearly understood. That leaves PI(4)P as the only reasonable candidate for maintaining PI(4,5)P2 levels, besides PI. PI(4)P is a substrate for the kinase PIP5KIγ and has a physiological concentration roughly similar to PI(4,5)P2 pool, i.e., around 10,000 molecules/µm2. However, it has been observed that PI(4,5)P2 levels can be maintained even with low PI(4)P levels11,26. Figure 3b shows the changes in PI(4)P and PI(4,5)P2 levels predicted by the model when different sources are perturbed.

Contribution of v0→45 to PI(4,5)P2 levels: The flux v0→45 represents the direct transformation of PI into PI(4,5)P2 by means of a ternary complex of proteins containing PI4K and PIP5KI20. The model suggests that v0→45 alone can maintain 80% of the basal level of PI(4,5)P2, thereby making it the main source of PI(4,5)P2 (Fig. 3b). This direct transformation of PI into PI(4,5)P2 should exist to ensure the stability of the PI(4,5)P2 pool, and reports in the literature20,32 seem to support this finding.

Contribution of PI4P influx to PI(4,5)P2 levels: The flux v→4 represents the amount of PI(4)P coming from the Golgi through vesicle trafficking or non-vesicle transfer (Fig. 2), which has been reported to constitute a sizeable contribution to the maintenance of plasma membrane PI(4)P, but contributes only moderately to the maintenance of PI(4,5)P232,33. Indeed, the model simulations show that v→4 by itself can maintain PI(4)P at 30% of its basal level and only generates a 9% increase in the PI(4,5)P2 pool (Fig. 3d).

Contribution of PI5P to PI(4,5)P2 levels: The flux v5→45 can maintain the PI(4,5)P2 pool at 34% of its basal level (Fig. 3b). However, the influence of this flux is highly dependent on v→3. If v→3 increases 25 times, which makes this input flux similar to the one for PI(4)P, v5→45 can sustain PI(4,5)P2 levels at 71%. If v→3 increases 50 times, v5→45 can sustain 100% of PI(4,5)P2. This result suggests that PI(5)P may have an influential role in the maintenance of PI(4,5)P2 levels and function as a means of channelling material from PI(3)P toward the linear pathway of PI(4)P, PI(4,5)P2 and PI(3,4,5)P3.

Taken together, these results suggest that the cell employs at least four mechanisms to maintain adequate PI(4,5)P2 levels. This level of redundancy highlights the importance of PI(4,5)P2. Indeed, PI(4,5)P2 is known as a characteristic component of the cell membrane11,26, and it is to be expected that down-regulation of PI(4,5)P2 levels would interfere with the proper functioning of the proteins in the membrane. Compromising these proteins, in turn, would have a negative impact on fundamental processes, such as cellular nutrient intake, information sensing, chemical messaging and the secretion of waste.

Therapeutic Targets for the Modulation of ENaC Activity in CF

The components of the phosphoinositide pathway, and PI(4,5)P2 in particular, are involved in numerous physiological processes, and our model has the potential to deepen our understanding in many of these areas. One specific motivation for us to develop this model was to explore the role of the phosphoinositide pathway in the modulation of the epithelial Na+ channel (ENaC) activity in the lung tissue of patients with CF. ENaC is a sodium and water channel whose activity is upregulated in CF. It is well established that PI(4,5)P2 promotes ENaC activity11,34. We have also previously identified the phosphoinositide pathway to be a key regulator of ENaC13. Indeed, performing an siRNA screen in the CF context using a microscopy-based live-cell assay, we identified 30 enzymes in the phosphoinositide pathway as significant modulators of ENaC activity. We performed independent siRNA knockdowns of phosphoinositide enzymes and re-evaluated ENaC activity with the same live-cell assay. Assuming that if a siRNA increases PI(4,5)P2 it will enhance ENaC activity, we compared ENaC activity results (Supplementary Table ST10) with model predictions of an siRNA effect on PI(4,5)P2.

Our model predictions are consistent with four out of five siRNA assays targeting phosphoinositide kinases. As these assays were not used to calibrate model parameters, this agreement of model predictions with experimental observations supports the validity of our model. Furthermore, model simulations allow us to check if the tested siRNA perturbations may have undesirable side effects on the steady-state profile of the pathway, which were not observable in the original experiments. The results for specific pathway perturbations are discussed in the next sections.

PIP5KI

The most direct and effective way to decrease PI(4,5)P2 levels is by decreasing PIP5KI (E4→45 and E0→45) or enhancing the 5-phosphatases of the SIOSS enzyme group that hydrolyse the fifth position of PI(4,5)P2 (E45→4). It is documented in the literature that decreasing PIP5KI will significantly affect PI(4,5)P2 levels27. The model predicts that a knock-out of PIP5KI will trigger a decrease in PI(4,5)P2 to 13% of its basal steady state (Fig. 7). The performed siRNAs validation tests corroborate the model prediction (Supplementary Table ST10).

Figure 7
figure 7

Predicted changes in PI(4,5)P2 and PI(3,4)P2 levels as a consequence of siRNA knockdown assays. Each kinase is inactivated, one at the time, and phosphatases are upregulated. The protein complex catalysing v0→45 is composed of PIP5KI and PI4K; therefore, when one of these enzymes is knocked out, the complex should be also knocked out. There is also the possibility of knocking out only the complex. Enzymes are coloured according to the classification in the siRNA screens: enzymes activating ENaC are marked red, those inhibiting ENaC are marked green and those exhibiting both effects are marked black. The graph was created in R42 and the x axis labels were added with PowerPoint.

Full size image

An alternative to trigger the decrease in PI(4,5)P2 levels is to increase the activity of 5-phosphatases in the SIOSS enzyme group. Doubling the activity of this phosphatase group in the model results in a 35% decrease in PI(4,5)P2 (Fig. 7). Both our previous dataset13 and results from the new siRNAs validation tests included in the present study (Supplementary Table ST10) are not as conclusive about the SIOSS phosphatases, as for most of the phosphatases tested, which could be a consequence of the unspecific activity that characterizes phosphatases. For example, synaptojanins catalyse several reactions in the pathway and perturbing them would probably cause unexpected side effects.

PI4K

A PI4K knockout affects the fluxes v0→4 and v0→45. It decreases PI(4,5)P2 in 59% (Fig. 7). Accordingly, the model suggests that PI4K should be classified as an ENaC activating gene, which is in line with our previous observations13. An undesirable side effect within these model predictions is the change of the PI(3,4)P2 concentration, a lipid involved in clathrin-coated vesicle formation and activation of AKT. According to the model, this perturbation would cause a 54% decrease in the level of PI(3,4)P2.

PI4K, PIP5KI and DVL protein complex

A knockout of the protein complex formed by PI4K, PIP5KI and segment polarity protein dishevelled homolog (DVL) (PI4K + PIP5KI + DVL, E0→45), which transforms PI directly into PI(4,5)P2, causes a 57% decrease in this lipid (Fig. 7). This perturbation also causes a 7% decrease in the levels of PI(3,4)P2. The PI4K + PIP5KI + DVL protein complex is formed upon wingless-Type MMTV Integration Site Family, Member 3 A (Wnt3a) stimulation. The possibility of targeting the segment polarity protein Dishevelled homolog DVL (DVL) to suppress the formation of the protein complex is interesting because it would avoid interfering with other reactions in the pathway.

PIP5KII and SYNJ/TMEM55

An increase of SYNJ/TMEM55 (E45→5) phosphatases and a decrease of the kinase PIP5KII (E5→45) could decrease PI(4,5)P2 (Fig. 7). The model predicts that SYNJ/TMEM55 has a negligible effect, which is consistent with the literature11 and our previous data13 concerning synaptojanins. However, no phosphatases belonging to the TMEM55 group were screened in the Almaça et al. study. Knocking out PIP5KII reduces the pool of PI(4,5)P2 by 16%. PIP5KII was classified as an ENaC enhancer in both our previous screens13 and in the present siRNA validation tests (Supplementary Table ST10), in agreement with the model prediction. Altering PIP5KII function only causes a 2% decrease in PI(3,4)P2, however perturbing PIP5KII activities could have unforeseen consequences since the role of PI(5)P is not clearly understood and the flux catalyzed by PIP5KII, v5→45, is the main efflux for the PI(5)P pool.

PI3KI and PI3KII

The model predicts that PI(4,5)P2 levels are insensitive to knockouts of PI3KI and PI3KII. We previously classified PI3KI as an ENaC activating gene13, and this is corroborated here by the siRNA validation tests. The PI3KI knockdown increases the level of PI(4,5)P2 if the model parameters are configured to reproduce a basolateral-like membrane composition (enriched in PI(3,4,5)P3). At the same time, this simulated PI3KI knockdown decreases PI(3,4,5)P3 levels, which is also known to control ENaC11. One should note that in polarized cells ENaC localizes to the apical part of the membrane which contains neither PI(3,4,5)P3 nor PI3KI. Therefore, the effect of PI3KI on ENaC may only be observable in non-polarized cells.

The model predicts a negligible influence of PI3KII on PI(4,5)P2 and PI(3,4,5)P3 but causes an almost complete depletion of PI(3,4)P2. We previously13 classified PI3KII as an ENaC inhibiting gene. If this is so, the model suggests that this inhibition could be caused by the depletion of PI(3,4)P2 or components not belonging to the phosphoinositide pathway.

Discussion

In this work, we developed a new mathematical model that captures the complex metabolic network of phosphoinositides. The proposed model successfully replicates the phosphoinositide metabolite levels in mammalian cells and reflects numerous observed phenomena. The model is also able to reproduce the differentiation of the cell membrane into apical and basolateral types.

Using model simulations we were able to dissect the control of the levels of PI(4,5)P2, for which the low abundant phosphoinositide PI(5)P seems to have a significant role as an alternative source. This finding was not detectable in previous models of the pathway, due to their simplifying assumptions.

The results obtained here are of potential interest for a variety of physiological conditions, because the different phosphoinositides play uncounted roles in lipid signalling and membrane dynamics. Of particular interest to us was the fact that the model was helpful in explaining observed effects in a siRNA screen of ENaC modulators in CF. Namely, the model suggested targeting the enzyme that catalyses v4→45, PIP5KI, as the most effective way to decrease the levels of PI(4,5)P2. Targeting PI4K would also reduce PI(4,5)P2 levels significantly, but model simulations point to a possible undesired side effect, namely, the simultaneous reduction of PI(3,4)P2 levels. Targeting the PI4K + PIP5KI + DVL protein complex does not significantly alter other lipids (Supplementary Fig. S2) while yielding a large PI(4,5)P2 reduction. Because this reduction is not as extensive as the one induced by PIP5KI targeting, it may moderate ENaC activity without drastic negative effects in the activity of other proteins regulated by PI(4,5)P2.

The model also suggests that, in order to replicate phenomena retrieved from the literature, v0→45 should be the main flux producing PI(4,5)P2. In particular, this flux may explain the maintenance of PI(4,5)P2 levels when the levels of PI(4)P are low. This result suggests the importance of a close functional relationship between PI4K and PIP5KI. This relationship does not imply that the two kinases must be in physical proximity through this particular protein complex20. They may also work in close proximity within lipid raft-like structures, for example.

The coupling of PI4K and PIP5KI activities may define two configurations of the system. One, where the two kinases are working together closely, in which case they are more sensitive to alterations in PI and in the levels of PI4K. The other configuration is more robust in terms of PI(4,5)P2 levels, where the bulk of this phosphoinositide is created through PI(4)P.

Of course, the model could be improved in the future when new experimental data regarding phosphatases and higher parameter precision are available, but this information is much scarcer than that of kinases. For instance, it is unclear how exactly phosphatases act on the system. Their versatility may suggest the existence of competitive inhibition among their substrates, but this competition could cause substrate coupling, when all substrates of a phosphatase are influenced by the alteration of a single substrate, especially if the phosphatase is saturated35.

Along the same lines, some kinases catalyse multiple reactions. It would thus make sense to consider substrate competition at a more general level. In preliminary studies, we already considered substrate competition, but did not detect significant differences in model behaviours.

Although the model is quite robust, it has few shortcomings. For instance, we estimated values for the parameters γ5→45, f5→45, γ3→35 and γ35→5, which differed somewhat from literature reports, in order to replicate the levels of PI(5)P and PI(3,5)P2. Also, not all phenomena were fully replicated by the model: PI(4,5)P2 did not decrease proportionally to PI, and when MTMR was reduced to 65% (simulating the knockout of MTMR2), PI(5)P did not drop to 20%, only to 98.97%. These discrepancies could be due to gaps in the information about the system. In particular, most of the quantitative data address the total cell and are not membrane specific. Also, in vitro experimental results used to parameterize the model may not truly replicate the system behaviour in vivo.

The current model does not incorporate some regulatory mechanisms which nevertheless may be implemented in future versions. For instance, Bulley et al.36 report the activation of PTEN and PI3KI by their own products, as well as activation of myotubularins and PTEN, and inhibition of SHIP by PI(5)P.

Finally, because the phosphoinositide pathway acts differently in different organelle membranes2, it could be interesting to model not only a cell membrane patch but also the membranes of the Golgi, nucleus and the endoplasmic reticulum with a multiple compartment model featuring lipid transport between them (Fig. 2).

In spite of these simplifications, the proposed model is the first to successfully replicate phosphoinositide metabolism in the mammalian cell membranes. In contrast to earlier models, the model accounts for all known phosphoinositide species and permits unprecedented explorations of the roles of those phosphoinositide’s that are physiologically present in small amounts. The current model, as it is designed, focuses on a single, very small membrane patch of one particular compartment, the plasma membrane, and is not really geared to describe analyses of multiple compartments. One reason is that we simply do not have sufficient metabolic information about fluxes between compartments. If this information were available, we could “multiply” our model several times, eliminate those reactions that are not present in any specific sub-model and use the inter-compartmental fluxes to connect these models.

The model was used to identify the best approaches to control PI(4,5)P2 levels with the goal of establishing new therapeutic targets in the context of CF. The model suggests that the most effective way to accomplish this goal is to decrease the activity of the enzyme PIP5KI (v4→45). Additionally, v0→45 was also found to be very important in the maintenance of PI(4,5)P2 levels. Targeting proteins that are part of the protein complex of PI4K, PIP5KI and DVL or contribute to its control should offer an effective way to control PI(4,5)P2 levels.

In this work, we tried to arrange the current knowledge on the phosphoinositide pathway into a coherent structure. This is an important tool into the understanding of a complex layer of cell regulation that is usually overlooked and can impact fields of study with great potential to improve the human well-being like CF and cancer.

Methods

Model Equations

A dynamical model of phosphoinositide metabolism was designed within the framework of Biochemical Systems Theory (BST)24,37,38,39,40,41, using ordinary differential equations (ODEs) in the format of a generalized mass action (GMA) system. In this approach, each ODE describes the dynamics of a dependent variable X i , which is formulated as a sum of all fluxes that are directly related to this variable; furthermore, each flux vi→j is formulated as a power law function, as indicated in equation (1).

$$\begin{array}{c}\frac{d{X}_{i}}{dt}=\sum _{s}{v}_{s\to i}-\sum _{p}{v}_{i\to p}\\ {v}_{i\to j}={\gamma }_{i\to j}\cdot {E}_{i\to j}\cdot {X}_{i}^{{f}_{i\to j}}\end{array}$$
(1)

The dependent variables (X i ) represent the actual numbers of phosphoinositide molecules (X3: PI(3)P; X4: PI(4)P; X5: PI(5)P; X34: PI(3,4)P2; X35: PI(3,5)P2; X45: PI(4,5)P2; and X345: PI(3,4,5)P3), and of PI (X0: PI) in a membrane patch of size 1 µm2. If more than one substrate contributes to the reaction, or if the reaction is modulated by other variables, the flux term in (1) contains these contributors as additional X’s with their own powers. Because all modelled reactions transform one molecule of some phosphoinositide species into one molecule of another phosphoinositide species, all stoichiometric coefficients are 1 and Eq. (1), therefore, does not explicitly show these coefficients. Kinase catalysed reactions consume ATP and produce ADP, while phosphatase catalysed reactions consume H20 and produce one phosphate ion. These four metabolites were considered to be available in sufficient quantities and not to affect reaction rates.

The model accounts for fluxes transporting PI (v→0), PI(4)P (v→4) and PI(3)P (v→3) into the membrane from the ER, Golgi and endosome, respectively. Additionally, all included species were allowed to be transported out of the membrane via fluxes vi→. We assume that these effluxes follow first-order kinetics (fi→ = 1) and share one common rate constant. The input flux values were restricted in order to allow 4.5% of the membrane phosphoinositides to recycle per minute (see Supplementary Information). Both influxes and effluxes represent transport of lipids that enters or exits the plasma membrane by vesicle- or non-vesicle-mediated transport. The latter can be mediated by specialized proteins like LTP’s or occur spontaneously at membrane contact sites.

Parameter Estimation

Rate constants (γl→j) and kinetic orders (f i→j) were derived from enzyme kinetic parameters obtained in BRENDA or in the literature, as detailed in the Supplementary Information. Enzyme activities (E i→j) and transport fluxes where manually set to approximate reported phosphoinositide steady-state values. This manually adjusted parameter set was used as an initial input for a genetic algorithm (detailed in Supplemental Information). This algorithm found a parameter set that minimized the deviations between model predictions and experimental observations and computing an adjustment score. The parameterized model was characterized through sensitivity and identifiability analysis, and the parameter space was explored with a Monte-Carlo approach (see Supplemental Information).

Model Implementation

The model was implemented in the programming language R v3.1.042 together with the package deSolve43. We used the ODE integration function with the LSODA method. Figures 1, 2 and 4 were created in MS PowerPoint, Fig. 5 and SF1 were created in R42 and finally 3, 6, 7 and SF2 were created in R42 and modified in WS PowerPoint.

Code availability

The R code is available in GITHUB at the following: https://github.com/dolivenca/MK15_phosphoinositide_pathway_model.

siRNA knockdown validity test

To confirm model predictions, selected phosphoinositide pathway hits identified in a large scale siRNA screen13 were validated with an independent round of siRNA knockdown assays. Human alveolar type II epithelial A549 cells (ATCC, Cat no. CCL-185) were transfected with 2 or 3 different siRNAs targeting phosphoinositide pathway enzymes. After transfection the FMP/Amiloride live-cell assay13 was applied to measure ENaC activity. Detailed methods and analysis are described in the Supplementary Information.

Data availability

All data generated or analysed during this study are included in this published article. Please see Supplementary Table ST10 in the Supplementary Information file.