- Research Article
- Open Access
A computational model of PKD and CERT interactions at the trans-Golgi network of mammalian cells
BMC Systems Biology volume 9, Article number: 9 (2015)
In mammalian cells protein-lipid interactions at the trans-Golgi network (TGN) determine the formation of vesicles, which transfer secretory proteins to the cellular membrane. This process is regulated by a complex molecular network including protein kinase D (PKD), which is directly involved in the fission of transport vesicles, and its interaction with the ceramide transfer protein CERT that transports ceramide from the endoplasmic reticulum to the TGN.
Here we present a novel quantitative kinetic model for the interactions of the key players PKD, phosphatidylinositol 4-kinase III beta (PI4KIII β) and CERT at the TGN membranes. We use sampling-based Bayesian analysis and perturbation experiments for model calibration and validation.
Our quantitative predictions of absolute molecular concentrations and reaction fluxes have major biological implications: Model comparison provides evidence that PKD and CERT interact in a cooperative manner to regulate ceramide transfer. Furthermore, we identify active PKD to be the dominant regulator of the network, especially of CERT-mediated ceramide transfer.
All transmembrane proteins and soluble proteins that are secreted from a cell into the extracellular space generally pass through the endoplasmic reticulum (ER) and Golgi complex (GC). In mammalian cells the GC constitutes a highly ordered organelle composed of cisternal stacks, the main function of which is the modification of proteins and sorting to their distinct membrane localizations . Cargo is received from the ER at the cis-Golgi complex, traverses the medial Golgi stacks, then reaches the trans-Golgi side where vesicles bud off from the trans-Golgi network (TGN) [2,3]. The formation of these vesicles is driven by a complex interaction of proteins and lipids and is still only partly understood. One of the key players of this network is protein kinase D, which comprises a family of three closely related serine-threonine kinases (PKD1, PKD2 and PKD3). PKD has an established role in the regulation of constitutive but also regulated vesicular transport processes at the TGN .
In the following, we shortly describe the molecular interaction network of PKD at the TGN (Figure 1 and reviewed in [5,6]). The first TGN-localized PKD substrate identified was PI4KIII β whose phosphorylation on serine 294 increased the activity of the lipid kinase, resulting in the enhanced production of its substrate PI4P . This phosphorylated lipid is especially important in defining Golgi membrane identity and providing a local signaling platform to which several PI4P-binding proteins are recruited, including the ceramide transfer protein CERT [8,9].
Ceramide is produced at the ER and then transported to the TGN in a non-vesicular manner by CERT, which is recruited by the ER-resident membrane protein VAP [10-12]. The mechanistic understanding is that CERT takes up a ceramide molecule into its START domain, which forms a hydrophobic pocket, followed by the release of ceramide at the TGN membrane . This transport is thought to occur at so-called membrane contact sites (MCS) where the ER and the TGN come into close contact . At the TGN, sphingomyelin synthases (SMS) catalyze the conversion of ceramide and phosphatidylcholine (PC) to sphingomyelin (SM) and diacylglycerol (DAG) [14-16]. DAG fulfills several central functions at the TGN by activating novel PKC isoforms which phosphorylate and activate PKD, by recruiting and activating PKD, and by directly impacting membrane curvature involved in the process of vesicle formation [3,17-19]. CERT, in turn, is phosphorylated by PKD at serine 132 , which serves as a priming site for multiple phosphorylations by CKI γ2 in the serine rich (SR) motif . This hyperphosphorylation of CERT negatively affects its affinity for PI4P and induces a conformational change that inhibits START domain function . At the ER, CERT is dephosphorylated by PP2C ε , enhancing the interaction of CERT with the Golgi membranes and also with VAP . Thereby, the proteins PKD, PI4KIII β and CERT are involved in a complex scenario of interrelated feedback loops . Despite the knowledge on the qualitative effects of phosphorylation on CERT membrane binding and function, it is currently still unclear how this complex molecular interaction network impacts the efficacy of ceramide transport at ER-Golgi MCS.
Only few studies have attempted to address aspects of TGN function by mathematical modeling. The existing models describe mechanisms for vesicle kinetics, membrane physics, and feedback in the SMS reaction [24-27]. Thus far, no quantitative model exists that describes the kinetics of key protein interactions at the TGN important for secretory transport. Such a quantitative dynamic model would contribute to an increased understanding of molecule interactions and organelle function and is therefore valuable to basic research. In addition, such a model has broader relevance, as the vesicular transport from the TGN to the plasma membrane drives the polarized migration of cancer cells. A further application are biotechnological interventions in the secretory pathway targeting the optimization of the production of therapeutic proteins in mammalian cells .
Here, we use mathematical modeling and perturbation experiments to address how the feedback loops within the PKD-CERT network interact to coordinate ceramide transport. We present a quantitative dynamical model for the molecular interactions at the TGN that is based on chemical reaction kinetics. Biochemical time series experiments after perturbation of the system and absolute protein quantification measurements are used for model calibration. As the information in the data is not sufficient to identify all model parameters uniquely, we use statistical Bayesian approaches, which are computationally demanding but particularly tailored to deal with this problem [29-31].
Our calibrated model is able to capture dynamic interactions between PKD, PI4KIII β and CERT on an average cellular level. Furthermore, validation experiments confirm that we are able to reliably predict different perturbation scenarios. Our model-based analysis thus provides insight into the regulatory network of key components underlying ceramide transfer at ER-Golgi MCS.
Plasmids and reagents
Expression plasmids encoding Flag-tagged CERT as well as EGFP-tagged CERT, PI4KIII β, and PKD2 are described elsewhere [7,20]. The antibodies specific for PKD1 autophosphorylation at serine 910, phosphorylated serine 294 in PI4KIII β and phosphorylated serine 132 in CERT have been described elsewhere [7,32,33]. Commercially available antibodies used were as follows: anti-GFP mouse monoclonal (Roche Diagnostics), anti-PI4KIII β rabbit polyclonal (Millipore), anti-CERT rabbit polyclonal antibody (Bethyl Laboratories), anti-PKD2 rabbit polyclonal antibody (Cell Signaling), anti-Flag M2 mouse monoclonal antibody (Sigma-Aldrich), anti- α-tubulin mouse monoclonal antibody (Millipore). Recombinant purified GFP and GST-PI4KIII β were obtained from Vector Laboratories and Life technologies, respectively.
HEK293T cells were maintained in RPMI 1640 medium supplemented with 10% FCS. For transient transfections, 3·105 HEK293T cells were seeded per well of a 6-well plate. The next day, cells were transfected with 2 μg of plasmid DNA using TransIT293 (Mirus Bio Corporation, Madison, WI) according to the manufacturer’s instructions and grown for 24 hours. Treatment of cells with kbNB142-70 (Tocris) and PDBu (Sigma-Aldrich) was at 10 μM and 1 μM, respectively.
Protein extraction of cells
Whole cell extracts were obtained by solubilizing cells in lysis buffer (20 mM Tris pH 7.4, 150 mM NaCl, 1% Triton X-100, 1 mM EDTA, 1 mM ethylene glycol tetraacetic acid (EGTA), plus Complete protease inhibitors and PhosSTOP (Roche Diagnostics)). Lysates were clarified by centrifugation at 16,000g for 15 min.
Proteins were separated on a precast 4-12% Bis-Tris polyacrylamide gel (Life Technologies) and blotted onto nitrocellulose membranes (Pall, Dreieich, Germany). After blocking with 0.5% blocking reagent (Roche Diagnostics), filters were probed with specific antibodies. Proteins were visualized with IRdye-coupled secondary antibodies on a Odyssey scanner followed by analysis with Odyssey software (LI-COR Biosciences).
Image data quantification
Blots have been scanned with the LI-COR ODYSSEE Infrared imaging system. Chanel gain was set to intensity value 5, software gain was varied from 3 to 8 depending on the antibody in use. Gain information is contained on all attached image files. The ’*.tif’ image data was exported with the LI-COR software. For quantification a MATLAB 2011b script and the software’s image analysis toolbox has been used. The script follows recent findings of densitometry research  like individual lane background correction and a unified sample tool geometry of one third of average signal width. The image analysis script can be supplied by the authors upon request.
All calculations are performed using MATLAB R2011b scripts. Models are implemented using MATLAB toolboxes SBtoolbox2 and SBPD package . The ODE solvers embedded in the toolboxes are adopted from the SUNDIALS suite . Relative and absolute solver accuracy has been set to 10−8. All MATLAB scripts are available from the authors upon request.
Bayesian inference strategy
Here we give a short overview of the implemented Bayesian inference strategy, for a more theoretical introduction in ODE modeling and Bayesian inference we refer to Section Background in the supplementary text see Additional file 1. We first perform several repeated minimizations using a bounded local minimization strategy with uniformly sampled initial conditions to gain a set of maximum a posteriori parameters estimates on a logarithmic scale. We then create a log uniform prior distribution for the Bayesian inference centered several orders of magnitude around the best found values (see Additional file 1: Section 1.2). Three repeated runs of a customized version of the parallel tempering MCMC algorithm presented in  are performed using ten temperature levels. In each run half a million samples are drawn from the posterior distribution for all temperatures, using MCMC starting values drawn from the best found parameter values from the initial optimizations. Negative results for non-convergence are successfully tested for all sampling runs based on intra-inter-chain variance . A representative parameter subsample (Additional file 2) is used to generate the model predictions. Validation data has not been used for model calibration. Bayes Factors for model comparisons (and their uncertainties) have been calculated for each of the three independent runs via thermodynamic integration .
Data-driven modeling of PKD, PI4KIII β and CERT at the TGN
The exact mechanism by which CERT transfers ceramide from the ER to the TGN is currently still unclear [6,10,13]. In particular, two models have been discussed (see e.g. Figure six in ). The ’short distance shuttle’ model suggests that CERT travels 10 nm or more through the cytosol, shuttling between ER and TGN membranes during the ceramide transfer process . By contrast, the ’neck-swinging’ model suggests that CERT is simultaneously bound to both membranes, while extracting ceramide from the ER membranes via its START domain and releasing it at the TGN membrane [2,39]. There are several regulatory differences underlying these distinct transport mechanisms, which we considered by building two different models, A and B: First, PKD-mediated phosphorylation affects CERT transport activity differently in both models. While it has a positive effect on CERT transport activity in the shuttle model, it decreases transport activity in the swinging neck model. Together with an overall positive effect of CERT transport on PKD activity, this constitutes a positive feedback loop for model A and a negative feedback loop for model B, respectively, as illustrated in Figures 1A and 1B. Second, CERT is regulated in a cyclic reaction scheme in the shuttle model, driven by active PKD and active PI4KIII β. By contrast, in model B, CERT transfer activity is switched on and off by PI4KIII β and PKD, respectively.
Figure 1 shows a schematic scheme of the differences between model A (top) and model B (bottom). Solid lines represent conversions of molecules or changes in localization, whereas dashed lines refer to regulatory or catalytic influences. For the sake of simplicity, some regulations are represented as effective direct regulations which summarize several intermediate reactions.
The figure also shows molecular targets of the perturbation experiments that are used for model calibration, denoted by inputs u (ectopic expression of PI4KIII β (u 3) and CERT (u 4)), as well as activation and inhibition of PKD (u 1 and u 2, respectively). Moreover, measurable outputs y are relative phosphorylation states of PKD (y 4), CERT (y 6) and PI4KIII β (y 5).
In both models, active TGN-bound PKD (PKDpDAG) phosphorylates and activates PI4KIII β (R1) and, by direct phosphorylation (R2), triggers the detachment of CERT from TGN membranes (R3). At the same time, active PI4KIII β produces PI4P (R4), recruiting CERT back to the TGN (R5). In turn, CERT mediated ceramide transfer has an activating influence upon PKD, which comprises SMS driven conversion of PC and ceramide into SM and DAG (R6) and subsequent DAG-mediated recruitment (R7) and activation (R8) of PKD.
In model A (Figure 1A), CERT is permanently associated with the ER via interaction with VAP-A. After delivery of a ceramide molecule to the TGN, PKD-mediated phosphorylation of CERT leads to dissociation from the TGN (R3) while association with the ER remains. Once dephosphorylated at the ER (R9), CERT can deliver a new ceramide molecule to the TGN (R10).
In model B (Figure 1B), the PKD-mediated phosphorylation of CERT leads to its detachment from both the TGN and ER membranes (R3). Thus, only the unphosphorylated CERT molecule is associated with the TGN and the ER and therefore transfer-active. Dephosphorylation of CERT is required for MCS re-recruitment.
Based on the current state of knowledge about these molecular interactions and on available measurements, we used chemical reaction kinetics to build differential equation models for the two model variants, which resulted in systems with seven state variables and 26 (model A) and 27 (model B) parameters, respectively. Model equations and details are given in the Additional file 1: Section 1.3 and Table S1.
Absolute measurements for model calibration
Calculation of absolute protein amounts per cell is a prerequisite to establish a quantitative model . To do so, we used an indirect method making use of a commercially available recombinant GFP standard. Decreasing numbers of recombinant GFP molecules were detected by Western Blot analysis using a GFP-specific antibody (Figure 2A, top). The band intensity was quantified to generate a standard curve connecting signal strength and protein amount (Figure 2A, bottom). Next, we transiently transfected HEK293T cells with a plasmid encoding GFP-tagged PI4KIII β. Cells were lysed 24 hours later and expression of GFP-tagged PI4KIII β was detected by quantitative Western Blot analysis using the GFP-specific antibody. In parallel, cell lysates were probed with a PI4KIII β-specific antibody, detecting endogenous and ectopically expressed PI4KIII β at the same time (Figure 2B, top). Based on the GFP signal, we determined the amount of GFP-PI4KIII β expressed in a defined number of cells by inverse regression using the GFP standard curve (for details see Section 2 in the Additional file 1). Subsequently, we calculated the absolute amount of endogenous PI4KIII β, which accounted to approximately two million molecules per cell (Figure 2C) by determining the signal ratio between endogenous and ectopically expressed protein from Figure 2B. We then validated this result with a commercially available purified recombinant GST-PI4KIII β protein. Similarly to the GFP standard, we detected decreasing amounts of GST-PI4KIII β by Western Blot analysis using a PI4KIII β-specific antibody (Figure 2B, bottom) and correlated the signal intensity with the number of molecules. Indeed, the results were in the same order of magnitude and in good agreement with respect to the error estimates of both methods thus confirming our calculation (Figure 2C). Likewise, we calculated the amount of endogenous and ectopically expressed PKD and CERT and integrated these absolute estimates together with the estimated standard deviation into our dataset. Our results show that a single cell contains roughly half a million PKD and CERT molecules each. The numbers of the ectopically expressed proteins ranged from ten (GFP-tagged PKD1) to hundred million (GFP-tagged CERT) molecules per cell (Figure 2C). Coefficients of variation are determined to be between 25 to 40%, accurate enough to estimate the order of magnitude of a protein.
Time series measurements for model calibration
In order to identify perturbations leading to a significant system response, we performed experiments in which we detected dynamic changes of the phosphorylation of PKD and its substrates PI4KIII β and CERT. Since the antibodies specific for phosphorylated CERT (pS132) and PI4KIII β (pS294) are not suited for the detection of the endogenous proteins we ectopically expressed CERT and PI4KIII β while additionally perturbing endogenous PKD activity by adding a respective activator and inhibitor to the cell cultures. First, we ectopically expressed GFP-tagged PI4KIII β in HEK293T cells for 24 hours and subsequently left cells untreated or stimulated the cells with Phorbol 12,13-dibutyrate (PDBu), a potent activator of PKD , see Figure 3A. Cells were harvested at different time points, and the level of active PKD was analyzed using an antibody specific for autophosphorylated PKD (pS910). In parallel, we detected the expression and the PKD-mediated phosphorylation of GFP-PI4KIII β using the GFP-specific and pS294-specific antibodies, respectively. Stimulation of cells with PDBu steadily increased PKD activity within one hour. Likewise, PKD-mediated phosphorylation of PI4KIII β was strongly enhanced during this time.
In a second experiment, we ectopically expressed Flag-tagged CERT, see Figure 3B. After 24 hours, cells were left untreated or treated with the selective PKD inhibitor kb-NB142-70  for the indicated time. Activation of PKD as well as phosphorylation and expression of Flag-CERT were analyzed using pS910-specific and pS132 and Flag-specific antibodies, respectively. Expression of Flag-CERT strongly increased endogenous PKD activity, as can be seen from comparison of the ’0’-lanes in Figure 3 B. This activation, however, was decreased upon treatment of cells with kb-NB142-70. The strongest reduction in PKD activation was detected after treating cells with kb-NB142-70 for 3 hours. Likewise, PKD-mediated basal phosphorylation of CERT was reduced in kb-NB142-70 treated cells compared to untreated cells, see Additional file 1: Section 3.1 and 3.2.
Bayesian parameter estimation
Both datasets (Additional file 3), absolute protein amounts and relative phosphorylation states after the perturbations described were used for model calibration in a statistical Bayesian framework. Therefore, we interpreted the measurements as random variates that were generated by a stochastic process, which is defined as the deterministic solution of the differential equation model perturbed by additive stochastic noise terms. The resulting stochastic model allows us to define the likelihood function, which assigns for each set of parameters θ the probability p(y|θ) to obtain the dataset y. Especially in a Bayesian approach parameters are interpreted as random variables as well. Prior knowledge about these are included via prior distributions p(θ) in the parameter space, which reflect the state of knowledge about the parameters before having seen the data. The objective in a Bayesian approach is the posterior distribution p(θ|y), the conditional distribution in the parameter space after having taken the data into account, which is according to Bayes’ Theorem proportional to the product of the prior distribution and the likelihood function, see Additional file 1: Section 1.1 and 1.2.
We investigated this distribution by generating samples via Markov Chain Monte Carlo sampling (Additional file 1: Section 4 and 5). These samples were subsequently used to estimate distributions for various quantities via Monte Carlo integration. Figure 4 shows a comparison of the calibrated model A with the set of experimental data that was used for calibration. Absolute endogenous protein amounts and respective distributions predicted by the model are shown in Figure 4A. As can be seen, both agree very well for all three proteins.
The data shown in Figures 4B-E correspond to time series measurements after system perturbations. Experimental data are represented as mean with standard deviations, while respective model predictions represent time-varying densities for the states of the system, which are illustrated by different colors (density is increasing from yellow to dark red). The maximum a posteriori estimate that maximizes the posterior distribution is depicted as solid black line. Outer lines represent the 99.5% confidence intervals estimated from the sample.
The data shown in Figure 4B correspond to the perturbation experiment in Figure 3A, i.e. ectopic expression of PI4KIII β and subsequent activation of PKD with PDBu. The absolute amount of PI4KIII β molecules per cell were quantified as described after 24 h vector expression, followed by time series experiments of phospho-PKD and phospho-PI4KIII β within a few hours after treatment with PDBu. The model suggests that PI4KIII β rapidly reaches a new steady state after a fast increase within 24 hours. Furthermore, PKD activity rapidly rises after treatment with PDBu, followed by PI4KIII β. Both active kinases reach a new steady state approximately 1-2 hours after treatment.
The model fit to the perturbation experiment in Figure 3B, ectopic expression of CERT and subsequent inhibition of PKD, is shown in Figure 4C. In contrast to Figure 4B, in which PI4KIII β was ectopically expressed, the model suggests that CERT still increases after 24 hours of vector expression. Inhibition of PKD causes a rapid decrease in phospho-PKD, followed by a decrease in phospho-CERT, which is well reflected in our model. Our simulations suggest that this is a transient decrease due to the continuing increase in the overall CERT amount.
Figure 4D depicts model-based refinement experiments, which were particularly chosen in order to reduce the variance in model predictions, as described in . These are PKD activation early after ectopic expression of PI4KIII β, and the long-term behavior of PKD 30 hours after treatment with PDBu.
Overall, Figure 4 shows a very good agreement of experiments and predictions of model A for all scenarios. The figure also shows that the model uncertainty is significantly reduced for time points near experimental data. Moreover, while some uncertainty remains for variables and time points far away from measurement instances, the flexibility of the model was significantly reduced by the data, which led us conclude that the model is not too complex and the data contain considerable information about model parameters.
The respective figures for model B are shown in Additional file 1: Figure S6 and Section 6.2. The data fit of model B shows two major differences. First, model B is not able to reproduce the experimentally observed early increase of PKD activity (Additional file 1: Figure S6D, in comparison to Figure 4). Second, model B suggests an overshoot of phosphorylated CERT at the beginning of the ectopic expression of CERT in Additional file 1: Figure S6C, which was not the case for model A (Figure 4C).
In addition to this visual inspection, we decided to use the Bayesian framework for a sophisticated model comparison that also takes differences in model complexity into account. For this purpose we used thermodynamic integration to calculate the marginal likelihoods of both models, which indicates the expected likelihood to observe the data, averaged over all possible parameterizations. Then we used the Bayes factor K A,B , which is the ratio of these two probabilities, to rank the two model variants. The logarithmic Bayes factor reads 2 logK A,B =12.3±0.4, with standard deviation that was estimated from three repeated sampling runs. Hence, this analysis renders model A superior to model B with very strong evidence, as described in  (see Additional file 1: Sections 4.3, 5 and 6.1).
Summarizing these results, our model study favors model A, in which CERT and PKD regulate CERT transport activity in a cooperative manner. Thus, we decided to focus primarily on the analysis of model A in the following.
Finally, Figure 4E shows an independent validation experiment that was not used for parameter estimation but conducted after the model was used for a prediction of this scenario. Here, PKD activity was inhibited 24 h after ectopically expressing PI4KIII β and phospho-PI4KIII β was measured 60 min after inhibition (see Additional file 1: Figure S1B). It can be seen that also in this scenario model prediction and experiment are in very good agreement. Thus, our model not only reproduces a set of different experiments, but can also be used for reliable predictions of new perturbation scenarios.
Hence we are equipped with an appropriate model that captures absolute protein levels and dynamic responses of the system to different perturbations.
The endogenous TGN: highly responsive regulatory couplings
We next use model A to predict and analyze the network in the endogenous state. Therefore, we simulated the system without any perturbations by setting all inputs to zero. The system reaches a dynamic equilibrium in which concentrations and fluxes are constant and ceramide is constantly extracted from the ER and released to the TGN membrane. We calculated distributions for the concentration of each molecular species, as well as distributions for fluxes and regulatory influences in the network. The result is illustrated in Figure 5. Expectation values for molecule numbers per cell and fluxes in molecules per hour are visualized by different colors on a logarithmic scale. Additionally, distributions for concentrations are shown next to each node along with expectation values and 99.5% confidence intervals.
In accordance with the experimental absolute quantification experiments (Figure 2), PI4KIII β is the most abundant protein with about 1 mio molecules per cell, while PKD and CERT are both present with about half a million molecules. Furthermore, most PKD and PI4KIII β molecules are present in their inactive non-phosphorylated state, while there remains large uncertainty for the active states.
Flux rates correlate positively with respective protein abundances. Accordingly, highest fluxes are encountered for phosphorylation and dephosphorylation of PI4KIII β, followed by lower rates for PKD and CERT. Basal phosphorylation and dephosphorylation rates of PKD are two orders of magnitude higher than ceramide transfer dependent activation. The model also predicts that the majority of CERT molecules is unphosphorylated. However, the rates of phosphorylation and dephosphosphorylation are rather similar, indicating that dephosphorylation at the ER is very efficient. An uncertainty of several orders of magnitude is predicted for TGN bound CERT (CERTaTGN), suggesting that our data agrees with different TGN resting times.
Regulatory influences that represent the feedback in the network are indicated by gray lines. They have been modeled with Michael Menten like kinetics that are upper bounded with maximal rates v max . Regulation strengths are given in terms of these v max values. As can be seen in Figure 5, all regulation strengths are far below their saturation values, indicating that they are highly sensitive to changes in molecular concentrations of the regulators.
Summarizing, we are equipped with a quantitative holistic picture of endogenous TGN functioning including also uncertainties. The overall picture suggests that the system operates in a state that is highly responsive to changes in the regulator concentrations.
Active PKD regulates CERT dependent ceramide transfer
The question arises whether this sensitivity of single edges also results in a highly responsive overall network response of the closed loop feedback system and which of the key players determines the operating state of the endogenous system. We therefore analyzed the mutual influence of the phosphorylation state of PKD and CERT transport activity. Simulation results are shown in Figure 6. Figure 6A shows the relative change in PKD activity in the case that ceramide transport activity of CERT is completely switched off. It can be seen that PKD activity is only slightly affected in this scenario, with an expected increase of less than 1% with a rather small variance. The model preserves PKD activity via the basal ceramide independent phosphorylation reaction, see Figure 5. Similarly, PKD activity shows an increase of about only 0.1% in response to an increase in ceramide transfer activity by 10% (Figure 6B). These results indicate that PKD activity is only marginally influenced by CERT transport activity and that CERT has a rather weak effect on PKD activity. By contrast, a 10% increase in PKD activity resulted in an expected increase in CERT transport activity of about 6% (Figure 6C), leading to the conclusion that PKD is a good candidate to regulate CERT transport activity. In summary, these results identify PKD as an important regulator of the systems operating point, which is in turn not much affected by CERT.
In this study we present for the first time a quantitative dynamic model describing the phosphorylation-dependent CERT-mediated regulation of ceramide transfer between ER and Golgi membranes. Our model comprises the complex interrelated feedback loops of interdependent TGN-localized protein and lipid kinases on the activity of the lipid transfer protein CERT. For model calibration we used statistical sampling-based Bayesian inference techniques, which yield beside point estimates also information about uncertainties in terms of probability distributions. Our computational model was able to reproduce the kinetic response of the molecular interaction network comprising PKD, PI4KIII β and CERT in terms of phosphorylation states to a set of different perturbation experiments (Figures 3 and 4). Furthermore, our model reliably predicted the behavior to new perturbations, which were experimentally verified, validating the predictive power of the model.
Using Bayesian model comparison techniques, we were able to rank model variants that differ in the impact of phosphorylation of CERT at the SR motif regarding its ceramide transfer activity. In model A (Figure 1), which is supported by our study, PKD-mediated phosphorylation of CERT constitutes an overall positive feedback loop between CERT transfer activity and PKD activity. By contrast, in model B, this feedback has a negative effect on CERT-mediated ceramide transfer. Our results thus suggest that CERT and PKD interact in a cooperative rather than a competitive manner to regulate ceramide transfer activity. An experiment that helped to distinguish between the two models was the ectopic expression of PIKIII β, followed by the measurement of the phosphorylation and thus activation state of PKD (Figure 4D). Model A was able to mimic the experimentally observed early increase in PKD phosphorylation, which was not the case for model B. In the latter model this is due to an immediate buffering of this positive effect by the strong negative feedback on CERT mediated ceramide transfer via PKD phosphorylation.
Although energy-consuming, the circular reaction scheme has the advantage that it can rapidly respond to changes at the level of the ER or Golgi to either accelerate or attenuate lipid metabolic pathways required for vesicle formation at the TGN. Indeed, our model predicts a high responsiveness for all regulations in the endogenous system to changes in the regulator concentrations.
Additionally, we predicted protein amounts and their phosphorylation states, as well as fluxes between these states and strengths of regulatory effects within the entire molecular network for the endogenous system (Figure 5). Predicted absolute protein amounts were in good agreement with corresponding measurements (Figure 2). In HEK293T cells, PI4KIII β is the most abundant protein with about 1 million molecules per cell, followed by PKD, while the average amount of CERT molecules per cell is about one order of magnitude lower. Moreover, our model predicts that the majority of both the PKD and PI4KIII β molecules are present in an inactive unphosphorylated state, while most CERT molecules are unphosphorylated and thus transfer-active. Unfortunately, there is currently limited literature about absolute protein abundances in HEK293T pathways, making a direct comparison of our quantification results difficult. However, a recent study of the Wnt pathway in HEK293T cells reports an average abundance of 1.6·105 copies per cell for these signaling molecules . This is in the same order of magnitude as our average abundance of the key proteins for the endogenous system, which we estimated to be 7·105 copies per cell.
Our simulations furthermore suggest that PKD has a strong regulatory effect on CERT transfer activity, while, in turn, PKD activity is only modestly affected by the amount of ceramide that is delivered to the TGN via CERT, rendering PKD a key regulator of the whole network state (Figure 6). This also suggests that at the TGN, the DAG pool produced as a consequence of CERT-mediated lipid transfer plays a minor role in PKD recruitment and/or activation. Therefore, alternative pathways, such as the hydrolysis of PI(4,5)P 2 by PLC β3  or dephosphorylation of phosphatidic acid by lipid-phosphate phosphatases [26,47] most likely provide DAG to locally regulate PKD at the TGN membranes.
One limitation of our current model is its calibration using data derived from ectopic expression experiments. The predictions made for the behavior of the endogenous system thus contain large uncertainties. However, early PKD responses and absolute protein measurements were obtained in minimally perturbed cells. The model is currently calibrated with data from HEK293T cells since this cell line is well suited for experimental interventions like transfections with expression constructs. However, multiple interactions in our model have also been identified in other mammalian epithelial cell lines like HeLa cells, and structural homology in the regulation of the TGN has been observed [4,18].
An important issue in Bayesian learning approaches is the impact of the prior distribution on the results. Since there was limited quantitative literature data that could directly be included into our model, we decided to work with a rather conservative approach. We used log uniform prior distributions for all parameters, which cover a wide range of orders of magnitude. Bounds were taken from the literature whenever available, such as for example in case of protein half lifes or protein amounts [48,49]. In addition, we took initial parameter optimization runs into account to set sufficiently conservative estimates for upper and lower bounds. Thus our priors are rather uninformative. Since the data contains quite some information about the model, as can be seen by the model fits, we believe that our results are not dominated by the prior distributions. However, as it is usually the case with Bayesian inference, we are aware that the results depend on the parametrization and thus might look different when using uniform priors after a transformation of the parameters.
In the future, our model is not only useful to predict perturbations at the protein level, but can also serve as a starting point to investigate further molecular interactions of interest via suitable model extensions. First, by combining the model developed here with a previous model describing the regulation of the SMS1 reaction at the TGN , it will be possible to predict lipid synthesis and transfer rates in vivo, which are difficult to measure experimentally.
Second, experimental results described in  provide ideas for further extensions regarding CERT transport activity. While our model is qualitatively able to reproduce the observed decrease in ceramide transport of a CERT-10E mutant, which mimics the hyperphosphorylated state, capturing the behavior of the S132A mutant that cannot be phosphorylated at the SRM, would require structural extensions. Furthermore, it was reported that the amount of sphingomyelin in the plasma membrane affects CERT activity via the regulation of its phosphorylation state, thereby constituting a negative feedback mechanism in cellular lipid homeostasis. Interestingly, phosphorylation at an additional site, namely serine 315, via a still unknown kinase was recently reported to enhance the affinity of CERT for VAP at the ER, thereby upregulating ER to Golgi trafficking of ceramide . The post-translational modification of CERT by phosphorylation thus appears to be a central general mechanism by which ceramide transfer at MCS is regulated.
Third, alternative metabolic pathways that contribute to local DAG generation at the TGN membranes can be integrated into our model, enabling the study of the dynamic behavior of PKD activation and/or recruitment .
Finally, it will be interesting to integrate into our model additional proteins such as the sterol transfer protein OSBP. Similar to CERT, OSBP contains an FFAT motif and a PH domain which target the protein to ER and Golgi membranes, respectively . Importantly, OSBP is also regulated by PKD-mediated phosphorylation and is thought to stabilize CERT localization at MCS, thereby coordinating sphingomyelin synthesis with sterol metabolism .
In sum, our computational analysis of ceramide transport between the ER and Golgi membranes reveals the cooperativity between PKD and CERT and identifies PKD as the critical regulator in the system. Based on these findings and its further implications, we believe that our quantitative computational model constitutes a valuable contribution towards the holistic mechanistic understanding of the molecular network supporting Golgi secretory function.
Availability of supporting data
The data sets supporting the results of this article are included within additional files.
Membrane contact site
Protein kinase D
- PI4KIII β :
Phosphatidylinositol 4-kinase III beta
Ceramide transfer protein
Vesicle associated membrane associated protein
- PP2C ε :
Protein phosphatase 2C epsilon
- CKI γ2:
Casein kinase I gamma-2
Steroidogenic acute regulatory protein
Two phenylalanine in an acidic tract
Human embryonic kidney cells 293
- PI(4,5)P 2 :
- PLC β3:
Phospholipase C beta 3
Markov Chain Monte Carlo
Mogelsvang S, Marsh BJ, Ladinsky MS, Howell KE. Predicting function from structure: 3D structure studies of the mammalian Golgi complex.Traffic. 2004; 5(5):338–45. doi:10.1111/j.1398-9219.2004.00186.x.
Hanada K, Kumagai K, Yasuda S, Miura Y, Kawano M, Fukasawa M, et al.Molecular machinery for non-vesicular trafficking of ceramide. Nature. 2003; 426(6968):803–9. doi:10.1038/nature02188.
Bard F, Malhotra V. The formation of TGN-to-Plasma-Membrane transport carriers. Annu Rev Cell Dev Biol. 2006; 22(1):439–55. doi:10.1146/annurev.cellbio.21.012704.133126. http://www.annualreviews.org/doi/pdf/10.1146/annurev.cellbio.21.012704.133126.
Malhotra V, Campelo F. PKD regulates membrane fission to generate TGN to cell surface transport carriers. Cold Spring Harb Perspect Biol. 2011; 3(2):a005280. doi:10.1101/cshperspect.a005280.
Olayioye MA, Hausser A. Integration of non-vesicular and vesicular transport processes at the Golgi complex by the PKD-CERT network. Biochim Biophys Acta. 2012; 1821(8):1096–103. doi:10.1016/j.bbalip.2011.12.005.
Hanada K. Intracellular trafficking of ceramide by ceramide transfer protein. Proc Jpn Acad Ser B Phys Biol Sci. 2010; 86:426–37.
Hausser A, Storz P, Märtens S, Link G, Toker A, Pfizenmaier K. Protein kinase D regulates vesicular transport by phosphorylating and activating phosphatidylinositol-4 kinase IIIbeta at the Golgi complex. Nat Cell Biol. 2005; 7(9):880–6. doi:10.1038/ncb1289.
Tóth B, Balla A, Ma H, Knight ZA, Shokat KM, Balla T. Phosphatidylinositol 4-kinase III β regulates the transport of ceramide between the endoplasmic reticulum and Golgi. J Biol Chem. 2006; 281(47):36369–77. doi:10.1074/jbc.M604935200.
D’Angelo G, Vicinanza M, Campli AD, Matteis MAD. The multiple roles of PtdIns(4)P – not just the precursor of PtdIns(4,5)P2. J Cell Sci. 2008; 121(Pt 12):1955–63. doi:10.1242/jcs.023630.
Hanada K. Discovery of the molecular machinery CERT for endoplasmic reticulum-to-Golgi trafficking of ceramide. Mol Cell Biochem. 2006; 286(1–2):23–31. doi:10.1007/s11010-005-9044-z.
Hanada K, Kumagai K, Tomishige N, Kawano M. CERT and intracellular trafficking of ceramide. Biochim Biophys Acta. 2007; 1771(6):644–53. doi:10.1016/j.bbalip.2007.01.009.
Hanada K, Kumagai K, Tomishige N, Yamaji T. CERT-mediated trafficking of ceramide. Biochim Biophys Acta. 2009; 1791(7):684–91. doi:10.1016/j.bbalip.2009.01.006.
Perry RJ, Ridgway ND. Molecular mechanisms and regulation of ceramide transport. Biochim Biophys Acta. 2005; 1734(3):220–34. doi:10.1016/j.bbalip.2005.04.001.
Huitema K, van den Dikkenberg J, Brouwers JFHM, Holthuis JCM. Identification of a family of animal sphingomyelin synthases. EMBO J. 2004; 23(1):33–44. doi:10.1038/sj.emboj.7600034.
Tafesse FG, Ternes P, Holthuis JCM. The multigenic sphingomyelin synthase family. J Biol Chem. 2006; 281(40):29421–25. doi:10.1074/jbc.R600021200.
Tafesse FG, Huitema K, Hermansson M, van der Poel S, van den Dikkenberg J, Uphoff A, et al.Both sphingomyelin synthases SMS1 and SMS2 are required for sphingomyelin homeostasis and growth in human HeLa cells. J Biol Chem. 2007; 282(24):17537–47. doi:10.1074/jbc.M702423200.
Maeda Y, Beznoussenko GV, Lint JV, Mironov AA, Malhotra V. Recruitment of protein kinase D to the trans-Golgi network via the first cysteine-rich domain. EMBO J. 2001; 20(21):5982–90. doi:10.1093/emboj/20.21.5982.
Baron CL, Malhotra V. Role of diacylglycerol in PKD recruitment to the TGN and protein transport to the plasma membrane. Science. 2002; 295(5553):325–8. doi:10.1126/science.1066759.
Añel AMD, Malhotra V. PKCeta is required for β1γ2/ β3γ2- and PKD-mediated transport to the cell surface and the organization of the Golgi apparatus. J Cell Biol. 2005; 169(1):83–91. doi:10.1083/jcb.200412089.
Fugmann T, Hausser A, Schöffler P, Schmid S, Pfizenmaier K, Olayioye MA. Regulation of secretory transport by protein kinase D-mediated phosphorylation of the Ceramide Transfer Protein. J Cell Biol. 2007; 178:15–22.
Tomishige N, Kumagai K, Kusuda J, Nishijima M, Hanada K. Casein kinase I gamma 2 down-regulates trafficking of ceramide in the synthesis of sphingomyelin. Mol Biol Cell. 2009; 20(1):348–57. doi:10.1091/mbc.E08-07-0669.
Kumagai K, Kawano M, Shinkai-Ouchi F, Nishijima M, Hanada K. Interorganelle trafficking of ceramide is regulated by phosphorylation-dependent cooperativity between the PH and START domains of CERT. J Biol Chem. 2007; 282(24):17758–66. doi:10.1074/jbc.M702291200.
Saito S, Matsui H, Kawano M, Kumagai K, Tomishige N, Hanada K, et al.Protein phosphatase 2C ε is an endoplasmic reticulum integral membrane protein that dephosphorylates the ceramide transport protein CERT to enhance its association with organelle membranes. J Biol Chem. 2008; 283(10):6584–93. doi:10.1074/jbc.M707691200.
Hirschberg K, Miller CM, Ellenberg J, Presley JF, Siggia ED, Phair RD, et al.Kinetic analysis of secretory protein traffic and characterization of Golgi to plasma membrane transport intermediates in living cells. J Cell Biol. 1998; 143(6):1485–503.
Corda D, Carcedo CH, Bonazzi M, Luini A, Spanò S. Molecular aspects of membrane fission in the secretory pathway. Cell Mol Life Sci. 2002; 59(11):1819–32.
Shemesh T, Luini A, Malhotra V, Burger KNJ, Kozlov MM. Prefission constriction of Golgi tubular carriers driven by local lipid metabolism: a theoretical model. Biophys J. 2003; 85(6):3813–27. doi:10.1016/S0006-3495(03)74796-1.
Thomaseth C, Weber P, Hamm T, Kashima K, Radde N. Modeling sphingomyelin synthase 1 driven reaction at the Golgi apparatus can explain data by inclusion of a positive feedback mechanism. J Theor Biol. 2013; 337:174–80. doi:10.1016/j.jtbi.2013.08.022.
Florin L, Pegel A, Becker E, Hausser A, Olayioye MA, Kaufmann H. Heterologous expression of the lipid transfer protein CERT increases therapeutic protein productivity of mammalian cells. J Biotechnol. 2009; 141:84–90. doi:10.1016/j.jbiotec.2009.02.014.
Gelman A, Carlin JB, Stern HS, Rubin DB. Bayesian Data Analysis, 2nd edn. Texts in Statistical Science. Boca Raton, London, NewYork, Washington D.C: Chapman & Hall, CRC; 2004.
Gutenkunst RN, Waterfall JJ, Casey FP, Brown KS, Myers CR, Sethna JP. Universally sloppy parameter sensitivities in systems biology models. PLoS Comp Biol. 2007; 3(10):1871–78.
Eydgahi H, Chen WW, Muhlich JL, Vitkup D, Tsitsiklis JN, Sorger PK. Properties of cell death models calibrated and compared using Bayesian approaches. Mol Syst Biol. 2013; 9(644):644.
Hausser A, Link G, Bamberg L, Burzlaff A, Lutz S, Pfizenmaier K, et al.Structural requirements for localization and activation of protein kinase C mu (PKC mu) at the Golgi compartment. J Cell Biol. 2002; 156(1):65–74. doi:10.1083/jcb.200110047.
Huck B, Kemkemer R, Franz-Wachtel M, Macek B, Hausser A, Olayioye MA. GIT1 phosphorylation on serine 46 by PKD3 regulates paxillin trafficking and cellular protrusive activity. J Biol Chem. 2012; 287(41):34604–13. doi:10.1074/jbc.M112.374652.
Gassmann M, Grenacher B, Rohde B, Vogel J. Quantifying western blots: pitfalls of densitometry. Electrophoresis. 2009; 30(11):1845–55. doi:10.1002/elps.200800720.
Schmidt H, Jirstrand M. Systems biology toolbox for matlab: a computational platform for research in systems biology. Bioinformatics. 2006; 22(4):514–5. doi:10.1093/bioinformatics/bti799. http://bioinformatics.oxfordjournals.org/cgi/reprint/22/4/514.pdf.
Cohen S, Hindmarsh C. Cvode, A Stiff/nonstiff Ode Solver In C. Comput Phys. 1996; 10(2):138–43.
Calderhead B, Girolami M. Estimating bayes factors via thermodynamic integration and population mcmc. Comput Stat Data Anal. 2009; 53:4028–45.
Gelman A, Rubin DB. Inference from iterative simulation using multiple sequences. Stat Sci. 1992; 7(4):457–72.
Kawano M, Kumagai K, Nishijima M, Hanada K. Efficient trafficking of ceramide from the endoplasmic reticulum to the Golgi apparatus requires a vamp-associated protein-interacting ffat motif of cert. J Biol Chem. 2006; 281(40):30279–88. doi:10.1074/jbc.M605032200.
Kreutz C, Timmer J. Systems biology: experimental design. FEBS J. 2009; 276(4):923–42. doi:10.1111/j.1742-4658.2008.06843.x.
Chiu T, Rozengurt E. PKD in intestinal epithelial cells: rapid activation by phorbol esters, LPA, and angiotensin through PKC. Am J Physiol Cell Physiol. 2001; 280(4):929–42.
Bravo-Altamirano K, George KM, Frantz M-C, Lavalle CR, Tandon M, Leimgruber S, et al.Synthesis and structure-activity relationships of benzothienothiazepinone inhibitors of protein kinase D. ACS Med Chem Lett. 2011; 2(2):154–9. doi:10.1021/ml100230n.
Weber P, Kramer A, Dingler C, Radde N. Trajectory-oriented bayesian experiment design versus Fisher A-optimal design: an in depth comparison study. Bioinformatics. 2012; 28(18):535–41. doi:10.1093/bioinformatics/bts377. http://bioinformatics.oxfordjournals.org/content/28/18/i535.full.pdf+html.
Kass RE, Raftery AE. Bayes factors. J Am Stat Assoc. 1995; 90(430):773–95.
Tan CW, Gardiner BS, Hirokawa Y, Layton MJ, Smith DW, Burgess AW. Wnt signalling pathway parameters for mammalian cells. PLoS One. 2012; 7(2):31882. doi:10.1371/journal.pone.0031882.
Díaz Añel AM. Phospholipase c beta3 is a key component in the gbetagamma/pkceta/pkd-mediated regulation of trans-Golgi network to plasma membrane transport. Biochem J. 2007; 406(1):157–65. doi:10.1042/BJ20070359.
Sarri E, Sicart A, Lázaro-Diéguez F, Egea G. Phospholipid synthesis participates in the regulation of diacylglycerol required for membrane trafficking at the Golgi complex. J Biol Chem. 2011; 286(32):28632–43. doi:10.1074/jbc.M111.267534.
Eden E, Geva-Zatorsky N, Issaeva I, Cohen A, Dekel E, Danon T, et al.Proteome half-life dynamics in living human cells. Science. 2011; 331(6018):764–8. doi:10.1126/science.1199784.
Finka A, Goloubinoff P. Proteomic data from human cell cultures refine mechanisms of chaperone-mediated protein homeostasis. Cell Stress Chaperones. 2013; 18(5):591–605. doi:10.1007/s12192-013-0413-3.
Kumagai K, Kawano-Kawada M, Hanada K. Phosphoregulation of the ceramide transport protein cert at serine 315 in the interaction with vamp-associated protein (vap) for inter-organelle trafficking of ceramide in mammalian cells. J Biol Chem. 2014; 289(15):10748–60. doi:10.1074/jbc.M113.528380.
Lev S. Lipid homoeostasis and Golgi secretory function. Biochem Soc Trans. 2006; 34(Pt 3):363–6. doi:10.1042/BST0340363.
Nhek S, Ngo M, Yang X, Ng MM, Field SJ, Asara JM, et al.Regulation of osbp Golgi localization through protein kinase d-mediated phosphorylation. Mol Biol Cell. 2010; 21:2327–2337. doi:10.1091/mbc.E10-02-0090.
We want to thank Thomas Hamm for important preliminary literature research throughout his diploma thesis. We also want to thank Gisela Link for important feedback to our experimental procedures. All authors acknowledge financial support from the German Research Foundation (DFG) (GZ:RA 1840/1-1). Furthermore, the authors NR and PW would like to thank the German Research Foundation (DFG) for financial support of the project within the Cluster of Excellence in Simulation Technology (EXC 310/2) at the University of Stuttgart.
The authors declare that they have no competing interests.
PW, MO, AH, NR developed the idea of the holistic model study and wrote the article. All authors planned and refined the wet lab experiments in regular lab meetings. PW and NR established the theoretical foundation. PW did the modeling and performed the simulation study. MH performed the wet lab experiments. All authors read and approved the final manuscript.
Supplementary text. Description: A pdf-file including a detailed theoretical introduction to ODE models, Bayesian Inference and the employed MCMC sampling algorithm. It includes tables of the model equations, further supporting experimental results and explanations for data normalization and error calculus.
Model files and parameters. Software needed: MATLAB R2011b (or higher), SBtoolbox2, SBML viewer, text editor. A compressed archive containing the simulation models A and B (SBML and SBtoolbox2 formats). It includes general and SBtoolbox2 specific simulation instructions (readme.txt, experiment.exp) and two MATLAB Data-files containing double precision variables of 5000 representative parameter samples.
Experimental data. Software needed: spreadsheet viewer (MS Excel, Libre Office, Open Office), Adobe Reader. A compressed archive containing a xls-sheet with readily quantified image data and two pdf-files with a summary of the raw image data in form of expert labeled Western Blots.