ENaC regulation by phospholipids and DGK explained through mathematical modeling

Cystic fibrosis is a condition caused by mutations in the cystic fibrosis transmembrane conductance regulator (CFTR). It is also thought to increase the activity of epithelial sodium channels (ENaC). The altered function of these ion channels is one of the causes of the thick dehydrated mucus that characterizes the disease and is partially responsible for recurrent pulmonary infections and inflammation events that ultimately destroy the lungs of affected subjects. Phosphoinositides are signaling lipids that regulate numerous cellular processes and membrane proteins, including ENaC. Inhibition of diacylglycerol kinase (DGK), an enzyme of the phosphoinositide pathway, reduces ENaC function. We propose a computational analysis that is based on the combination of two existing mathematical models: one representing the dynamics of phosphoinositides and the other explaining how phosphatidylinositol 4,5-bisphosphate (PI(4,5)P2) influences ENaC activity and, consequently, airway surface liquid. This integrated model permits, for the first time, a detailed assessment of the intricate interactions between DGK and ENaC and is consistent with available literature data. In particular, the computational approach allows comparisons of two competing hypotheses regarding the regulation of ENaC. The results strongly suggest that the regulation of ENaC is primarily exerted through the control of PI(4,5)P2 production by type-I phosphatidylinositol-4-phosphate 5-kinase (PIP5KI), which in turn is controlled by phosphatidic acid (PA), the product of the DGK reaction.


Motivation
Information pertaining to our original phosphoinositide pathway model ( Figure S1) was presented in 1 . The model mainly addresses concentrations of phosphoinositides and changes in their amounts in response to numerous perturbations of the pathway. These perturbations usually affected enzymes and their activators and inhibitors or consisted of alterations to influxes or effluxes of the system. The results of these experiments are summarized in the supplementary information of 1 . It reproduces several phenomena described in the literature that characterize the pathway, even for phosphoinositides present in small quantities. Notably, however, this model does not account for ENaC or other elements that are needed to study the regulation of this channel. The fluxes V0→45 and V45→0 were originally unknown but included to allow the model to replicate the observation that PI (4,5)P2 levels can be maintained even with low PI(4)P levels (Balla, 2013;Hammond et al., 2012). However, there is now biological evidence for the existence of both, as explained in our original paper about the phosphoinositide pathway model (Olivença et al., 2018).
The new additions to the model result in data fits of a similar quality as for the simpler version of the model, but allow us to assess a host of new features. In particular, they permit assessments of the effects of alterations in phosphoinositides on the dynamics of ENaC, which provide further details on the interactions between PI, PI(4)P, PI(5)P, PI(4,5)P2, PI4K and PIP5KI. Thus, we connect this model to a module containing all biochemical components that are necessary for studying ENaC regulation by PI (4,5)P2. Specifically, we add phospholipase C (PLC), which cleaves PI(4,5)P2 into DAG and inositol triphosphate (IP3), and allows for competitive inhibition of this phospholipase by PI(4)P and PS 2 . In addition to DAG and IP3, we explicitly define the enzyme DGK, which facilitates the transformation of DAG into PA. We also add phosphatidate phosphatases (LLPs) that hydrolyze PA back to DAG.
A new input flux of material to the PA pool now accounts for the transformation of PC into PA by phospholipase D (PLD). We opted for this implementation because PC exists in abundance in the cell membrane, with approximately 3,000,000 molecules/μm 2 , which corresponds to approximately 25% of the cell membrane 3 , and will therefore never be significantly depleted by the action of PLD. It also has a structural role in the cell membrane, and severe depletion of PC causes the cell to lyse. The flux rate constant is adjusted to include the contribution of PC.
We also include the observed activation of PLD by PI(4,5)P2 4,5 by letting the levels of phosphoinositide influence the conversion of PC into PA. This activation creates a positive feedback between PI(4,5)P2 and the production of PA by PLD, which enables the system very elegantly to increase the sensitivity of PI(4,5)P2 production to alterations in PA when PLC is active ( Figure 6). This mechanism appears to do the following: PA has two main sources. It may be created from PC by PLD, or it may be derived from the cleavage of PI(4,5)P2 by PLC. When PLC displays low activity, a good amount of PA is being produced by PLD. However, when PLC is activated, more PA is produced via degradation of PI(4,5)P2. Less PI(4,5)P2 reduces the activation of PLD and leads to lower PA production from PC. Not only does this dual mechanism contribute to the stability of PA levels, but it also alters the ratio between PA coming from PLC and PA coming from PLD.
The combined model moreover accounts for the activation of PIP5KI by PA (Figures 1, 2 and 6 in the main text), a regulation frequently described in the literature [6][7][8] . This regulation renders PI(4,5)P2 production dependent on DGK and allows us to study the activity of ENaC under different DGK levels of activation.
In addition to including these mechanisms, we introduce three minor refinements. First, we separate some enzymes catalyzing the same reaction into different fluxes when there is sufficient information to do so. For instance, we define two groups of the myotubularin phosphatase family: MTMR_1_6_14, which corresponds to myotubularins 1 to 6 plus 14, and MTMR78, which contains myotubularins 7 and 8. We cannot really separate these groups further into individual enzymes because no information is available for specific differences among MTMR's. For the same reason, ORCL1, INPP5 B/J, SKIP and SAC2 are grouped into O_I_SK_SA2. Outside these cases, the enzyme separation allows a better understanding of the mechanics in the model and a more accurate replication of experimental perturbations to the pathway, especially with respect to explorations that involve the inhibition or activation of enzymes.
Finally, we account for substrate competition in accordance to observations from wet lab experiments that reported the existence of multiple substrates, especially in the case of phosphatases 9 . The regulation among PTEN, PI3KI, PI(4,5)P2 and PI (3,4,5)P3 is implemented by creating pools of active (PTENa and PI3KIa) and inactive (PTENc and PI3KIc) PI3KI and PTEN and allowing the products of these enzymes to activate the enzymes. Parameters for this part of the model were found by fitting the model behavior to data reported by Feng et al. 10 , which were not used in the earlier phosphoinositide model.
As for the earlier models, rate constants and kinetic orders are derived from enzyme kinetic information provided by the BRENDA database 11 or from the literature, after appropriate unit transformations. In this model, KM values are given in molecules/µm 3 , specific activities in molecules/min/mg of enzyme and enzymes in mg. Details can be found in Table ST4 of our previous paper 1 . The kinetic parameters of enzymes used only in the combined model are presented in Table S2. Enzyme activities and transport fluxes were manually set to approximate reported phosphoinositide steady-state values. In some cases, the rate constants were estimated based on data fits. The parameters were derived from different data sources. Initially, we manually adjusted the parameters to fit the data coarsely. Next, for each dataset, a local, general-purpose optimization function implemented in the R language was used to search for the best fit; specifically, the Nelder-Mead method was used. In the final step, a hybrid genetic algorithm was implemented in the GA package to search for the global minimum.

PIP5KI activation by PA
PIP5KI is the enzyme responsible for phosphorylating the fifth carbon on the PI(4)P inositol ring and transforming it into PI(4,5)P2. It has two identified binding sites for PA in the C-terminal domain 8 . PIP5KI also has two binding sites for its substrate, PI(4)P, and these have different functions. One allows inhibition of the enzyme, whereas the other is a catalytic site that phosphorylates PI(4)P into PI(4,5)P2. According to Jarquin-Pardo et al., 8 , the catalytic site has a KM of 134 μM in the absence of PA, which is much higher than the inhibitory site with 2.4 μM. Consequently, the enzyme is inactive or in a reduced activity state if the concentration of PA is low. In the presence of PA, the situation is reversed: the KM for the catalytic site is 2 μM and the KM for the inhibitory site is 4.2 μM. The values of the kcat for both these cases are 16.4 and 17 min -1 . They suggest that PA will mainly alter the affinity of the catalytic site for PI(4)P.
To identify a function that describes the PIP5KI activation by PA, we consulted three papers containing data on PI(4,5)P2 production under different levels of PA 6,8,12 . These articles were described before, and we focus here specifically on PIP5KI activation. Unfortunately, all three studies used Triton X-100, which Jones et al. 13 describe as greatly inhibiting the basal lipid kinase activity. We explored different means of compensating for this inhibition. In Moritz' case, we altered the basal activity of PIP5KI from 1 to 1.5 molecules/min/unit of PIP5KI and scaled the data in such a manner that the maximum activation was 3, as reported by Jenkins and colleagues 7 . We achieved this scaling by multiplying the values of PI(4,5)P2 production with 3/20. For Jarquin-Pardo's data, we performed two transformations of the data. First, we altered the basal activity of PIP5KI from 0.2 to 1.5 molecules/min/unit of PIP5KI. Second, we scaled the data so that the maximum activation was again 3. In this case, we achieved this scaling by multiplying the values of PI(4,5)P2 production with 3/30. Finally, for Jenkins' data, we altered the basal activity of PIP5KI to 1.5 molecules/min/unit of PIP5KI and scaled the data so that the maximum activation was again 3. In this case, we multiplied the values of PI(4,5)P2 production with 2.2/13.

Creating a PIP5KI PA activation function
Because data presented by Moritz and Jarquin-Pardo both suggest a sigmoid function for the activation of PIP5KI by PA, we chose a shifted Hill function. This function (Eq. 1) has a KM of 10,000, a Hill coefficient of 2 and a maximal fold-activation of PIP5KI of 2; it is furthermore shifted by 1.
These settings ensure that the maximum activation of the PIP5KI enzyme is 3-fold (4.5 molecules/min/unit of PIP5KI), which is in line with the range of 2.2 and 3-fold activation reported by Jones et al. 13 and Jenkins et al. 7 . When PA has a concentration of 10,000 molecules/μm 3 , the enzyme is activated at half the maximum rate, and this value is within the physiological range for PA 3,6,14 . The Hill coefficient of 2 agrees well with the two reported binding sites for PA 8 . Model results are shown as lines and points represent data. Black symbols represent data used to estimate the sigmoid function in the model, which is shown in blue. Grey points represent Jarquin-Pardo's data with the estimated correction for the use of Triton X-100.

Moritz' data
Moritz et al. 6 studied PIP5KI activation as a function of PA concentration ( Figure S2) using Triton X-100 which, as Jones et al. 13 report, greatly inhibits the basal PIP5KI activity. This inhibition causes the reported basal activity to be lower and the apparent fold change activation to be increased. The authors report a basal activity of 1 pmol/min and a maximum activation of the enzyme of 20-fold. We adjusted for the effects of the detergent, as described below.
Jarquin-Pardo's data Jarquin-Pardo and colleagues 8 reported on an impressive amount of work addressing PIP5KI kinetics and its activation by PA. Unfortunately, they also used Triton X-100. Their data on PIP5KI activation as a function of PA concentration can be seen in Figure S2. The authors report a basal activity of 0.2 molecule/min/unit of PIP5KI and a maximum activation of the enzyme of 30-fold.
To determine the range of PA concentrations in Jarquin-Pardo's experiments, we estimated the total concentration of Triton X-100 they used. Our calculations suggested that a pure Triton X-100 solution at 0.1% volume will have a concentration of 1,600 μM. Because the Triton X-100 was in a mixture, we searched for a number that would produce a better agreement between Jarquin-Pardo's and Moritz' data.
We found that a Triton X-100 concentration of 35,000 μM would produce good agreement between the two datasets and was close to the calculated value. Furthermore, although we do not know the exact figure for PA levels, different sources report them as close to the levels of PI(4)P and PI(4,5)P2 3,6,14 , which is in good agreement with Moritz' data.

Jenkins' data
Jenkins et al. 12 studied PIP5KI activation by PA with and without Triton X-100. Unfortunately, they only showed the data for the experiment with the detergent. The authors report a maximum activation of the enzyme of 8, 15 and in some cases even 50-fold with Triton X-100. For the experiment without the detergent they report that the maximum activation was around 3-fold, which agrees with reports by Jones and colleagues. The data points attributed to Jerkins et al. in Figure S2 are the data points of the experiment with Triton X-100 scaled down such that the maximum activation of PIP5KI corresponds to the maximum activation observed in the experiment without Triton X-100.
Other pertinent data Pochynyuk et al. 15 provided information about the concentrations of PI(4,5)P2 under different levels of PLC activity and measured 15 minutes after the perturbations occurred. We used this information, as well as the levels of DAG, PA and IP3, for the parameterization of PLC activity in the model.
The paper by Sampaio et al. 3 contains values for PA and DAG levels that correspond to 0.5% -1% of total cell lipids in polarizing MDCK cells. This level makes them similarly rare as PI(4)P and PI(4,5)P2, which are found at around 10,000 molecules/μm 2 in our models, which replicate a 1 μm 2 cell membrane patch 1,16 . PA and DAG values were set in the model to 5,600 and 7,800 molecules/μm 2 , respectively.
Studies on platelets 14,17 suggest that the levels of IP3 are between 300 to 700 molecules/μm 3 ; the model uses a value of 500 molecules/μm 3 . This value was achieved by enabling fast removal of IP3 from the system, which is in agreement with the fact that IP3 is a soluble molecule.
Some components of the system display a large dispersion of parameter values, especially in specific enzyme activity. For instance, according to James et al. 18 , PLC is inhibited by many detergents commonly used in laboratory experiments, specially Triton X-100. James et al. used dodecylmaltoside, which was shown not to affect the activity of PLC. They report values between 0.11 and 0.18 μM for the KM and between 31.3 and 38.9 μmol/min/mg for the specific activity. We decided to use these parameter values but deemed the activity of PLC too strong. As a correction, we considered that PI(4)P competes for PLC, as shown by Ginger, Seifert and others 19,20 and that PLC is inhibited by PS 20 . Values for PS were taken from Sampaio et al. 3 .
Parameter values from BRENDA Several papers cited in the BRENDA database 11  For our model calibration, we ensured that all pertinent values fell into the ranges provided by BRENDA and obtained specific values through various fitting procedures.

Validation of the Expanded Phosphoinositide Pathway Model
Several datasets were retrieved from the literature and used to fit and validate the extended phosphoinositide pathway module. The same tests used in the previous version of the phosphoinositide pathway model 1 were also applied to the new version and the results are summarized in the first section of Results of the main text and in Figure 3.

PI(4,5)P2 dynamics in PLC perturbations
Xu et al. 21 , Gericke et al. 22 and Purvis et al. 14 state that PI(4,5)P2 recovers rapidly from an initially steep decrease in response to activation of PLC. The activation of PIP5KI by PA contributes to this recovery 6-8 . These reports do not quantify the magnitude of the recovery or the magnitude of PLC activation that allows a recovery, but it is known that a strong activation of PLC causes PI(4,5)P2 depletion 23 .
The model confirms PI(4,5)P2 recovery right after PLC activation. Figure S3 exhibits the time course of PI(4,5)P2 after a two-fold activation of PLC and the subsequent recovery. The original parameter set allows only a modest recovery while PLC is active, but a combination of stronger PA / PIP5KI activation and an adequate activation of PLC produces very clear PI(4,5)P2 recovery.  The recovery of PI(4,5)P2 after termination of PLC activation is much slower than reported in the literature. Specifically, Figure S3 suggests that the recovery of PI(4,5)P2 after the end of a two-fold activation of PLC takes approximately 10 minutes (reaching half of the recovery in 2.5 minutes), according to our model. To resolve the discrepancy between the model and the data, we followed Falkenburger et al. 2 and Xu et al. 21 , who accelerated the production of PI(4,5)P2, without actually mentioning a particular mechanism or rationale. Nevertheless, this speed-up resolved the discrepancy in time scale in their models.
The creation of PI(4,5)P2 can follow alternative routes, and we chose to accelerate the direct transformation of PI into PI(4,5)P2. To maintain the level of the lipid at the observed steady state value, we compensated by very slightly increasing the efflux of the PI(4,5)P2 pool.
In the top panel in Figure S4, with the previously mentioned changes, the model can now simulate PI(4,5)P2 recovery during PLC activation (black dots in the top panel in Figure S4 represent data from Chang et al. 24 ), with fast recovery after the end of PLC activation. Even with this accelerated PI(4,5)P2 production, the model faithfully replicates the data from Almaça et al., namely, the normalization of ENaC activity when DGK is inhibited ( Figure S4 bottom).
Unfortunately, these parameter alterations compromise the fit of model simulations with other experimental observations. PLC consumes PI(4,5)P2, so our intuition assumes that if we activate PLC, PI(4,5)P2 should decrease. It is possible to set PLC activations that are small enough so that the consumed PI(4,5)P2 will be less than the PI(4,5)P2 produced by the increase in PA and its activation of PIP5KI. As a consequence, small PLC activations cause the level of PI(4,5)P2 to increase. Because of this secondary effect, we did not use the PI(4,5)P2 production acceleration. Importantly, the speed of recovery is not important for the main focus of our analysis. To the best of our knowledge, there are no reports in the literature about this behavior. It is possible that the model actually predicts the effects of PLC on PI(4,5)P2 correctly but that researchers did not report them due to the small scale of the lipid fluctuations, considered them as artifacts produced by their experimental setup, or hesitated to mention them because the observations were contradictory to their intuition. Also, IP3 causes the release of Ca 2+ from the ER reservoir, which in turn activates PLC, which secondarily increases the production of IP3. We did not add Ca 2+ to the model because it would bring considerable complexity to this work and distract from our main message. Yet, it is possible that the missing components exert regulation suppressing these artifacts.
Interestingly, Antonescu et al. 25 reported a similar behavior between PLD and PA. They explain it by suggesting a negative feedback mechanism where the PA produced by PLD could negatively regulate the more efficient DGK PA production. They also highlight the fact that certain isoforms of DGK are negatively regulated by PA.
Both versions of the model replicate Almaça's observations of ENaC activity with DGK inhibition. Therefore, this discrepancy in the time scale of PI(4,5)P2 recovery after a PLC activity pulse does not affect our conclusions about the relevance of the PA/PIP5KI activation for the impact of DGK inhibition on ENaC activity.
If the activation of PIP5KI by PA is sufficiently strong, the model system becomes bistable, and the levels of the dependent variables do not return to the initial steady state, even if PLC activation is returned to its initial level. Instead, the system becomes locked into a state of high PA levels that increase PI(4,5)P2 production so much that the basal activity of PLC produces enough PA to keep the high PA levels and maintain the new steady state. This model result could be an artifact because PLC saturation would prevent high DAG and PA production only as a consequence of a great increase in PI(4,5)P2. However, if this bistability could be demonstrated in reality, it could be an effective way to maintain the signal after the initial stimulus has ceased. Interestingly, Purvis et al. 14 encountered the same phenomenon with their model. Pochynyuk et al. 15 measured the levels of PI(4,5)P2 for the first 15 minutes after altering PLC activity. The data and model results can be seen in Figure S5a. The model presents a good fit to these data. Furthermore, while no data exist for later time points, the simulation can be extended and shows that, after about 500 minutes, PI(4,5)P2 reaches a steady state that is increased to 1.6-fold for PLC inhibition and decreased to 0.4-fold for a 4-fold activation of PLC activity in comparison to the basal level.
Degradation of PI(4)P and PI(4,5)P2 Várnai et al. 26 presented data using reporters for PI(4)P and PI(4,5)P2 that derived from two experiments: one where SAC1 4-phosphatase was used to degrade PI(4)P and another where the INPP5E 5-phosphatase was used to degrade PI(4,5)P2. For the SAC1 experiment, our model produced a good fit ( Figure S5b) with an initial 8-fold increase in the model's SAC1 activity followed by attenuation to a 2-fold activation after 4 minutes.
Our model was initially unable to reproduce the INPP5E results when we used perturbations of a reasonable magnitude. This failure was probably due to the fact that we had underestimated the amount of enzyme in the cell membrane or the enzyme activity. Raising the activity of the enzyme 30-fold, we obtain noticeable alterations to the PI(4)P and PI(4,5)P2 levels. Furthermore, the shape of the data seems to suggest that the chimera INPP5E used by Várnai et al. was, in addition to its 5-phosphatase ability, able to hydrolyze the 4-phosphate of PI(4,5)P2. To test this speculation, we increased the activity of the SYNJ phosphatase, which is able to hydrolyze the 4 th and 5 th position of the inositol ring in PI(4,5)P2. With this alteration, the model led to a qualitatively reasonable time course ( Figure S5c). Zhong et al. 27 showed that Sac1 and other phosphates can have an allosteric effect which can essentially increase/decrease the Sac1 activity based on substrate concentration. The model results suggest that this may also be the case for INPP5E.
Perturbations in PTEN, PI3KI, PI(4,5)P2 and PI(3,4,5)P3 For completeness, we discuss fitting the earlier phosphoinositide model to data reported by Feng et al. 10 and Pochynyuk et al. 28 , which indicated that the specific activity of PTEN must be greater than the value retrieved from the enzyme database BRENDA 11 . This discrepancy is in line not only with Feng's conclusions, but also with reports by McConnachie et al. 29 and by Johnston and Raines 30 . We proceeded by using the values presented by Johnston and Raines 30 and also adjusted the amounts of total and active PI3KI to obtain an adequate quantity of PI(3,4,5)P3.
The data fits from our model are not particularly good, which is at least partially due to the fact that the system we are interested in is the apical part of the plasma membrane of bronchial epithelial cells, which is characterized by a lack of PI(3,4,5)P3. Nonetheless, we show the results of this simulation.
Feng et al. 10 presented two time courses that resulted from first activating PI3KI with rCD1 at time point t=5 min and then inhibiting the enzyme with FK506 at time point t=40 min. We simulated these perturbations with a 28-fold increase in the activation flux of PI3KI at t=5 min and a return of PI3KI to its steady-state values at t=40 min ( Figure S6a). Figure S6b shows results of blocking PTEN activity with H2O2 at t=2 min and inhibiting PI3KI with FK506 at t=8 min. We simulated these perturbations with a decrease of 50% in PTEN at t=2 min and a decrease of PI3KI of 30% values at t=8 min. The simulation results reflect the observations quite well.
Pochynyuk et al. 31 present a PI (3,4,5)P3 time course where PI3KI is inhibited with LY294002. We simulated this perturbation with a 90% decrease in PI3KI ( Figure S6 c). Again, our model was able to reproduce this experimental observation. While the fit in Figure S6a is not particularly good, Figures S6b and 6c are much better. It is true that the data suggest the existence of a mechanism that limits the increase in concentration of PI(3,4,5)P3 when it reaches about 120% of its normal level, but we are unaware of what this mechanism could be. Even so, we considered it pertinent to show Figure S6a to convey an impression of the limitations of the model. One should note that PI(3,4,5)P3 is not fundamental for our system which addresses the apical part of the plasma membrane of bronchial epithelial cells, which is characterized by a lack of PI(3,4,5)P3. In a future work, we intend to pay attention to the study of the regulation of ENaC by PI(3,4,5)P3, which is expected to improve the fit of the model to these particular data.
1.4. Previous modelling efforts to simulate the Phosphoinositide pathway.
Other models have been proposed to simulate the phosphoinositide pathway or at least parts of it. Falkenburger, Dickson and Hille 2,23 studied the kinetics of PI(4)P, PI(4,5)P2, PLC, DAG, IP3 and calcium. They showed that an acceleration of PI(4,5)P2 production during PLC activation is necessary. If PI(4,5)P2 production is not accelerated, PLC action will fade or stop. This study did not account for the activation of PIP5KI by PA, although it considered an increase in PI4K function when PLC was activated, which has a similar effect.
Narang et al. 32 developed a mathematical model that takes the activation of PIP5KI by PA into account. However, this model oversimplifies the phosphoinositide pathway by representing PI, PI(4)P, PI(4,5)P2 and PI(3,4,5)P3 with just one variable and implementing PA activation of PIP5KI with an auto-activation of the phosphoinositides.
Purvis et al. 14 built a model of the pathway in platelets. While generally impressive, this model does not account for all pertinent phosphoinositide species and does not include the activation of PIP5KI by PA.
PI is phosphorylated twice to become PI(4,5)P2: first by PI4K and afterwards by PIP5KI. The limiting rate of this process seems to be PI4K phosphorylation, which is 20 times slower than the phosphorylation of PIP5KI. In an attempt to explain this difference, Falkenburger, Dickson and Hille 2,23 suggested that an acceleration of the production of PI(4,5)P2, which would compensate for the depletion caused by PLC, should impact PI4K. We cannot fully agree with this rationale because of the following observations. PIP5KI is inhibited by its substrate, PI(4)P, which Falkenburger's argument ignores. Specifically, Jarquin-Pardo et al. 8 state that PIP5KI has an inhibitor as well as an activator binding site for PI(4)P and that, in the absence of PA, the KM for the inhibitor site is 50-fold smaller than the KM for the activator site. These authors furthermore state that "the addition of PA increases the affinity with which the active site of mPIP5K-Ib binds PI(4)P by 67-fold." Furthermore, Moritz et al. 6 report that the PIP5KI activation by PA is inhibited by PI(4,5)P2. These two pieces of evidence suggest that PIP5KI only functions fully in the presence of PA and with reduced levels of PI(4,5)P2. Finally, PI4K has PI as its substrate, which is 300 times more abundant than PI(4)P, the substrate for PIP5KI. These observations suggest that PIP5KI is the rate limiting component of PI(4,5)P2 production.
Falkenburger et al. 2 used in their model an amount of PI(4)P that is less than a third of the quantity of PI(4,5)P2 and assumed that two thirds of PI(4,5)P2 are in a bound state while all PI(4)P is free. We kept the amounts of these two phosphoinositides similar, in accordance with what is commonly accepted 33,34 .

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

ENaC control by DGK
Almaça et al. 35 performed siRNA screens to identify modulators of ENaC activity and noticed that DGK inhibition reduces ENaC activity to WT basal levels in CF F508del cells, but does not have a significant effect in healthy (WT) cells. Also, in CF, DGK has no effect on ENaC when PLC is activated or inhibited.
To explore these findings with our model, we started by studying the effects of SPLUNC1 and DGK perturbations. We simulated CF conditions by setting SPLUNC1 and the CFTR dependent flux (V5) equal to zero. As a result, the number of ENaC channels (N) increases 2.25 times relatively to the WT, which is close to the value reported by Tarran's group 36 . Concerning ENaC activity, ENaC is projected to be 2.25 more active in CF than in WT ( Figure 4). It should be noted that, in the model, the increase of ENaC activity under CF conditions is due to an increase in the number of channels. If we inhibit DGK in the CF setup, N is not affected but, due to the drop in Po, the activity of ENaC decreases to about half. This finding suggests that the moderation of ENaC by DGK inhibition is obtained by a different mechanism, such a s alteration of the channel-open probability, rather than the one being affected in the disease, namely the number of channels.

Comparison between WT conditions with and without DGK inhibition
Almaça et al. 35 showed that there are no significant differences between WT conditions with and without DGK inhibition.
To estimate the magnitude of DGK inhibition we used the information in the first paragraph of Almaça el al. 35 . On page 1394 they say: "First, we showed that a chemical blocker of DGK (DGKinh, 25 mM) led to 50% ENaC inhibition in A549 cells (Figures 4A and 4B)." We used this information to estimate the level of DGK inhibition and tested it in WT simulations of our model. The results can be seen in Fig. 4 of the main text. We choose 25% as a good substitute for values with several decimal places that would have yielded perfect 50% inhibition of ENaC.
The model, by contrast, indicates a clear drop when DGK is inhibited, which reduces the activity of ENaC by roughly 40%. Even so, these model predictions fall well within the confidence intervals of the data ( Figure 4). This effect is caused by a decrease in Po and a slight increase in N; the estimated drop in ENaC's activity is mainly due to the reduction in Po.

Comparison between WT and CF with DGK inhibition
Almaça et al. 35 report that inhibiting DGK in CF normalizes ENaC activity, i.e., the ENaC activity in these situations is similar to WT. The model yields some differences between WT, WT \ DGK-and CF \ DGKcases that should have roughly the same value ( Figure 4). However, the original data also show a nonsignificant difference between WT \ DGK-and the other two cases and all results are well within the confidence intervals in the data, so that we may consider these model results as an adequate reflection of the data. Overall, the model predicts a DGK-induced reduction of ENaC activity under CF conditions that is similar to the basal activity in WT.
Comparison between ENaC activity with and without DGK inhibition in CF, while accounting for inhibition of PLC Almaça et al. 35 found no significant differences for these cases. Confirming this result, the corresponding model settings produce the same value for the two situations ( Figure 4). Specifically, the model replicates a 1 μm 2 cell membrane patch where PI(4,5)P2 basal levels are around 10,000 molecules/μm 2 1,16 . When PLC is inhibited, the estimated levels of PI(4,5)P2 are around 16,000 molecules/μm 2 , which corresponds to an increase of approximately 50% over the basal level. The levels of PA, and even more so of DAG, are greatly reduced to 2,600 and 42 molecules/μm 2 , which corresponds to 33% and 1% of their basal levels, respectively. The severe reduction of DAG explains the lack of effect of DGK's inhibition over ENaC activity. With such small amount of substrate, the level of enzyme activity becomes irrelevant.
In the data from Almaça et al., the value for the PLC-scenario has the same magnitude as the unperturbed CF state (CF). This discrepancy is presumably due to the fact that Almaça et al. used a fluorescence essay where the marker was probably saturated for values higher than in the case of CF.
Some discrepancies remain when PLC is perturbed. In the case of PLC inhibition (CF\PLC-), Almaça's data indicate the same magnitude as in the unperturbed CF state (CF), while the model predicts a much higher value. This discrepancy is again presumably due to the fact that Almaça et al. used a fluorescence essay where the marker was probably saturated for values around 100% of ENaC activity. Nonetheless, the model is able to replicate the data associated with PLC inhibition at least qualitatively.

Comparison between CF with and without inhibition of DGK, while accounting for activation of PLC
According to Almaça's results, CF / DGK-, CF / PLC+ and CF / PLC+ / DGK-have similar values of about 50% ENaC activity in comparison with the unperturbed case of CF. These results are similar to those from the model (Figure 4). The model is able to replicate the CF / PLC+ case quantitatively; however, when the kinase is inhibited (CF / PLC+ / DGK-), a slightly more pronounced reduction in ENaC activity is predicted by the model. This reduction is probably due to missing regulation that is not yet described in the literature and therefore not implemented in the model, which suggests the need of more detailed studies of PLC and DGK. Notwithstanding this slight discrepancy, the model captures the qualitative behavior of the biological system quite well.
2.2. The combined model with phosphoinositide recycling, but without PIP5KI activation by PA, does not replicate the effects of a DGK knockdown.
The model allows us to test the hypothesis of Almaça and colleagues that the moderation of ENaC activity by DGK is caused by a general decrease in phosphoinositides, due to their reduced recycling. To assess this hypothesis, we alter the extended phosphoinositide model in two ways. First, we allow 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. This additional, simplified step in a crude way closes the circle of phosphoinositide recycling. Because the efflux of PA is only 5% of the influx of PI, we maintain the already present influx into the PI pool, which is now reduced by an amount equal to the PA efflux. The flux V➔0 in this model extension can be interpreted biologically as the amount of PI produced in the ER from PA that did not originate in the plasma membrane; expressed differently, it represents PA created de novo. Second, we remove the regulation of PIP5KI by PA.
With these presumably reasonable settings, the model does not replicate the observation by Almaça and colleagues: if DGK is inhibited, it no longer affects ENaC activity, whether in WT or CF (Figure S7 a). The reason is that the DGK independent V➔0 flux is much greater than that DGK dependent flux VPA➔.
To remedy the situation, it would be necessary to balance the input into PI and the efflux out of PA. Specifically, the input into PI would have to be an exact multiple of the efflux VPA➔ out of PA such that it would be equal to the original PI influx. Such a coupling would make the PI pool dependent on the PA efflux. In a simulation with these settings, the ENaC activity is more sensitive to DGK inhibition, but the decrease in the activity is very small (Figure S7 b).
If the inhibition of DGK in the model is made stronger, this configuration of the model can actually replicate the data reported by Almaça and colleagues ( Figure S7c). However, to achieve this configuration, PI influx would have to change proportionally to the PA efflux, which is 20 times smaller. Another possible explanation could be that the levels of PA and PI in the plasma membrane would somehow have to be tightly coordinated. There is no evidence of this degree of coupling between the PI influx and the PA efflux, and while there are no exact measurements for these lipids in the plasma membrane, all existing cell measurements suggest that PI is 10 times more abundant than PA 3,37 . Moreover, it is known that PA is produced de novo in the ER, and this production affects the sensitivity of PI to plasma membrane PA. Thus, the model does not support the hypothesis of Almaça et al. 35 that inhibiting DGK brings the recycling of the phosphoinositides to a halt, which in turn decreases PI(4,5)P2 and PI(3,4,5)P3 and causes reduced ENaC activity. , and additional 75% inhibition by DGK instead of the usual 25%, the model does replicate the observation of Almaça and colleagues. However, this degree of sensitivity of PI to plasma membrane PA efflux is contrary to empirical evidence.

Sensitivity Analysis
The results of a comprehensive sensitivity analysis are presented in Table S3. The box at the lower right corner of this  16 . As expected, these values are similar because this model was adopted without alterations.
Many of the high sensitivities are associated with the new additions to the model (PLC, DAG , PA), which is consistent with the observation that signaling systems may present high sensitivities 38 .
As in Olivenca et al. (2018), we found high sensitivities with respect to the V0➔45 flux parameters which, by affecting PI(4,5)P2, affect all other variables that are heavily dependent on lipids, including ENaC and ASL. In this new model, but not in the previous version, we also found high sensitivities with respect to the V45➔0 parameters.
The presented results are not sensitive to alterations in the initial conditions. The steady state is stable and given enough time, the system returns to this steady state. ENaC, ASL and PI (3,4,5)P3 are slower to recover.

Considerations about PI(4)P accumulation in Suratekar et al. model
In the fourth section of Results in the main text we link our ENaC / ASL model to the model of the recycling of phosphoinositides created by Suratekar et al. When DGK was inhibited, it caused a sizable increase in PI(4)P. To mitigate this we created a small efflux in the PI(4)P pool that function much like as escape valve, and prevent the dramatic increase. Here we discuss this situation further.
There was the possibility that a PIP5KI isoform insensitive to PA could function as the PI(4)P escape valve. But according to Sasaki et al. (2009), all PIP5KI isoforms are activated by PA. The phosphatidylinositol phosphate kinases type II (PIP4KII) transform PI(5)P into PI(4,5)P2. It is possible that a PA insensitive phosphorylation of PI(4)P occurs but there is no evidence for this insensitivity in the literature.
However, as a thought experiment, let's imagine the suggested PA-insensitive phosphorylation of PI(4)P does exist. To simulate this situation, we transform the PI(4)P escape valve flux into a flux from PI(4)P to PI(4,5)P2 that is insensitive to PA, thereby creating two PI(4,5)P2 influxes, one sensitive and the other insensitive to PA. The steady state is very slightly altered but the difference is at the order of 10^-4. Let's perturb the system. Reducing the Vmax of DGK to 20% does yield a similar dynamics in comparison to the system with the PI(4)P escape valve. But can this system show a dynamics like the one reported by Almaça et al.. Namely, will the value of PI(4,5)P2 drop and inhibit ENaC activity, when DGK is inhibited? The answer is no ( Figure S8). The extra PI(4,5)P2 produced by the PA-insensitive phosphorylation of PI(4)P is enough to maintain ENaC activity. Thus, a system with PA-insensitive phosphorylation of PI(4)P does not replicate the data from Almaça et al.

Figure S8. Creating a parallel PI(4)P to PI(4,5)P2 flux that is insensitive to PA will not reproduce the data reported by Almaça et al. that ENaC activity will decrease when DGK is inhibited.
Nonetheless, these results made us revisit our implementation of Suratekar's model. We realized that the escape valve flux works so well because it follows a first-order kinetics, and therefore does not saturate. What happens if the process saturates? Assuming a Michaelis-Menten representation, the answer depends on the Vmax of that flux. If it is high enough it prevents the system from being flooded with PI(4)P. Because an escape valve with an inadequate Vmax would not accomplish its task, we can assume that Vmax is relatively high.
What happens if the PA-insensitive phosphorylation of PI(4)P saturates as well? Let's assume the PAinsensitive phosphorylation of PI4P flux is similar to the PI(4)P escape valve at steady state, that a Michaelis-Menten model is a good representation of its dynamics and that the Km is equal to that of the PA sensitive phosphorylation flux of PI(4)P. If so, the results are again a very similar to the steady states of the original system. As before, the extra PI(4,5)P2 produced is sufficient to maintain the levels of ENaC activity when DGK is inhibited.
These results are not favorable to the hypothesis of the existence of a PA-insensitive phosphorylation of PI(4)P. However, we should be aware that other possibilities can cause the PI(4)P overload, like Sac1 activity or additional regulators of PIP5KI, just to name a few.  Figure S5 a, blue line and points.

2
In the apical part of the cell membrane, PLC activation will decrease PI(4,5)P2 by 50% in the first 15 min. 15 Figure S5 a, red line and points. 3 Intense activation of PLC can lead to 75% 22 or 90% 39 of PI(4,5)P2 depletion.
Increasing 7 times the basal activity of PLC will deplete PI(4,5)P2 by 73%, but PLC must be increased 20 times to deplete PI(4,5)P2 by 90%. 4 Inhibition of class I DGK in BSC-1 cells results in a 50% reduction in total PA levels, indicating the majority of cellular PA is synthesized by type I DGKs. 25 Inhibiting DGK by 75% will reduce PA in 52%.