Introduction

Mutations in the cystic fibrosis transmembrane conductance regulator (CFTR), a chloride and bicarbonate membrane channel, can cause problems in several organs and lead to cystic fibrosis (CF). In the lungs, the production of thick, dehydrated mucus associated with these mutations leads to recurrent infections and frequent inflammation events that eventually compromise organ function. Life expectancy for subjects with CF has improved considerably1, and promising new drugs like lumacaftor and ivacaftor that correct the mutated protein were recently brought to the market2. Despite these advancements, a complete cure for CF has not yet been achieved, in part due to mutations that are not treatable with the available drugs.

While CF is directly caused by mutations in the CFTR protein, other ion channels, like the Epithelium Sodium Channel (ENaC), are affected as well. It has been hypothesized that the lack of CFTR in CF lungs causes ENaC function to increase. Consequently, large amounts of sodium and water are absorbed, implying that this channel may contribute to the thick dehydrated mucus that characterizes the disease.

The reasons for ENaC’s upregulation in CF are not clear, but there is no shortage of hypotheses3,4,5,6,7,8,9. Among these, Tarran and colleagues10,11,12 advanced the idea that short palate lung and nasal epithelial clone 1 (SPLUNC1) is deactivated in the absence of CFTR, either by increased acidity in the specific environment of the CF lung or due to an increased presence of proteases13. The inactive SPLUNC1 does not fulfill its role of inducing the channel disassembly and removal of ENaC α and γ subunits from the plasma membrane9. As a consequence, ENaC channels have a lower probability of being removed from the plasma membrane and, ultimately, ENaC function is up-regulated. We adopted this view because it was at least partially validated by results obtained under physiologically and clinically relevant conditions14.

Phosphoinositides are relatively rare membrane lipids with various signaling functions. They are characteristic for different types of cell membranes and play a dynamic role in cell process control as second messengers and precursors of other messengers. As a result, phosphoinositides are important in a myriad of cell functions like cytoskeleton formation, chemotaxis, cell polarization, T cell activation and cytokinesis15,16. Here, we are particularly interested in their function as docking sites for proteins to the cell membrane and as membrane protein regulators.

Several studies have shown that two of these lipids, phosphatidylinositol 4,5-bisphosphate (PI(4,5)P2) and phosphatidylinositol 3,4,5-triphosphate (PI(3,4,5)P3), have an effect on ENaC17,18,19. Their key precursor, phosphatidylinositol (PI), is created in the ER from phosphatidic acid (PA) and transported to the plasma membrane, where it is phosphorylated into other phosphoinositide species (Fig. 1, upper box). PI(4,5)P2 is cleaved by phospholipase C (PLC) into inositol triphosphate (IP3) and diacylglycerol (DAG) and transformed into PA, which is transported back to the ER, closing a functional cycle.

Figure 1
figure 1

Diagram of the integrated model. Upper box: Phosphoinositide model: Many components were already represented in our previous model24, while the components associated with PA and DAG (left) are new extensions. The arrow from PI(4,5)P2 to ENaC represents the control that the lipid exerts over the channel and constitutes the link between the two component models. Thick black arrows represent influxes and effluxes of material entering and leaving the system. Red and blue arrows represent fluxes of phosphorylation and hydrolysis, respectively. Thin black arrows represent regulations. For each flux, the name (vi→j) and the group of enzymes that catalyze the reaction are shown. Orange arrows represent phospholipase fluxes. PTEN and PI3KI have an active (a) and inactive (i) state. O_I_SK_SA2 is a group of phosphatases, consisting of OCRL1, INPP5 B/J, SAC2 and SKIP. PI4K + PIP5KI + DVL denotes a complex formed by the three proteins. Proteins separated by commas catalyze the same reaction. Lower box: ENaC—ASL model diagram. Green arrows represent influx and efflux of ENaC channels. Cyan arrows represent influxes and effluxes of material for ASL. Black arrows depict regulation. The value of PI(4,5)P2 regulates the open probability of ENaC. For equations and parameter values please see our previous paper29. Abbreviations: V1—ENaC influx. V2—ASL independent ENaC efflux. V3—ASL dependent ENaC efflux. V4—CFTR independent ASL influx. V5—CFTR dependent ASL influx. V6—ENaC dependent ASL efflux. V7—ENaC independent ASL efflux.

Kota et al.20 offered an explanation linking phosphoinositides to ENaC control. They found that when the intracellular N-termini of ENaC connect to phosphoinositides, a conformational change occurs that exposes the extracellular loops of the channel to proteases. Severing these loops leads to an increase in the probability Po that ENaC is open. Thus, if the levels of PI(4,5)P2 or PI(3,4,5)P3 are decreased, proteases do not cut ENaC’s extracellular loops as often, and the channel activity is reduced.

Diacylglycerol kinase (DGK) transforms DAG into PA. Almaça et al.21 found that inhibiting DGK causes a moderation of ENaC activity and normalizes the increased sodium channel activity in CF. The authors hypothesized that inhibiting DGK might bring the recycling of phosphoinositides to a halt, which in turn would decrease the levels of PI(4,5)P2 and PI(3,4,5)P3 and cause the observed ENaC moderation. However, DGK is active in the plasma membrane, while phosphoinositide synthesis occurs in the ER. Transport of lipids between membranes of different compartments is typically mediated by vesicles or specialized proteins, but quantitative details about the dynamics of this transport in the CF lung, if it occurs, are lacking. Thus, many open questions remain unanswered. Crucially, the dynamics of the control of ENaC by DGK is still unclear, and it remains to be explored whether and how DGK could potentially be used as a therapy in situations where ENaC function is dysregulated.

In this work, we investigate ENaC control by DGK. As an alternative to Almaça’s hypothesis, we proffer that ENaC is regulated by PI(4,5)P2, which is produced by type-I phosphatidylinositol-4-phosphate 5-kinase (PIP5KI) under the control of PA, which in turn is produced from DAG under the control of DGK (Fig. 2).

Figure 2
figure 2

Consequences of inhibiting DGK for phosphoinositide metabolism, ENaC and ASL, visualized with sizes of boxes and arrows. Panel a: Normal state. Panel b: DGK inhibition reduces the production of PA which, in turn, reduces the production of PI(4,5)P2, which is catalyzed by PIP5KI. Low levels of PI(4,5)P2 reduce ENaC activity and the absorption of ASL stimulated by this channel, which consequently leads to an increase in ASL thickness.

We furthermore suggest that PI(4,5)P2 specifically influences the probability Po that ENaC is open (“open-probability”) rather than the number N of channels in the membrane, in agreement with Ma et al.22. N itself is influenced by ubiquitination by NEDD4-2, the Neural precursor cell Expressed Developmentally Down-regulated protein 4–223 and an interaction with SPLUNC1, the protein Short Palate Lung and Nasal epithelial Clone 1, which is abundant in the airways9.

Our core objective in this work is to test these hypotheses. Our strategy of computational modeling is specifically designed to yield a deeper understanding of how DGK and phosphoinositides control ENaC activity. The proposed model consists of two modules, adapted from our previous work, that are embedded in an appropriate context and explain the regulation of ENaC function. One module addresses the phosphoinositide pathway, while the other captures the regulation of ENaC and the airway surface liquid (ASL) (Fig. 1). The merging of these modules allows, for the first time, a detailed assessment of the dynamics of ENaC regulation by phosphoinositides. Our results show that the combined model matches data from several sources remarkably well. In particular, testing our model against the observations of Almaça et al. yields good agreement. This agreement suggests that the hypothesis of ENaC regulation by DGK, accomplished though PA activation of PIP5KI, is valid and can be used to explore therapeutic interventions in clinical conditions where ENaC is thought to be relevant.

Results

Model expansion, parameterization and validation

The phosphoinositide pathway module was adapted from our recent work24 and expanded with processes that are fundamental to the regulation of the pathway. These expansions allow us to account for: (1) competition among enzymes for the same substrates25; (2) regulation among phosphatase and the tensin homolog (PTEN), phosphoinositide 3-kinase (PI3KI), PI(4,5)P2 and PI(3,4,5)P3; (3) the cleavage of PI(4,5)P2 into DAG and IP3 by phospholipase C (PLC); (4) the production of PA by the phosphorylation of DAG by DGK; (5) hydrolysis of PA back to DAG by phosphatidate phosphatases (LLP’s); (6) replenishment of the PA pool from phosphatidylcholine (PC) by phospholipase D (PLD); and (7) activation of PIP5KI, the enzyme that transforms phosphatidylinositol 4-phosphate (PI(4)P) into PI(4,5)P2, by PA, as previously described in the literature26,27,28.

The map for the expanded phosphoinositide pathway module is depicted in the upper box of Fig. 1. Fluxes and equations are presented in Table 1 and parameters and initial values are given in Table 2. References supporting the parameter values of the model can be found in Table S2 and in the publications underlying the original models24,29.

Table 1 Details of the integrated model.
Table 2 Model parameters and initial values.

The earlier phosphoinositide pathway model was successfully tested against a long list of phenomena reported in the literature24. Corresponding results for the extended phosphoinositide model are summarized in Fig. 3. As is to be expected, the new modules and extensions slightly affect the model fits in comparison to the previous sub-model, but the combined phosphoinositide-ENaC model yields fits to the data that are similar. In addition, the combined model generates genuinely new results, which are summarized in Table S1 and detailed in the following, along with reports from the literature.

Figure 3
figure 3

Perturbations to the phosphoinositide pathway. Blue lines represent experimental observations and bars represent model predictions. (a) Perturbation of PI levels and PI4K and PI5KI activities and resulting effects on PI(4,5)P2 and PI(4)P. γ→0 is reduced to 50% to trigger a decrease of 50% in PI. (b) Perturbation of input fluxes into the pools 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 fluxes 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 behavior of phosphoinositides with small pools: PI5P, phosphatidylinositol 3,5-bisphosphate (PI(3,5)P2) and phosphatidylinositol 3-phosphate (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. γa→b denotes the rate constant of the flux that transforms a into b. γ→0 denotes an influx of material into the system, in this case into the PI pool. γ4→ denotes an efflux of material out of the system, in this case from the PI(4)P pool. Some examples: γ→0 denotes the rate constant of the influx of material to the PI pool from the exterior of the system, γ0→4 denotes the rate constant of the flux that represents the transformation of PI into PI(4)P, γ5→45 denotes the rate constant of the flux that represents the transformation of PI(5)P into PI(4,5)P2.

Figure 3a, c show model results as bar plots, while the blue lines represent data. When phosphatidylinositol (PI) is depleted, the decreases in PI(4)P and PI(4,5)P2 are now slightly closer to 50%, whereas the results of perturbations in PI4K and PI4P5KI are somewhat inferior (Fig. 3a). Perturbations that affect the lipids to a smaller degree produce similar results (Fig. 3c).

Figure 3b indicates that we can still create conditions where PI(4,5)P2 levels are maintained with low levels of PI(4)P, but the PI(4,5)P2 pool is much more dependent on the V0→45 flux. In the new model, the contribution of phosphatidylinositol 5-phosphate (PI(5)P) to the PI(4,5)P2 pool is small. Figure 3c, d show that PI(4)P is also more dependent on the V0→4 flux.

Several other datasets were retrieved from the literature and used to fit and validate the extended phosphoinositide pathway module. These results are described in the Supplements (Sections 1 and 3), where also further details about the phosphoinositide model expansion can be found.

The second module captures the dynamics of ENaC and ASL regulation by SPLUNC1 and PI(4,5)P2 (Fig. 1, lower box); it was fully described in our previous work29. This module simulates the number of ENaCs, which are regulated by SPLUNC1 and ASL thickness. ASL thickness, in turn, is regulated by CFTR and ENaC activities. The latter is the product of the number of ENaCs and their collective open probability (Po), which is controlled by PI(4,5)P2 levels. The ENaC—ASL model simulates both healthy and CF lungs through the change of CFTR activity.

The phosphoinositide pathway and ENaC—ASL modules were constructed to simulate the dynamics of regulation in a 1 µm thick patch of cell membrane in human bronchial epithelial cells. Fluxes and equations are presented in Table 1, and parameters and initial values are given in Table 2.

The key connection between the two modules is PI(4,5)P2. The phosphoinositide pathway module determines the PI(4,5)P2 levels, while the ENaC—ASL module translates these levels into ENaC activity. Therefore, the combined model permits novel explorations of the interactions between phosphoinositides and the dynamics of ENaC and ASL (Fig. 1).

The combined model with PIP5KI activation by PA replicates the effects of a DGK knockdown

One hypothesis to be tested with the proposed model is that the DGK knockdown effect on ENaC activity is mediated by PI(4,5)P2 through PA control of PIP5KI. The work of Antonescu et al.30 demonstrated that inhibition of DGK causes a decrease of about 50% in the level of PA. If our hypothesis regarding the role of DGK is valid, the reduction in PA should decrease PI(4,5)P2 production by PIP5KI. Indeed, the model exhibits a decrease of 52% in PA when DGK is inhibited 75%. The same perturbation reduces PI(4,5)P2 by 28%.

Almaça et al.21 performed siRNA screens to identify modulators of ENaC activity and observed that DGK inhibition reduces ENaC activity to WT basal levels in CF F508del cells, but has no significant effect in WT cells. In CF, when PLC is inhibited, DGK inhibition has no effect on ENaC. When PLC is activated, ENaC function drops to WT basal levels, and inhibiting DGK causes an additional small decrease in ENaC function.

A comparison between model results and data is illustrated in Fig. 4. Given the good qualitative (and semi-quantitative) agreement between the data and the results from our model analyses, we can cautiously conclude that our hypothesis regarding the mechanisms of ENaC regulation provides a good explanation for the activity of ENaC in WT and CF. Our combined model tests the effect of the activation of PIP5KI by PA, but, in this implementation, cannot test the PI recycling hypothesis, as it focuses on a plasma membrane patch. Further details about a pertinent model simulation can be found in the Supplements (Section 2.1).

Figure 4
figure 4

Model results demonstrating consequences of DGK inhibition on ENaC activity in healthy subjects (WT) and patients with CF, CF with PLC inhibition, CF with PLC activation, and CF with PI3KI inhibition of ENaC activity. Bars and confidence intervals represent data, as reported by Almaça et al.21. The dark grey bars correspond to basal levels of DGK activity, the light grey bars correspond to 25% inhibition of DGK (DGK−). Blue lines represent model results. (1st to 4th bars) DGK inhibition in CF (CF/DGK−) reduces ENaC’s activity to a similar level as the channel activity in WT, as described by Almaça et al.21. (5th and 6th bars) With PLC inhibition, a decrease in PI(4,5)P2 caused by DGK inhibition does not affect ENaC’s action, as proposed by Almaça et al.21. The discrepancy between the data and the model results are probably due to saturation of the florescence reporter used. (7th and 8th bars) When PLC is activated, ENaC activity is reduced. DGK inhibition further reduces the channel activity. The data provided by Almaça et al. suggest a decrease after DGK inhibition, which however is not significant. *Data for WT and WT/DGK− were collected using transepithelial voltage measurements; otherwise, a voltage dependent FMP assay was used (cf. Almaça et al.21 for further details).

Some discrepancies persist when PLC is perturbed. In the case of PLC inhibition (CF\PLC-), Almaça’s data indicate the same magnitude of ENaC activity as in the unperturbed CF state (CF), while the model predicts a much higher value. This discrepancy is presumably due to the fact that Almaça et al. used a fluorescence essay where the marker was saturated for values around 100% of ENaC activity. In the case of PLC activation (CF\PLC+), the model replicates the case without DGK perturbation but when the kinase is inhibited, a slightly more pronounced reduction in ENaC activity is predicted by the model.

The combined model with phosphoinositide recycling, but without PIP5KI activation by PA, does not replicate the effects of a DGK knockdown

Phosphoinositides are synthesized in the ER from PA and transported to the plasma membrane where they are phosphorylated, cleaved by PLC into IP3 and DAG and transformed back to PA, which is transported back to the ER, thereby closing a functional cycle. Almaça et al.21 hypothesized that inhibiting DGK would bring the recycling of phosphoinositides to a halt, which in turn would decrease PI(4,5)P2 and consequently reduce ENaC action.

To test the hypothesis of Almaça and colleagues, we temporarily altered our model by removing the regulation of PIP5KI by PA and allowing the efflux from the PA pool, VPA→, to supply the PI pool with material, which simulates the transformation of PA into PI in the ER. When DGK is inhibited in this configuration, it no longer influences ENaC activity (Figure S7a). This result happens because the efflux of PA is only 5% of the influx of PI. Furthermore, the influx of PI, which is independent of the plasma membrane PA efflux (V→0) that represents PA created de novo in the ER, is enough to sustain the PI pool when DGK is inhibited.

ENaC can be made sensitive to DGK if the much larger PI influx is made exactly proportional to the PA efflux, which however is unrealistic (details in the Supplements (Section 2.2) and Figure S7c).

Suratekar’s phosphoinositide cycle model replicates DGK knockdown effects only if PIP5KI is activated by PA

Suratekar et al.31 recently modeled the phosphoinositide cycle in the plasma membrane and in the ER. Coupling our ENaC—ASL module with their model allowed us to test Almaça’s hypothesis in a different, almost independent manner.

We used the version of the Suratekar model that the authors showed to be consistent with all data available. This version uses Michaelis–Menten kinetics to represent an open cycle with influx of PA into the ER and efflux out of DAG into the plasma membrane; it does not include regulation of PIP5KI by PA. This representation turned out to be inconsistent with Almaça’s observations of ENaC moderation under DGK inhibition (Fig. 5a). Suratekar’s model contains an influx of PA into the ER that is independent of plasma membrane PA. When we inhibit DGK, which severely reduces PA in the plasma membrane, this influx suffices to maintain the necessary phosphoinositide levels.

Figure 5
figure 5

Suratekar’s model of phosphoinositide recycling by itself does not replicate Almaça’s observations of ENaC moderation when DGK is inhibited. The results presented here were obtained with a model combining the Suratekar model with our ENaC—ASL model. Panels (a) and (b) represent ENaC activity under healthy (WT) and CF conditions. (a) Data from the Suratekar model fail to replicate experimental observations. (b) If Suratekar’s model of phosphoinositide recycling is modified to include PA regulation of PIP5KI and an efflux from the PI(4)P pool, the resulting model does replicate Suratekar’s data and Almaça’s observations of DGK control of ENaC. DGK inhibition causes a moderation of ENaC’s action. In CF, this moderation brings the ENaC activity close to the WT ENaC activity. (c) Plots represent simulated lipid ratios in Drosophila melanogaster photoreceptor cell mutants relative to WT lipid ratios presented in panel d.

When DGK is inhibited, plasma membrane PA decreases by 90%. This effect suggests that the inclusion of regulation of PIP5KI by PA in the plasma membrane might enable Suratekar’s model to replicate Almaça’s observations. We implemented this regulation in two ways: first, by changing the Vmax of PIP5KI from a constant into a Hill function of plasma membrane PA and, second, by making the KM of PIP5KI inversely proportional to plasma membrane PA. In both cases, the inhibition of DGK reduces PA in the membrane, and this reduction is transmitted to PIP5KI, resulting in a decreased rate of PI(4,5)P2 production. But, as PIP5KI loses efficiency, PI(4)P accumulates rapidly. This accumulation continues until the amount of PI(4)P compensates for the reduced efficiency of PIP5KI, thereby establishing a new steady state with a very elevated PI(4)P level that restores the levels of PI(4,5)P2.

Assigning a small rate constant (0.08) to the efflux from the pool of PI(4)P creates an “escape valve” that prevents PI(4)P from unduly accumulating when DGK is inhibited. This setting very slightly alters the steady state, but is easily balanced by increasing the endoplasmic reticulum PA source by 1.1%, which compensates for the new loss of material from the system. These settings lead to a model that is consistent with Suratekar’s data and replicates Almaça’s observations (Fig. 5b). More details on this matter can be found in Section 2.4 of the supplements.

Figures 5c, d permit an easy comparison with the analogous figures in the paper of Suratekar et al.31 (namely Figures S6, S8 and S9 in the supplements and Table 2 in the main text of their paper). Thus, one can readily confirm that the altered model satisfies the conditions that were required to validate the original model.

Taking all results together gives us confidence that our hypothesis of DGK moderating ENaC through PA control of PIP5KI convincingly explains the available data.

Discussion and Conclusions

The proposed work combines two models: one representing the dynamics of phosphoinositides and the other accounting for the regulation of ENaC and ASL. This combination allowed a detailed study of the intricate interactions between DGK and ENaC and yielded good consistency with available data.

Several results in the literature support our hypothesis that PIP5KI regulation by PA explains the moderation of ENaC activity induced by DGK knockdown, as Almaça et al.21 had observed. From the work of Moritz26, Jarquin-Pardo28 and Jenkins27 we know that reducing PA downregulates PIP5KI in vitro. The inhibition of DGK affecting PIP5KI is supported by in vitro data presented in our previous paper24 and by Luo, Prescot and Topham32 in Zebrafish, although this finding was contested by Jones, Sanjuan and Mérida33.

To the best of our knowledge, only Sandefur and colleagues34 made an attempt to study EnaC regulation by phosphoinositides. However, these authors did not consider SPLUNC1 and practically ignored phosphoinositides, only referring to them as mediators of P2Y2 purinoreceptor signaling, which is activated by extracellular adenosine triphosphate (ATP). This simplification has the crucial disadvantage that it becomes difficult to study the regulation of PIP5KI by DGK, which could be a promising drug target.

As mentioned before, several studies showed that PI(4,5)P2 and PI(3,4,5)P3 both affect ENaC, but we only focused on PI(4,5)P2. The rationale is that we are interested in human pulmonary epithelial cells, which are polarized, with a cell membrane that has different characteristics in different locations. So far, the data of highest quality point to a uniform distribution of PI(4,5)P2 throughout the plasma membrane and absence of PI(3,4,5)P3 in the apical part of polarized cells35, where CFTR and ENaC are located.

The past years have witnessed numerous discoveries regarding lipid transfer proteins (LPTs) (e.g.,36) and it is conceivable that a so-far unknown cellular mechanism could create the high sensitivity of PA production in the endoplasmic reticulum to the levels of PA in the plasma membrane, which would validate the hypothesis of Almaça et al. However, extensive, targeted research over several decades has not identified such a mechanism. It therefore appears that our hypothesis regarding ENaC regulation by DGK has a higher likelihood of being correct than the earlier hypothesis of Almaça and colleagues. According to this new hypothesis, which is supported by our computational analyses, the regulation of ENaC is primarily exerted through the control of PI(4,5)P2 production by PIP5KI, which in turn is controlled by PA, the product of the DGK reaction.

ENaC is composed of three subunits, which are usually called α, β and γ37. Channels with alternate stoichiometries have been reported to have very low activity38, so that some ENaC channels could be considered constitutively closed. Also, not all ENaC channels are necessarily equal, due to variability introduced by alternative splicing, alternative folding, glycosylation and ubiquitination. Thus, one might expect a range of ENaCs where, at one end, some channels are open even if ENaC controllers are signaling a closed configuration, and where the opposite is true at the other end. If there are indeed constitutively open and closed ENaC channels, the open probability function would have a higher minimum than the reported value, which is about 0.0239. By the same token, as a result of constitutively closed channels, the maximum should be lower than reported. To explore this situation, we tested an alternative ENaC Po function with a minimum of 0.12, a basal ENaC Po of 0.22, and maximum of 0.72. With this Po function, the model exhibit results that are similarly good as previous findings but yields much improved results with respect to perturbations in PLC activation. This alteration would likely make the activity of ENaC less sensitive to DGK inhibition and yield model results that are more similar to observations by Almaça and colleagues. Unfortunately, there are no data supporting this strategy.

It would be beneficial to study in greater depth the positive feedback loop between PI(4,5)P2, PLC, PA and PLD (Fig. 6), which was previously identified by van der Bout and Divecha40. Not only would it be interesting to see the regulatory capabilities of this functional arrangement, but one could also study intriguing observations like the one made by Antonescu et al.30, where the knock-down of an PLD isoform led to increased PA levels. At current count, there are ten DGK isoforms41,42, thirteen PLC isoforms43, two PLD isoforms44 and nine protein kinase C (PKC) isoforms. We did not account for this detail, but the diversity could be important for shedding light on some of the so far unexplained observations in the field. At the same time, the subcellular locations of PLD are still a matter of controversy and PLD may not be abundant in the plasma membrane45. Also, PLC activation is currently represented in a rather simplified manner, because the model does not include calcium (Ca2+) release. Thus, discrepancies concerning PLC are likely caused by a missing or so-far unknown mechanism that accounts for Ca2+ release and its feedback effect on PLC. Also, the time scale of our model is in minutes and much of the pertinent phenomena are better visualized in seconds, as can be seen in supplement Figure S4 and in other studies46,47. This aspect could be the topic of a possible refinement to be taken into account in the future. PI(4,5)P2 recovery after PLC activation is a complicated subject that has caused other researchers similar difficulties48,49, and the present model is also not able to explain it convincingly. However, when the PLC activation in our model is artificially sped up in a manner similar to what was proposed in the literature48,49, discrepancies with respect to PLC activation disappear. Further details are presented in the Supplement subsection on PI(4,5)P2 dynamics in PLC perturbations (1.3). Importantly, with or without this speed-up, the discrepancy has no real bearing on the main focus of our work, the regulation of ENaC by DGK.

Figure 6
figure 6

Hypothetical positive feedback regulation of PI(4,5)P2 and PA. In the plasma membrane, PA may be created from PC by PLD and by cleavage of PI(4,5)P2 by PLC. When PLC displays low activity, PA is mainly produced by PLD. When PLC activity increases, more PA is produced via degradation of PI(4,5)P2, which lowers the pool of PI(4,5)P2, secondarily reduces the activation of PLD, and ultimately leads to less PA production from PC. This dual activation mechanism could alter the ratio between PA coming from PLC and from PLD and contribute to the stability of PA and PI(4,5)P2 levels; however, confirmation of this mechanism will require further experimental investigation.

Our computational results point to the conclusion that DGK inhibition and the consequent decrease in PI(4,5)P2 levels moderate ENaC gain of function in CF by compensating for the higher number of ENaC channels induced by ASL dysregulation. Also, there is strong indication that this regulatory mechanism is mediated through the regulation of PIP5KI by PA. While altering the levels of PI(4,5)P2 could seem to be an interesting therapeutic target, caution is necessary as such alterations will probably have unpredictable and possibly undesirable side effects. For instance, Balla’s work15 shows that PI(4,5)P2 influences many proteins in the cell. In addition, alterations in the levels of PI(4,5)P2 obviously lead to changes in other phosphoinositide levels, such as PI(3,4,5)P3, which could have further ramifications, such as alterations in the AKT signaling pathway, to name just one. As an alternative, our work suggests that controlling PIP5KI activity or PA levels could achieve therapeutic benefits.

A particularly promising aspect of this research is the recent evidence that phosphoinositides act locally15,50 and are present in different pools with different functions35,51. This fact could possibly be exploited to create safe and effective therapies for CF and other diseases. Of course, before new therapeutic candidates are pursued, we need to understand far better how these local activities interact and how they are controlled. Such a deeper understanding will probably be gained through targeted combinations of future experimentation and modeling.

Methods

The equations and parameters of the ENaC—ASL module used here are exactly the same as those presented in29. The way in which the phosphoinositide model was designed is flexible enough to simulate membrane patches from different compartments if supporting data are available. For example, it is easy to block or greatly slow down a given reaction a priori if we know that it is not present in a particular compartment. Alternately, one would obtain very low enzyme activities for some reactions as a direct result of the parameter optimization to fit experimental observations. The expansion of the phosphoinositide pathway model introduced some new parameters that were not used in the earlier model24,29. For some of these new parameters it was possible to retrieve values from public databases (Table S2). The remaining parameters were adjusted to fit phosphoinositide steady-state values and other phenomena observed experimentally (described in24 and Table S1). The adjustment was performed gradually, first with manual adjustments, followed by Nelder-Mead optimization for each dataset separately, and finally using a hybrid genetic algorithm to identify the parameter set with the best fit for all calibration datasets [details in Supplements (Section 1)]. The extended phosphoinositide model was validated by comparing simulations with different datasets [Supplements (Secttion 1); Figures S3, S4, S5 and S6]. The combined model was validated using data from Almaça et al.21 on the activity of ENaC with or without DGK inhibition (Fig. 4).

Mathematical framework

A dynamical model of phosphoinositide metabolism was recently designed within the framework of Biochemical Systems Theory (BST)52,53,54,55,56,57,58, 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 Xi, 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 shown in Eq. (1).

$$ \begin{aligned} \frac{{dX_{i} }}{dt} = \sum\limits_{s} {n_{s \to i} .v_{s \to i} - \sum\limits_{p} {m_{i \to p} .v_{i \to p} } } \hfill \\ v_{i \to j} = \gamma_{i \to j} \cdot E_{i \to j} \cdot X_{i}^{{f_{i \to j} }} \hfill \\ \end{aligned} $$
(1)

Each quantity γij or Eij represents the rate constant or enzyme activity for a given flux, respectively, fi→j is the kinetic order, and nsi and mip are the stoichiometric coefficients for the influxes and effluxes. Each function v could contain additional modulators, again as power-law functions.

Model design, equations and parameters estimation

The model proposed here is a functional integration of two sub-models. The first in an extension of a phosphoinositide pathway model that was recently published, along with pertinent information regarding equations and parameter values24. The phosphoinositide model was mostly calibrated with information from the BRENDA59 and GENECARDS60 databases. As is typical, not all parameter values were available in these databases or the literature, and the missing values were therefore determined with a genetic algorithm such that the optimized model replicated data retrieved from the literature (cf. complete description in24). Influxes were associated with pools that have known sources outside the plasma membrane. For instance, PI is supplied from the ER via Nir261, and PI(4)P is delivered by vesicles from the Golgi and PI3P from vesicles from endosomes. The extension of the model furthermore includes influxes into PA from PLD and DAG. In the case of DAG, the PLC independent influx was necessary to maintain the levels of this lipid when PLC was not active.

The second sub-model captures the dynamics of ENaC—ASL and is described in29. All but two parameters of the ENaC—ASL regulation model where found by constraining the system using biological data from the literature. The two remaining undetermined parameters were determined through optimization1. Both models were designed within the framework of Biochemical Systems Theory (BST)52,54,55,56,62,63.

The main coupling point between the two sub-models is PI(4,5)P2, which functionally connects the phosphoinositide pathway with the dynamics of ENaC. The coupling suggests slight modifications to the previous phosphoinositide model24 which, in the final form, is depicted in the upper box of Fig. 1; the fluxes and equations are presented in Table 1, and parameters and initial values are given in Table 2. Table S1 of the Supplements summarizes experimental findings, which we used for parameter estimation, and references for parameters can be found in Table S2, as well as in our previous papers24,29.

In total, the two combined models have 15 dependent variables and 231 parameters and independent variables. We account for 25 enzymes, including three additions to the model (PLC, DGK and phosphatidate phosphatase). Phosphatidylserine (PS) was added as an independent variable. These new additions added 33 new parameters. Four new parameters were added to account for PIP5KI regulation (f4→45_PA, f4→45_PI45P2, f0→45_PA, f0→45_PI45P2). The ENaC/ASL model accounts for 2 dependent variables and 8 parameters.

Modified Suretekar’s phosphoinositide pathway model

To explore the feasibility of Almaça’s hypothesis regarding the regulation of ENaC by DGK, we used a model that in some sense represents a coarse alternative to our combined model. This model was proposed by Suratekar and colleagues31 and uses data from photoreceptor cells of Drosophila melanogaster. One must caution that this model, although addressing the same phosphoinositide pathway, may have features that are not entirely representative of human cells.

Suratekar and colleagues tested many versions of their phosphoinositide pathway model and ultimately decided on one that was in accordance with all data available to them. It contains an open cycle with influx of PA into the ER and efflux of DAG from the plasma membrane, but does not include PA regulation of PIP5KI. We used this version, implemented with Michaelis–Menten kinetics, as proposed by the authors.

In the original publication31, Suratekar’s model was considered to be sufficient to explain experimental observations if it satisfied the following five criteria: 1. in the wild-type steady state, the levels of the lipid intermediates (relative to PItotal) lie within 15% of the values given in Fig. 5d; 2. in the case of ten-fold under-expression of DAGK, corresponding to the rdgA3 mutant, the steady-state DAG/PItotal ratio should change by less than 15% compared to the wildtype value of this ratio; 3. under ten-fold underexpression of LAZA, corresponding to the laza22 mutant, the steady-state DAG/PItotal ratio should change by less than 15% compared to the wildtype value of this ratio; 4. under ten-fold underexpression of DAGK, the steady-state PAtotal/PItotal ratio should change by less than 15% compared to the wild-type value of this ratio; 5. under ten-fold underexpression of LAZA, the steady-state PAtotal/PItotal ratio should increase by 2.5-fold (± 15%) compared to the wild-type value of this ratio.

The steady-state levels of our model and Suratekar’s do not completely agree. In order to successfully link Suratekar’s and our ENaC-ASL model, we therefore divided the PI(4,5)P2 level by its steady-state value and multiplied it by 10,000. In this way, the steady state of PI(4,5)P2 becomes similar to the one considered by our model.

For the model to be able to replicate Almaça’s observations of DGK regulation of ENaC, we must make three alterations. The first is an implementation of the enzyme activity of PA where Vmax for PIP5KI is substituted with a Hill function (Eq. 2) which, in effect, makes the production of PI(4,5)P2 dependent on the levels of PA.

$$ \begin{aligned} V_{\max ,pip5KI} & = 0.3 + \frac{{V_{\max ,pip5KI*} \times PMPA^{2} }}{{K_{{M,pip5KI^{*} }}^{2} + PMPA^{2} }} \\ & = 0.3 + \frac{{1.51 \times PMPA^{2} }}{{{0}{\text{.01672982}}^{2} + PMPA^{2} }} \\ \end{aligned} $$
(2)

The shift constant 0.3 guarantees that PIP5KI is still active without PA, KM,pip5KI* corresponds to the steady-state level of PMPA (Plasma Membrane PA) and Vmax,pip5KI* was calculated for Vmax,pip5KI to have the same value as in the original model at the steady state of PMPA. We substituted Vmax,pip5KI in the flux equation for PIP5KI by this new expression.

Because this alteration would cause an explosive increase of PI(4)P when DGK is inhibited, we added an efflux from the PI(4)P pool. The efflux is shown in Eq. (3) and included in the differential equation for PI(4)P as a negative term.

$$ {\text{V}}_{PI4P\_exit} = 0.08 \times{\text{PI}}\left( {4} \right){\text{P}} $$
(3)

The value for the rate constant in this very small flux was obtained by trial and error until the model exhibited a behavior consistent with experimental observations. We are not aware of any direct biological evidence of a significant efflux of this type but it seems reasonable to assume that every phosphoinositide pool should have some efflux representing the phospholipids that exit the plasma membrane by vesicle or non-vesicle transport. It is not clear how relevant this postulated efflux is for the phosphoinositide pools, but we have implemented effluxes in our phosphoinositide pathway model.

Finally, although the steady-state levels of the model are not drastically perturbed by these alterations, we compensated for the new exit of material from the system by increasing the ER PA source flux by 1.1%. Again, trial and error were used to determine this parameter.

One should note that we are not trying to find an optimized set of parameter values for the model. Our objective is solely to test whether Suratekar’s model can qualitatively replicate Almaça’s observations when regulation of PIP5KI by PA is included.

Model implementation

The model was implemented in the programming language R v3.1.064 together with the package deSolve65. We used the ODE integration function with the LSODA method. The code used to implement the models is available in the Github repository with the URL “https://github.com/dolivenca/PI_ENaC_model”.