Mechanistic models of PLC/PKC signaling implicate phosphatidic acid as a key amplifier of chemotactic gradient sensing.

Chemotaxis of fibroblasts and other mesenchymal cells is critical for embryonic development and wound healing. Fibroblast chemotaxis directed by a gradient of platelet-derived growth factor (PDGF) requires signaling through the phospholipase C (PLC)/protein kinase C (PKC) pathway. Diacylglycerol (DAG), the lipid product of PLC that activates conventional PKCs, is focally enriched at the up-gradient leading edge of fibroblasts responding to a shallow gradient of PDGF, signifying polarization. To explain the underlying mechanisms, we formulated reaction-diffusion models including as many as three putative feedback loops based on known biochemistry. These include the previously analyzed mechanism of substrate-buffering by myristoylated alanine-rich C kinase substrate (MARCKS) and two newly considered feedback loops involving the lipid, phosphatidic acid (PA). DAG kinases and phospholipase D, the enzymes that produce PA, are identified as key regulators in the models. Paradoxically, increasing DAG kinase activity can enhance the robustness of DAG/active PKC polarization with respect to chemoattractant concentration while decreasing their whole-cell levels. Finally, in simulations of wound invasion, efficient collective migration is achieved with thresholds for chemotaxis matching those of polarization in the reaction-diffusion models. This multi-scale modeling framework offers testable predictions to guide further study of signal transduction and cell behavior that affect mesenchymal chemotaxis.

Introduction

Chemotaxis, the bias of cell movement towards soluble chemical cues (chemoattractants), is critical for embryonic development, angiogenesis, the immune response, and wound healing in metazoans [1]. Fibroblasts, the cells directly responsible for regenerating wounded tissue, respond chemotactically to platelet-derived growth factor (PDGF) as a cue to invade the wound; this proliferative phase of wound healing typically unfolds over the course of several days [2-5]. Chemoattractant ligands such as PDGF bind to cognate receptors on the cell surface, and in eukaryotic cells this signal is sensed spatially, relying on a gradient of receptor occupancy and activation. The activated receptors interface with a network of intracellular signaling pathways that modulate the dynamics of the cytoskeleton and thus cell motility. This modulation can be achieved by spatially varying the rate of F-actin polymerization or the mechanical influence of Myosin II contractility [6,7]. In fibroblasts exposed to a steady PDGF gradient, many of the prominent signaling pathways that enhance the rate of F-actin polymerization have been found to be dispensable for chemotaxis [8-10], whereas regulation of Myosin IIA by phospholipase C (PLC)/protein kinase (PKC) signaling, a well-studied pathway activated by many receptors, is essential [11]. Another key finding in that study was that diacylglycerol (DAG), the lipid product of PLC that activates most PKC isoforms, is sharply concentrated in the fibroblasts’ protrusions (lamellipodia) exposed to the highest concentration of PDGF [11]. Given that the external gradients in such experiments are characteristically shallow (typically, ~ 5% across a cell’s length), the PLC/PKC signaling circuit must be locally amplified to explain the observed polarization of DAG production. What are the biochemical and biophysical mechanisms that cause the pathway to polarize?

To partially address this question, a reaction-diffusion model of the PLC/PKC signaling pathway identified phosphorylation of myristoylated alanine-rich C kinase substrate (MARCKS) by PKC, which increases the availability of the PLC substrate (PIP2) [12-15], as a positive feedback loop (PFL); this mechanism is sufficient for polarization in response to an abnormally steep external gradient [16]. By itself, however, the MARCKS feedback was unable to polarize signaling in response to shallow external gradients (≤ 10%), and the system lacked robustness to modest changes in the midpoint concentration of chemoattractant. In this work, we address these issues through formulation of more mechanistic, partial differential equation models of the PLC/PKC pathway. These models consider two additional PFLs supported by literature evidence and thus introduce a key molecular player: phosphatidic acid (PA), a lipid intermediate in the metabolism of DAG.

PA is recognized as a signaling molecule affecting a number of cellular functions including cell growth and proliferation, vesicular trafficking, and cytoskeletal rearrangement [17,18]. In the present models, we include reactions by which PA is produced by phosphorylation of DAG by DAG kinases or from hydrolysis of phosphatidylcholine by phospholipase D (PLD) [19-22]. Feedback loops incorporating PA were added to the model based on published evidence that: 1) PA binds PLCγ and increases the rate of PIP2 hydrolysis in vitro [23]; and 2) active PKC can enhance the activity of PLD for increased production of PA [24-26]. Model simulations show that the MARCKS feedback mechanism synergizes with these new feedback loops to polarize PLC/PKC signaling in response to shallow gradients of receptor occupancy and over an appreciable range of midpoint occupancy. Subtle asymmetry of the cell geometry can also polarize signaling. Focusing on the molecular details, simulations suggest that DAG kinases, the enzymes responsible for turnover of DAG, exert a critical and surprisingly positive influence on the responsiveness of the circuit. Finally, we applied the predicted receptor occupancy conditions for polarization to a two-state model of cell movement in a hybrid (stochastic/continuum) model of wound invasion. This model predicts a hierarchy of chemotactic waves that drive more efficient collective invasion than a single chemotactic front. With this framework, we can proceed to link signal transduction mechanisms at the molecular level to individual and collective cell movements directed by chemoattractant gradients in tissues.

Results

New models of the PLC/PKC pathway based on putative feedback mechanisms indicated in the literature

We formulated a biochemically realistic description of the chemotactic sensing circuit in fibroblasts, retaining the differential PIP2 buffering by MARCKS from the Mohan model [16] and considering two additional PFLs based on published evidence (; see also and ). The first, PFL 1, considers the modulation of PLC recruitment by phosphatidic acid (PA), the lipid produced by phosphorylation of DAG by DAG kinases. Using a detergent-phospholipid mixed micelle assay system, it was shown that inclusion of PA enhanced PIP2 hydrolysis catalyzed by either unphosphorylated or tyrosine-phosphorylated PLCγ1, by reducing the apparent K for the reaction [23]. The details of how PA affects PLC activity are not completely understood, but the reduction of the apparent K is consistent with PA-mediated stabilization of PLCγ1 association with the membrane, akin to the effect of the non-catalytic interaction of PLCδ with PIP2 [27]. Therefore, we modeled the effect of PA as an increased lifetime of the receptor-PLCγ1 complex at the membrane ().

PFL 2 considers the effect of active PKC on the activity of phospholipase D1 (PLD1), which produces PA by hydrolyzing the abundant lipid, phosphatidylcholine [19,21]. In murine fibroblasts stimulated with phorbol ester, PKCα interacts with PLD1 and increases the rate of phosphatidylcholine hydrolysis. PA produced by this reaction can be dephosphorylated to yield DAG, and thus PFL 2 exerts an influence on PLC/PKC signaling independent of PFL 1. We model the influence of PLD as a distinct source term for generation of PA, which increases according to a Hill function of the active PKC density at the membrane ().

To simplify and generalize the handling of receptor dynamics in this model, and considering that receptor activation is patterned by an external ligand, we assume a steady gradient of occupied/active receptors, r, that is linear in the direction of the cell’s long axis (on both the top and bottom of the cell) rather than along the cell’s contour (). In the expression for r, rfrac is the average receptor occupancy, expressed as a fraction of a characteristic receptor density of 130 μm-2, or 105/cell. The parameter rsteep is the relative steepness of the receptor occupancy gradient across the cell; for example, a value of rsteep = 0.1 corresponds to a 10% difference between the front and back of the cell.

Stabilization of PLC recruitment by phosphatidic acid (PFL 1), combined with neutralization of MARCKS by PKC, promotes sensitive and robust gradient sensing

Considering the network depicted in , with PFL 1 but not PFL 2, we evaluated the ability of receptor occupancy gradients (characterized by midpoint occupancy and % steepness) to polarize DAG and active PKC. For each of five gradient steepness values, ranging from 0% (uniform stimulation) up to 67% (2-fold) difference across the cell, simulations were run varying the value of the midpoint receptor occupancy, rfrac. All simulations were run sufficiently long to allow a steady state to be achieved, followed by an equally long period with the gradient reversed to check the stability of the spatial pattern. The steady-state concentrations of active PKC at the front and back of the cell are plotted versus rfrac (). In some simulations, the spatial pattern oscillated, in which case the maximum and minimum values of the oscillation at each end of the cell are plotted. For examples of the simulated time courses, showing the transient behavior with sustained oscillations where applicable, see The PKC activity pattern was reversible for all but one of these simulated conditions. The exception is marked with an asterisk in , signifying lack of reversibility here and in subsequent figures of this paper.

Fig 1

Model of the PLC/PKC network including phosphatidic acid (PA).

(A) Model schematic depicting the interactions and reactions among signaling proteins and plasma membrane lipids. Dashed lines ending in a filled circle indicate that the species enhances the associated process. The reactions and interactions shown in red are associated with the generation and influence of PA in positive feedback loops (PFLs) labeled (1) and (2). (B) Base model geometry and orientation of the receptor occupancy gradient, which is linear in the direction of the cell’s long axis, x. (C) Plots illustrating the linear profile of active receptor density imposed across the 40 μm length of the cell for varying values of rfrac (relative midpoint density) and rsteep (relative steepness, expressed here as a percentage difference across the cell).

Gradient amplification by PFL 1 combined with regulation of MARCKS.

(A) The concentration of active PKC molecules at the front (red circle) and back (blue circle) of the cell are plotted as a function of the mean fractional occupancy of receptors, rfrac, for varying values of gradient steepness. When the simulations produced oscillations, the maxima (upward-pointing triangles) and minima (downward-pointing triangles) of the oscillations are plotted; the front and back are still denoted by shades of red and blue, respectively. The simulation achieving steady state with the maximum front/back ratio is denoted by the dashed vertical line. For these simulations, the direction of the gradient was reversed after 20,000 s. If the active PKC pattern failed to reverse in response, the simulation is marked with an asterisk. (B) Table showing the maximum front/back ratio for each gradient steepness in simulations run with PFL1 and either with or without MARCKS protein.

For simulations that achieved steady state, the degree of polarization is quantified as the ratio of active PKC at the front of the cell over that of the back of the cell, and the value of rfrac where this ratio is maximum is marked on each of the plots in . For receptor occupancy gradients as shallow as 3%, we found maximum front/back ratios exceeding 10, indicating strong and sensitive polarization of the PLC/PKC pathway. Another quality of the gradient sensing response is its robustness, of which several features might be considered. For example, across the range of % steepness values, the greatest degree of polarization is consistently achieved at approximately the same value of rfrac (≈ 0.1), and polarized patterns exhibit a consistent density of active PKC at the front of the cell. The aspect of greatest interest to us, however, is the range of rfrac values that elicit polarization for a given % steepness, which we refer to as dose-response robustness. As might be expected, dose-response robustness is enhanced as gradient steepness is increased (); at a modest steepness of 10% (typical in magnitude of chemotaxis experiments), the range of rfrac values for which steady polarization would be evident spans a factor of 2. The range doubles if oscillatory simulations with consistently high front/back ratio are included.

It should be noted that the general inability to polarize the gradient sensing circuit at high rfrac is not caused by saturation of receptor occupancy, as would be expected at high concentrations of chemoattractant; in this model, the relative steepness of receptor activation is fixed. Rather, high receptor activation promotes ignition of the positive feedback at the back of the cell as well as at the front.

In previous work [16], it was shown that a single positive feedback was not sufficient for a demonstratively amplified response, i.e., without differential buffering of PIP2 by MARCKS. Accordingly, simulations with PFL 1 but no MARCKS show similar dose responsiveness with respect to rfrac but inconsequential polarization (). These results identify a promising gradient sensing circuit that combines differential buffering of PIP2 by MARCKS, and the amplification of PLC recruitment by PA, regulation mechanisms supported by evidence in the literature.

Analysis of the mechanisms driving PLC/PKC polarization

The differential buffering of PIP2 by MARCKS was described in detail previously [16], and it is important to understand this concept in the context of the present model as well. Phosphorylation of MARCKS by active, DAG-bound PKC liberates PIP2 locally, further enhancing DAG production by positive feedback. A key aspect of differential buffering is the maintenance of certain intracellular gradients. In the cytosol, PLC and PKC are close to uniform, whereas there is a gradient of phosphorylated versus unphosphorylated MARCKS in the cytosol. Accordingly, a substantial increase in the diffusivities of cytosolic MARCKS species breaks the polarization of the present model (). A large decrease in the diffusivities also decreases polarization, but modestly so. The optimum with respect to diffusivity has been attributed to the ability of dephosphorylated MARCKS to diffuse to the rear of the cell before re-binding PIP2 [16]. The other important gradients are of total and free PIP2 in the plasma membrane. MARCKS-bound PIP2 is protected from hydrolysis, and therefore its spatial range by diffusion is substantial. Thus, net diffusion of MARCKS-bound PIP2 from back to front supplies the front of the cell with substrate, allowing the free PIP2 density to be much higher at the front despite the much higher PLC activity there. Accordingly, reducing the diffusivities of total PIP2 and membrane-associated MARCKS prevents polarization, whereas increasing those diffusivities enhances the extent of polarization ().

Analysis of the mechanisms driving PLC/PKC polarization.

(A) Sensitivity of the results in , 10% gradient steepness, to the indicated fold-changes in MARCKS (M and M) cytosolic diffusivities (1 = base case). (B) Sensitivity of the results in , 10% gradient steepness, to the indicated fold-changes in total PIP2 and membrane-associated MARCKS (p and m) membrane diffusivities (1 = base case). (C) Plot of the d- and e-nullclines () evaluated at the front and back of the cell for the base-case parameters as in , with 10% gradient steepness and rfrac = 0.1.

To further explain the polarization of the present gradient sensing circuit, we developed a steady-state analysis (), with the simplifying assumption that diffusion of active PLC and DAG-containing species are negligible relative to associated reaction terms. We derived nullcline expressions for the membrane-recruited PLC and for the sum of free and PKC-bound DAG, which we refer to as the e- and d-nullclines. The intersections of these two curves are fixed points. To construct these curves on a (d, e) phase plane, certain other variables must be specified. The e-nullcline includes the local activated receptor density, r, and the cytosolic PLC concentration, E. The former is the input to the simulation, with different values at the front and back of the cell, whereas the latter is close to spatially uniform and determined from the simulation. The d-nullcline is a straight line with a slope that is inversely proportional to the free PIP2 concentration, with different values at the front and back of the cell determined from the simulation. The phase-plane plot for the case of 10% gradient, rfrac = 0.1 from illustrates how PFL 1 synergizes with the regulation of MARCKS (). PFL 1 is directly responsible for the positive slope of the e-nullcline; without PFL 1, the e-nullcline has zero slope. In conjunction, regulation of MARCKS is directly responsible for the PIP2 density being higher (lower slope of the d-nullcline) at the front versus the back of the cell; without MARCKS, free PIP2 can only be lower, not higher, at the front [16]. This analysis also shows how the situation changes for lower or higher rfrac, lending insight into dose-response robustness ().

Subtle asymmetry in cell morphology polarizes the PLC/PKC network with PFL 1 and influences external gradient sensing

Given the demonstrated sensitivity of the present gradient sensing circuit and its dependence on diffusion of MARCKS in the cytosol, we reasoned that a slight asymmetry in cell morphology might be a sufficiently strong spatial cue to polarize signaling in response to uniform receptor occupancy. To test this, the cell geometry was altered such that the back end of the cell was blunted compared to the front end (). With uniform stimulation and the same parameter set used in , the shape asymmetry polarized the signaling network, disfavoring the blunted end, for a range of rfrac values (). As expected, this polarization requires a gradient of phosphorylated MARCKS, with a dependence on cytosolic MARCKS diffusivities similar to polarization induced by an external gradient (, compare to ).

Subtle asymmetry in cell morphology causes spontaneous polarization of the proposed PLC/PKC network and influences external gradient sensing.

(A) Comparison of the symmetric cell geometry versus the asymmetric cell geometry. (B) Results for the asymmetric geometry with uniform receptor occupancy (0% gradient steepness) show spontaneous polarization of the system for a particular range of receptor occupancy values. Symbols have the same meanings as in . (C) Sensitivity of the results in B to the indicated fold-changes in MARCKS (M and M) cytosolic diffusivities (1 = base case). (D) Table showing the maximum front/back active PKC ratio for each gradient steepness indicated using the asymmetric geometry. Consistent with the rest of the paper, the front and back of the cell refer to the ends with the highest and lowest receptor activation, respectively. For the front-end blunted simulations where the ratio is less than 1 (internal gradient opposite the external gradient), the inverse of the ratio is indicated.

The tendency of the asymmetric cell geometry to polarize was also tested with receptor occupancy gradients in either direction (). Consistent with the rest of the paper, the front of the cell refers here to the right end, with the highest receptor occupancy, and the back of the cell is at the left end, with the lowest receptor occupancy. With the blunt end of the cell at the back, the maximum degree of active PKC polarization is only modestly affected by receptor occupancy gradients of 3% and 10%, relative to uniform stimulation, whereas 30% and 67% gradients elicit substantially greater polarization similar to the symmetric geometry (compare and ). Consistent with that trend, when the blunt end of the cell was at the front (right-to-left gradient), the shape effect opposes the influence of the receptor occupancy gradient and dominates for 3% and 10% gradients, whereas the steeper gradients were sufficient to overcome the geometry effect. These results indicate that the cell’s local morphology can readily alter or play a dominant role in the polarization of PLC/PKC signaling.

The model predicts a critical role of DAG kinases affecting the responsiveness of the gradient sensing network through PFL 1

Having characterized and explained how PFL 1 influences polarization of the signaling network, we considered the robustness of this system to substantial changes in rate parameters. Each of the rate constants and affinity parameters was increased by 3X and reduced to 0.3X to yield a full order-of-magnitude range, and the analysis shown in was repeated for each. A plot of the maximum front/back ratios of activated PKC shows that polarization is fairly robust to a 3-fold change of each parameter in at least one direction (). The changes that break the polarization are those that weaken PLC recruitment, alter MARCKS phosphorylation by PKC (as explored in [16]), or alter lipid metabolism. The three parameters associated with the latter (k, k, and k) are pseudo-first-order rate constants that reflect enzymatic metabolism of DAG and PA. Reasoning that the interconversion of DAG and PA is a pivotal aspect that needs to be understood, we devised a simple steady-state analysis assuming negligible diffusion of DAG and PA (). This analysis predicts a consistent proportional relationship between the concentrations of DAG and PA. To verify this conclusion, the steady-state concentration of PA at the front of the cell was plotted versus that of DAG for all of the simulations analyzed in ; each of these points lies approximately on the line predicted by the analytical expression ().

Variation of PA/DAG ratio modulates the sensitivity and dose response of the pathway.

(A) The final steady-state concentration of PA (d) is plotted versus the concentration of DAG (d) at the front of the cell for each simulation presented in . These points fall along the line predicted by the equation shown, for which net diffusion of lipids is assumed to be slow. (B) Relative to the parameter set associated with , the parameter k was taken at 0.3 or 3 times its base value, at 10% gradient steepness. Colors and symbols have the same meanings as in . The maximum front/back ratio is shown underneath the symbols. (C) Bar plots comparing the DAG and PA steady-state concentrations at the front and back of the cell for 1x and 3x k simulations, with 10% steepness and rfrac = 0.1. (D) Same as B, but with 3% gradient steepness. For B-D, the direction of the gradient was reversed after 20,000 s. If the active PKC pattern failed to reverse in response, the simulation is marked with an asterisk. (E) The parameter scaling factor refers to the fold-change by which the parameters k, k, k, and K were decreased, while the parameter k was divided by the scaling factor to increase its value. These parameter changes increase DAG and PA levels at the front of the cell systematically, while maintaining approximately the same PA/DAG ratio and comparable effects of PA and DAG on PFL 1 and PKC recruitment, respectively.

The analysis shows that the activity of DAG kinases, reflected in the value of the rate constant, k, directly alters the PA/DAG ratio. Intuitively, diminishing DAG kinase activity is expected to result in higher DAG and active PKC levels; however, in the model with PFL 1 and a nominal 10% gradient of receptor occupancy, a reduction of k to 0.3 times its base value ablates polarization, whereas an increase of k to 3 times its base value substantially enhances the degree of polarization and dose-response robustness (). Analysis of DAG and PA levels for 1x and 3x k () shows that as the parameter is increased, the density of DAG at the cell front is decreased as expected, but only modestly so (~30% reduction) because of the substantially increased abundance of PA; at the rear of the cell, the impact on DAG is relatively far greater (~70% reduction). These effects of increased k extend to receptor-activation gradients with 3% steepness ().

To further test this analysis, we predicted that we could increase the steady-state levels of both DAG and PA, while maintaining the same PA/DAG ratio, by proportional reduction of all three of k, k, and k. To normalize for these changes at the level of PFL 1 and PKC activity, we applied the same reduction factor to the PA-PLC affinity parameter, K, and to the DAG-PKC affinity (by increasing k). With these changes, we were able to maintain polarization while increasing DAG and PA levels by almost an order of magnitude ().

PKC-mediated activation of PLD (PFL 2) confers less responsiveness than PFL 1, but the two feedbacks can synergize in the polarization of PLC/PKC signaling

The other putative source of feedback we consider is the ability of active PKC to amplify the hydrolysis of phosphatidylcholine by PLD, producing PA and additional DAG upon dephosphorylation of PA (PFL 2). Parameters characterizing PFL 2 are a saturation constant, K, gain parameter, γ, and Hill coefficient, n. Extensive parameter sweeps of K and γ were performed with n fixed at either 1 (Michaelean sensitivity) or 2 (modestly ultrasensitive). At a 10% gradient in receptor activation, PFL 2 along with the regulation of MARCKS was unable to polarize DAG and active PKC with n = 1, whereas with n = 2, a particular combination of K and γ values yielded strong polarization; however, the degree of polarization and dose-response robustness are less than those with PFL 1 instead (). That said, ultrasensitive PFL 2 shows synergy when combined with PFL 1. Firstly, the inclusion of both PFLs (with re-optimization of K and γ for PFL 2) allowed for polarization without MARCKS (). When all three feedbacks (PFLs 1 and 2 and regulation of MARCKS) were present, both the degree of polarization and the dose-response robustness were dramatically enhanced, at the expense of reversibility (). Consistent with these conclusions, polarization of the system with PFL 2 is broken by 3-fold changes in parameters that influence PA levels (k or K) except when all three feedbacks are included ().

PFL 2 confers less responsiveness than PFL 1, but the two can synergize in the polarization of PLC/PKC signaling.

The simulations shown are all at a 10% gradient of receptor occupancy. All colors, symbols, and values have the same meanings as in previous figures. (A) PFL 2 and MARCKS regulation but no PFL 1. For PFL 2, the parameters K and γV are set at 1 and 1, respectively. (B) PFL 1 and PFL 2 but no MARCKS. The parameters K and γV are set at 0.1 and 10, respectively. (C) All three feedbacks are included, and K and γV are set as in B. For all of these simulations, the direction of the gradient was reversed after 20,000 s. If the active PKC pattern failed to reverse in response, the simulation is marked with an asterisk.

Thresholds for chemotactic migration matching those of DAG/PKC polarization yield efficient collective invasion and chemotactic wavelets in simulations of wound invasion

Having characterized the polarization of the proposed PLC/PKC network, we asked how it might influence directed cell migration in a physiological context. To address this, we adapted a hybrid simulation of wound invasion, in which the concentration of PDGF evolves as a continuum in one dimension, and fibroblasts are treated as motile line segments [28]. A key prediction of the previous wound invasion model, which did not consider a polarizable gradient sensing mechanism, is that the cells can collectively generate a PDGF gradient through receptor-mediated endocytosis and lysosomal degradation [29], a concept that has since been verified in experimental systems [30-33]. Thus, fibroblast invasion of the simulated wound is guided by a chemotactic wave. This model was modified such that chemotaxis is switched on or off according to the receptor activation conditions that yield polarization of PKC activity for the modestly robust model and base-case parameters used to generate the results shown in . When chemotaxis is switched off, the cell moves in a random direction (left or right). The cell migration parameters were chosen to be consistent with experimentally measured migration speed and chemotactic (forward migration) index values [10,11].

The movements of the cells in the simulation may be animated to visualize the collective behavior ( and ), and the progress of the cell population shows substantially faster invasion relative to simulations in which chemotaxis was not allowed (random migration only) (). This was anticipated because of the chemotactic behavior, but inspection of the simulations revealed an emergent property of the system. While the major chemotactic wave at the leading front of the cell population, located within the steepest gradient in PDGF concentration, remains a primary feature (), the present model predicts the existence of chemotactic wavelets that arise and propagate within the plateau region of the PDGF concentration profile (). Thus, in the present model, significant chemotaxis occurs throughout the wound (), as cells that enter the simulated wound at later times are able to create and follow the mesoscopic waves to explore more of the space (). These results show that even modest dose-response robustness of the circuit is sufficient to affect chemotaxis in tissues, where chemoattractant gradients are dynamic. Indeed, it is the dynamic nature of the gradients that makes this possible.

Simulation of wound invasion with chemotactic switching based on DAG/PKC polarization.

(A) Depiction of individual cells in the hybrid simulation, in which the concentration of PDGF is modeled as a continuum. The cells are initially seeded in the adjacent dermis; position along the vertical dimension is for visualization only. As time elapses, the cells invade and populate the clot region through a combination of directed migration and proliferation (see also ). (B) Depth of penetration into the wound at t = 10 days for random migration only or with chemotaxis allowed (mean ± s.d., n = 10 simulations each). (C) Spatial profiles of cell densities (randomly or chemotactically migrating) and dimensionless PDGF concentration. (D) Zoomed-in view of the cell density and PDGF concentration profiles for 1 mm length of the clot region, showing chemotactic wavelets.

Discussion

Novel ways that lipid signaling might enhance chemotactic gradient sensing, corroborated by modeling

While the bulk of the signal transduction literature focuses on proteins and their regulation by phosphorylation, it is easy to forget the significance of lipid second messengers. In the context of spatial sensing and directed cell migration, phosphoinositides and other plasma membrane lipids are uniquely positioned, owing to their peripheral location, slow diffusion, relative abundance, and diversity of lipid-protein interactions [14,34]. For more than two decades, research on lipid-mediated mechanisms of chemotactic gradient sensing has emphasized signaling through class I phosphoinositide 3-kinases (PI3Ks), which phosphorylate PIP2 to produce phosphatidylinositol (3,4,5)-trisphosphate [35,36]. That lipid and its breakdown product, phosphatidylinositol (3,4)-bisphosphate, interact with a host of signaling proteins that locally enhance F-actin polymerization by affecting activation of the Arp2/3 complex [37]. However, in most cell types tested, PI3K signaling is not absolutely required for chemotaxis, and, in fibroblasts and macrophages, even the Arp2/3 complex is dispensable for chemotaxis [10,11]. For chemotaxis of fibroblasts and other mesenchymal cells, PLC/PKC signaling is required, and so increased focus on this pathway and on DAG as a lipid second messenger is warranted.

The present models consider the dynamics of lipid signaling beyond DAG. Plasma-membrane DAG is phosphorylated by DAG kinases to produce PA, a lipid second messenger in its own right [17,18]. With PFL 1, we considered the ability of PA to enhance the observed affinity of PLCγ for PIP2 [23], presumably by increasing the lifetime of enzyme binding to membranes. We note that PLCβ1 activity, which is fostered by G protein-coupled receptor signaling, is also promoted by PA [38]. In the model, this putative feedback locally enhances PLC recruitment by activated receptors and the hydrolysis of PIP2. This effect of PA at the cell’s leading edge offsets the action of DAG kinases in metabolizing DAG. As a single feedback, PFL 1 was not sufficient to polarize signaling in the model, requiring also the recruitment and activation of PKC by DAG. The necessary action of PKC could involve either or both of two mechanisms. The more potent of the two is the neutralizing phosphorylation of MARCKS by PKC, which directly synergizes with the enhancement of PLC recruitment by enhancing the supply of PIP2 substrate. The less direct route is PFL 2, by which PKC enhances the activity of PLD; PLD hydrolyzes phosphatidylcholine to produce PA and thus presents a parallel pathway for DAG generation. Highlighting the weaker influence of PFL 2 on polarization in the model, it was necessary to add positive cooperativity with respect to active PKC and to fine-tune the values of the PFL 2-associated parameters. As one might expect, including all three feedback mechanisms in the model, with the provisions for PFL 2, yielded the most sensitive and robust gradient sensing circuit.

Robustness and reversibility of polarization

The key challenges met to varying degrees by our models is to achieve polarization in shallow gradients and for an appreciable range of chemoattractant concentration (dose-response robustness). These two aspects go hand in hand for our models: greater sensitivity and degree of polarization were generally accompanied by greater dose-response robustness. These positive aspects could be enhanced, through changes to the parameter values or addition of a feedback, but with loss of reversibility; the system is polarized and tends to remain locked in place even when the gradient is reversed. This observation is consistent with analyses of other gradient sensing models that show a tradeoff between amplification and self-locking behavior [39-42]. Self-locking is generally considered undesirable, because it precludes a prompt response to dynamic changes in the chemoattractant gradient. To adjust the directionality of movement, a cell with its chemotactic signaling polarity locked would have to execute a turn. For a gradual response such as wound invasion, such a system might be adequate or even advantageous. Obviously, with a 1D wound geometry, persistent chemotaxis would yield the maximum rate of invasion. But even with instantaneous switching between polarization states, our hybrid invasion model showed efficient invasion, with chemotactic waves arising throughout the volume of the wound occupied by cells.

As suggested previously, a gradient sensing circuit capable of polarizing in response to shallow gradients can also polarize in response to an asymmetric cell geometry [16]. This is not a generality, however. In the context of this model, it is particular to the proposed regulation of MARCKS, where the local ratio of plasma membrane area to cytosolic volume affects the dynamics. This follows basic ideas discussed in detail elsewhere [43-47]. If sensitivity to morphology were significant, it could act as an additional feedback mechanism. On the other hand, from the perspective of sensing external gradients, any such intrinsic cue might be considered a detriment.

Testable hypotheses guided by the models

Analysis of the present models offers a guide for experiments designed to test certain predictions, with polarization of DAG in chemotaxing cells as the essential readout. At a phenomenological level, a basic prediction of the models is that there is a minimum threshold of receptor activation, determined by the chemoattractant concentration, below which the pathway cannot polarize. Considering the heterogeneity in a cell population, for a suitable gradient one should expect to find subpopulations that do or do not show DAG polarization. For the shallow gradients typically achieved in chemotaxis chambers, an upper bound on receptor activation is also predicted. Whether or not either the minimum threshold or upper bound can be exceeded would depend on a cell’s expression level of the cognate receptor.

Regarding molecular mechanisms, the most basic prediction is the requirement for PKC activity, which may be readily tested using isoform-selective inhibitors. The regulation of MARCKS requires PKC activity, whereas modulation of PLD requires interaction with membrane-associated PKC but not PKC kinase activity [48,49]. From there, perturbations affecting MARCKS and PLD may be pursued. The line of experimentation outlined above does not address PFL1, which does not rely on PKC. The prediction most directly related to PFL 1 concerns the role of DAG kinases, which consume DAG but produce PA in the process. It is predicted that even a partial inhibition of DAG kinase activity (achieved either pharmacologically [50,51] or by depleting DAG kinase isoforms), while increasing DAG on a whole-cell level, would abrogate polarization for all doses of chemoattractant. Conversely, enhancement of DAG kinase activity, which might be achieved by overexpressing one or more of the isoforms, is predicted to enhance the dose-response robustness of polarization by increasing the upper bound of receptor activation. If true, more cells would be able to polarize at higher concentrations of chemoattractant.

Finally, the effect of cell geometry could be tested by co-expression of or labeling with a fluorescent, cytosolic volume marker. Imaged by epifluorescence, the marker would quantify the local height of the cytoplasm; the lower this quantity in a lamellipod, the higher the membrane area/volume, and the more likely the DAG is to be enriched there if the predicted geometry effect is significant.

Limitations of the models

The purpose of these models is to provisionally test the plausibility of the proposed feedback mechanisms. Modeling of stochastic effects and three-dimensional cell geometries with moving boundaries would make the simulations more realistic, with far greater computational expense. The topology of the signaling network is also by no means complete. As already noted, certain details of the mechanisms considered are incompletely understood, and other mechanisms not considered here might also contribute. Regarding the latter, one might consider other possible roles of the lipids involved. Whereas we focused on the putative influence of PA on PLC, PA has also been reported to enhance phosphatidylinositol 4-phosphate 5- kinase activity [52], which could boost the rate of PIP2 generation. In principle, we would expect this effect to amplify PFL 1. PIP2 also functions in other capacities that could impact PLC/PKC signaling and chemotaxis. For one, PIP2 is the preferred substrate of class I PI3Ks. Although PI3K signaling is not absolutely required for mesenchymal chemotaxisFormatting… please wait, PI3K competes with PLC for the common substrate and could be influenced by MARCKS [53]. PIP2 is also bound by the actin-modifying proteins cofilin and profilin at the plasma membrane [54]. Local modulation of free PIP2 and its hydrolysis by active PLC is expected to impact cofilin and profilin functions in chemotaxis [55,56]. Finally, we note that PIP2 is a cofactor for PLD, influencing the enzyme’s membrane localization and activity [57].

One of the key experimental observations that we have yet to adequately explain is the ability of phorbol ester (a DAG mimic) to elicit fibroblast chemotaxis when presented as a shallow gradient, even in PLCγ1-null cells [11]. This scenario is akin to the model with PFL 2 only and no MARCKS, for which we could not find any set of parameter values that yielded polarization of PKC activity in response to a shallow (≤ 10%) gradient. Whether or not a phorbol ester gradient polarizes DAG/PKC is presently unknown; it is possible that chemotaxis occurs regardless. If polarization does occur in PLC-null cells executing phorbol ester chemotaxis, the nature of the underlying mechanisms will need to be clarified.

Materials and methods

Models of PLC/PKC signaling

The reaction-diffusion models are composed of partial differential equations and associated boundary and initial conditions. The species in the models and their interactions are described here; mathematical details and justifications for chosen parameter values are given in . Within that document, model species, their diffusivity values, and initial conditions are listed in S1 Table, and rate equations and base-case values of rate constants are listed in S2 Table. Cytosolic species, signified by capital letters, have local concentrations in μM, whereas membrane species, signified by lowercase letters, have local densities in #/μm2. Rather than explicitly model the ligand-receptor dynamics, in this work we assume a linear profile of active receptors (r) as the model input, calculated from .

As explained under Results, the prefactor determines the midpoint value of r; when the dimensionless rfrac = 1, the midpoint value of r = 130/μm2, corresponding to approximately 1 x 105/cell. The distance from the midpoint, x, is in μm, and it is scaled by the cell length of 40 μm. Therefore, rsteep is the fractional gradient steepness; a value of rsteep = 0.1, for example, corresponds to a 10% difference between the values of r at the extreme front and back of the cell.

Referring to the diagram shown in , activated receptors recruit inactive PLC enzyme (E) from the cytosol to form active PLC at the membrane (e). Active PLC hydrolyzes PIP2 (free-density p). In the absence of stimulation, the density of PIP2 is maintained through basal synthesis and turnover. PIP2 hydrolysis by PLC generates the lipid second messenger, DAG (d). DAG is phosphorylated by DAG kinases to produce PA (d), which can be dephosphorylated by phosphatidic acid phosphatases to recover DAG [22,58]; these are modeled as pseudo-first-order reactions. Like PIP2, PA is also subject to basal generation and consumption; however, basal generation of PA is intentionally very low. DAG recruits catalytically competent but inactive PKC (C) from the cytosol by reversible binding of the tandem C1 domain of PKCs, forming active PKC (c*) at the membrane [59]. Membrane-bound PKC is sensitive to dephosphorylation, and so the active PKC is converted to an inactive, membrane-bound form (c); this species is either autophosphorylated to regenerate active PKC or it dissociates to join the cytosolic pool. MARCKS is an abundant substrate of PKC that is present in both cytosolic and membrane-bound forms [14]. The unphosphorylated, cytosolic form (M) inserts into the plasma membrane via its myristoyl group and interacts with PIP2 via its effector domain; the latter is a high-avidity, electrostatic interaction, and thus, MARCKS sequesters a substantial fraction of the intracellular PIP2 [13,60]. The membrane-bound forms of MARCKS (m) are phosphorylated by active PKC, causing loss of affinity for the plasma membrane and liberation of PIP2; the phosphorylated form of MARCKS (M) is cytosolic [12]. MARCKS is dephosphorylated in the cytosol by a pseudo-first-order reaction to complete the cycle.

Imposed on this basic network structure are two additional PFLs. As described under Results, PFL 1 considers that the PA engages PLC in complex with activated receptors and thus extends the enzyme’s mean lifetime (or, reduces its effective off-rate) at the plasma membrane. The net rate of PLC recruitment, V, is given by .

In this equation, the parameters characterizing PFL 1 are K, the equilibrium constant of PA-PLC interaction, and ε, a dimensionless escape probability (see ). With either K set to zero or ε set to 1, PFL 1 is turned off. PFL 2 considers that active PKC engages and enhances the activity of PLD. This influence on the rate of PA synthesis, V, is modeled as a Hill function in .

In this equation, the parameters characterizing PFL 2 are K, a saturation constant, γ, a dimensionless gain parameter, and n, the Hill coefficient. With either K set to zero or γ set to 1, PFL 2 is turned off.

Implementation of PLC/PKC models

The partial differential equation models are implemented in Virtual Cell (http://www.vcell.org), a computational environment for modeling and simulation in cell biology [61]. The Biomodel and primary simulations are publicly available in Virtual Cell under user name jnosbis, Biomodel name ‘Nosbisch chemotaxis 2020’. As in Mohan et al. [16], the base two-dimensional geometry is an ellipse with a major axis length of 40 μm and minor axis length of 1.8 μm (); these dimensions were selected according to the important dimensions of a migrating fibroblast. For the asymmetric geometry considered in , the left end of the ellipse was blunted by changing the y-component of the ellipse equation to have an exponent of 6, rather than 2, for negative values of x.

Virtual Cell uses a finite-volume method to numerically solve the reaction-diffusion equations, and the Fully-Implicit Finite Volume (variable time step) solver was used, with voxel dimensions of Δx = Δy = 0.1 μm and a maximum time step of 0.1 s. We confirmed that the computational results were not significantly affected by modest changes to those values. The initial conditions are such that the system is stationary in the absence of stimulation. Simulations were run for 20,000 s to ensure either a steady state or sustained oscillations. For all simulations, the reversibility of the spatial pattern was assessed by reversing the gradient of PDGF receptor activation across the cell after a steady state was reached and extending the simulation for another 20,000 s. If the spatial pattern achieved a steady state but failed to reverse, then we considered that simulation to display locking behavior, and the simulation is marked with an asterisk.

Wound healing model

The model of collective cell migration is adapted from the hybrid simulation strategy of Monine and Haugh, described in detail previously [28]. Briefly, the dimensionless PDGF concentration, u, is treated as a continuum, subject to synthesis (only in the clot), interstitial diffusion, intrinsic degradation, and degradation by cells (receptor-mediated endocytosis). Local receptor activation, r, is given by the algebraic function, r = u2/(1 + u + u2), based on a quasi-steady-state approximation [29]; therefore, u = 1 corresponds to r = 1/3. The local value of r, along with the calculated cell density, determines the cell-mediated degradation of PDGF and the net rate of cell proliferation or death. The original code was adapted from C to MATLAB; the version that was implemented to generate the results shown in is provided in the Supporting Information (. The difference in the new model is that the migration behavior, rather than depending on intracellular signaling in an analog fashion, switches between random and chemotactically biased states. In the random state, the cells move with a diffusivity D = 3x10-4 mm2/h. Considering a persistence time of ~ 0.5–1 h for fibroblasts, this corresponds to a cell speed of ~ 0.03 mm/h (~ 0.5 μm/min) [62]. In the chemotactic state, the cell movement is governed by diffusion and convection. For the base case, the chemotactic (convective) velocity in the direction of the PDGF gradient, S, was chosen as 0.03 mm/h to match the cell speed estimate given above. The random component of the chemotactic cells’ movements was set to a lower value of D = 1x10-4 mm2/h to represent noise in the cell’s movement up-gradient. It is envisioned that the cells’ movements in the other two orthogonal dimensions of the tissue would be comparable regardless of migration state, and therefore a rough estimate of the associated chemotactic (forward migration) index is ~ 1/3, a modest value considering that measured population estimates of this forward migration index are ≈ 0.2 and include tracks of chemotaxing and non-chemotaxing cells [11].

The switching between random and chemotactic migration were taken from the approximate cut-offs between non-polarizing and polarizing conditions indicated in , supplemented with additional data for rsteep = 0.2, 0.4, and 0.5. For each value of rsteep, the lower and upper bounds of rfrac were approximated (at rsteep = 0.67, it was found that the upper bound, if it exists, is >> 1), with oscillations considered non-polarizing. These data were used to construct a relationship between the value of rfrac (the variable r calculated in the wound healing model) and the corresponding minimum value of rsteep required for polarization; in the wound healing model, rsteep was equated to Δr/r, where Δr is the difference in r values evaluated +/- 20 μm from the midpoint (matching the cell length in the PLC/PKC models). The relationship was fit to , reminiscent of an inter-molecular potential.

The first term, 0.01, is an homage to the notion that the limit of spatial gradient sensing is a ~ 1% difference across a cell’s length [63]. The other constant parameter values are a = 1.7, r = 0.12, n = 2, and m = 1.5. A plot of Eq 4 with these parameter values is shown in

Modeling details.

(PDF)

Click here for additional data file.

MATLAB code used to simulate wound invasion (Nosbisch_et_al_wound_sim.m).

(M)

Click here for additional data file.

Time courses of all model species for Fig 2A, 10% steepness and rfrac = 0.1.

(PDF)

Click here for additional data file.

Transition from oscillations to stable pattern as rfrac is increased (10% steepness).

(PDF)

Click here for additional data file.

Phase plots for different rfrac values.

These correspond to the simulations analyzed in with 10% gradient steepness. The plot with rfrac = 0.1 is the same as in . When rfrac = 0.0278, DAG is low at both the front and back of the cell; with rfrac = 0.278, DAG is high at both the front and back of the cell.

(PDF)

Click here for additional data file.

Systematic variation of kinetic parameters (no PFL 2).

(PDF)

Click here for additional data file.

Variation of key parameters with PFL 2 included in the system.

The parameters targeted were k, which affects the conversion of DAG to PA; K, which affects the PA-PLC affinity in PFL 1; and K, which affects the sensitivity of PFL 2 to active PKC. Each parameter was decreased to 0.3X and increased to 3X. Each plot shows active PKC density vs. rfrac for 10% gradient steepness, as in .

(PDF)

Click here for additional data file.

Plot of Eq 4.

For each cell i in the simulation, its mean receptor activation r and the difference in receptor activation across its length Δr were calculated. Cells with receptor activation states above the curve engaged in chemotaxis, whereas those with states below the curve engaged in random migration only.

(PDF)

Click here for additional data file.

Animation of the simulation depicted in Fig 7A.

(AVI)

Click here for additional data file.

13 Sep 2019

Dear Dr Haugh,

Thank you very much for submitting your manuscript 'Mechanistic models of PLC/PKC signaling implicate phosphatidic acid as a key amplifier of chemotactic gradient sensing' for review by PLOS Computational Biology. Your manuscript has been fully evaluated by the PLOS Computational Biology editorial team and in this case also by independent peer reviewers.

The reviewers appreciated the attention to an important problem and appreciate the soundness of the model and the insights gained from it. However, they raised substantial concerns about the manuscript as it currently stands. While your manuscript cannot be accepted in its present form, we are willing to consider a revised version in which the issues raised by the reviewers have been adequately addressed.  We cannot, of course, promise publication at that time.

In addition to addressing the reviewers' points aimed at clarifying the results, you can more prominently convince readers of the originality of the current model in the context of your previous work.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

Your revisions should address the specific points made by each reviewer. Please return the revised version within the next 60 days. If you anticipate any delay in its return, we ask that you let us know the expected resubmission date by email at ploscompbiol@plos.org. Revised manuscripts received beyond 60 days may require evaluation and peer review similar to that applied to newly submitted manuscripts.

In addition, when you are ready to resubmit, please be prepared to provide the following:

(1) A detailed list of your responses to the review comments and the changes you have made in the manuscript. We require a file of this nature before your manuscript is passed back to the editors.

(2) A copy of your manuscript with the changes highlighted (encouraged). We encourage authors, if possible to show clearly where changes have been made to their manuscript e.g. by highlighting text.

(3) A striking still image to accompany your article (optional). If the image is judged to be suitable by the editors, it may be featured on our website and might be chosen as the issue image for that month. These square, high-quality images should be accompanied by a short caption. Please note as well that there should be no copyright restrictions on the use of the image, so that it can be published under the Open-Access license and be subject only to appropriate attribution.

Before you resubmit your manuscript, please consult our Submission Checklist to ensure your manuscript is formatted correctly for PLOS Computational Biology: http://www.ploscompbiol.org/static/checklist.action. Some key points to remember are:

- Figures uploaded separately as TIFF or EPS files (if you wish, your figures may remain in your main manuscript file in addition).

- Supporting Information uploaded as separate files, titled Dataset, Figure, Table, Text, Protocol, Audio, or Video.

- Funding information in the 'Financial Disclosure' box in the online system.

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see here.

We are sorry that we cannot be more positive about your manuscript at this stage, but if you have any concerns or questions, please do not hesitate to contact us.

Sincerely,

Stacey Finley, Ph.D.

Associate Editor

PLOS Computational Biology

Daniel Beard

Deputy Editor

PLOS Computational Biology

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors present an expanded version of a model published previously. This version explores the role of two additional positive feedback loops not considered in the old model.

The model analysis and simulations generated several interesting predictions that can be experimentally validated.

Major comments:

1- In the present form, the model is not validated, as the authors did not perform any comparison of model simulations with experimental data. I understand that ideal experimental data to validate the model may not be available, but the authors should, at least, compare PIP2, DAG and PA steady state levels in resting conditions (without chemotactic gradient) with available estimates (even if in different cells) to check that the model is working around a feasible range of concentrations.

2- The authors should also add a sensitivity analysis, so that we can have an idea of the influence of parameter values in model behaviour. The authors have studied the impact of some parameter changes that are more relevant for the studied PFLs, but a overall sensitivity analysis would allow us to identify parameters that should be estimated with higher precision, or detect problems of identifiability. I suggest the authors should compute the relative sensitivity of "Maximum front/back PKC activity ratio" and "dose-response robustness" to small relative changes in each parameter value.

3- The authors have studied a) PFL1 parameter effects and b) geometric effects in the absence of PFL2. When they study PFL2, they identify some parameter conditions where PFL2 does (or does not) contribute to PKC activity polarization. It would be nice to check (even if in supplementary materials only) if the observations made in the absence of PFL2 (both for a) and b) effects) are also seen in both PFL2 scenarios (PFL2 parameters that influence or not PKC activity polarization). If the observations are different, the differences can be a hint to identify the presence and type of PFL2 in a particular cellular system.

4- PA is also known to activate PI4P5K (https://onlinelibrary.wiley.com/doi/full/10.1002/jcb.21027), an enzyme that produces PIP2. So, PA may activate PIP2 synthesis rate. How would this regulatory interaction interfere with the remaining feedback loops in your system?

Minor comments:

1- The meaning of the asterisk (identifying irreversibility of the PKC polarization) should me mentioned in all the figure legends where it appears, and not only in the main text

2- Is it reasonable to have the same diffusivity for DAG and DAG-PKC?

3- In figure 3 C, in the second column, some lines have two numbers (1/16.7). Explain the meaning of the two numbers in the figure legend.

Reviewer #2: In this manuscript, the authors build on two previous model (Refs 17 and 29) to investigate the role of two new feedback loops associated with phosphatidic acid. They show that components of this feedback loops are important regulators of the model and can enhance the robustness of DAG/active PKC polarization. Using this network in a model of wound-healing, they show that efficient collective migration can be achieved with these polarization thresholds. This is an important area of research but the execution of the modeling was lacking.

Major comments (these comments are in the order of my reading the text):

1) Originality: It is very hard to distinguish the substantial new novel contributions from this model as compared to Mohan et al. 2017 [ref 17]. In fact, Figure 1A in very similar in both papers except for the addition of the red reactions in the current manuscript. It is not immediately obvious that the addition of the two new feedback loops as modeled by the reactions in red is a significant advance for publication in Plos Comp Biol. For a reader to understand this manuscript requires a thorough reading of ref 17 making it hard to be a standalone publication.

2) There are many abbreviations and assumptions throughout the text limiting accessibility. For example, in Results paragraph 1, the reasoning underlying how PA affects PLC activity were modeled are really assumptions of the model that must be carefully stated in a separate section. Additionally, these assumptions must be revisited in the limitations section of the discussion. The large number of abbreviations and notational issues also make the manuscript hard to follow. For example, using 'e' for active PLC enzyme is a weird notation. For the longest time, I interpreted to be Euler's number and struggled to follow the terms in the supplementary table. A list of abbreviations would also be helpful, although I'm sure the two feedback loops can be named better than PFL1 and PFL2.

3) Figure 1B: Where is the gradient of the receptors in the spatial geometry? If it is a 2D spatial model (as shown in Figure 3A), then is the receptor gradient simply the difference of receptor at x = 20 and x =-20 divided by 40? Or is it along the arc length? In which case, x doesn't span -20 to +20 but whatever the arc length corresponding to the locations is. Given that it is a 2D spatial model, and the membrane is only the arc length of the curve, the physicality of the imposed receptor gradient doesn't make sense to me.

4) Page 7: This conclusion held true across across simulations with extensive sweeps of the PFL parameters. Where are these shown?

5) Figure 2: The authors should consider using different symbol shapes and not just colors for front and back (this is a minor point). In panel 2A: do the front and back curves overlap throughout? What is most interesting is that the region close to rfrac of 0.1 is where many oscillatory solutions seem to lie. However, these oscillations are not fully explored. Does this oscillation translate into a spatial pattern? These regimes should be carefully analyzed for stability features since these oscillations are observed in simulations and the underlying cause identified. The issue of oscillations was completely glossed over in the manuscript. Could perturbations in rfrac result in loss of steady state or the maximum front/back ratio? Could stochasticity be important (likely yes, but how?)? Addressing these questions would complete the study and significantly elevate the manuscript.

6) Figure 2: Another issue is that membrane density of PKC front to back is compared but what does this mean in a spatial model (see comment regarding Figure 1B, the same issue lies here and throughout the manuscript.)

7) Figure 3 and associated results: This the first place where the geometry is shown and one asymmetry is introduced. The logic for introducing this asymmetry is not clear. But what becomes clear immediately is that shape effects are important. However, this effect is also not fully explored. Given the 2D model (1D membrane) and the role of shape in signaling (see the many references in ref 46, the authors just gloss over this important issue and move on to other things. This significantly diminishes the potential impact of the work. Also, the difference between the two columns in figure 3C is not clear.

8) Figure 4 and 5 and associated results: The authors identify a series of kinetic parameters that are important in determining the ratio of PKC front/back activity. Then they vary these parameters and show that they are important. This is a circular argument. These parameters are important because of the way the model was constructed. Additionally, it is not clear under what physiological conditions these parameters can be measured or manipulated in experiments. Varying copy number of enzymes is a much more meaningful experimental correlate. What are the *s in Figure 4? They are not defined in the caption of this or previous figures.

9) sensitivity analysis: sensitivity analysis of the kinetic parameters and initial conditions should be included, particularly in light of the oscillatory regimes. What if the diffusion constants were altered? Then the results across the entire manuscript would change I would think? It would be important to report to what extent the model is robust to these changes. Plus the diffusion constants must be reported from literature where available.

10) if receptors were allowed to diffuse freely rather than maintain an active gradient across the membrane (after correcting for the geometry issues raised above) would a steady state receptor gradient exist? That seems unlikely given that studies have shown that receptors diffuse in the plane of the membrane and diffusion serves to dissipate the gradient. This should be explored a carefully because this gradient is the input to the model.

11) Testable hypotheses: Regarding the experiments proposed by the authors to test molecular mechanisms, it is not as trivial as they propose to use isoform selective inhibitors. First, the model predicts that kinetic parameters are important, which are near impossible to alter in situ. Interestingly, this section has the least references to plausible experimental modalities, so it is not clear if these methods are even possible.

11) References: the fraction of self citation is rather high. I appreciate that the authors have done a substantial amount of work in this field but they have also missed many fundamental citations in cell shape regulation of signaling (see comment about the references in 46) and receptor diffusion.

Reviewer #3: Please see the attached word document.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: None

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

Submitted filename: review.docx

Click here for additional data file.

29 Nov 2019

Submitted filename: Comments to the Authors_PLOSCB_JN_2019_JH.docx

Click here for additional data file.

6 Jan 2020

Dear Dr Haugh,

Thank you very much for submitting your manuscript, 'Mechanistic models of PLC/PKC signaling implicate phosphatidic acid as a key amplifier of chemotactic gradient sensing', to PLOS Computational Biology. As with all papers submitted to the journal, yours was fully evaluated by the PLOS Computational Biology editorial team, and in this case, by independent peer reviewers. The reviewers appreciated the attention to an important topic, and the edits made in the first round of revision. There remains one issue to be addressed - further clarifying what is meant by the use of a linear gradient for the cell surface receptors, particularly in the context of your spatial model.

We would therefore like to ask you to modify the manuscript according to the review recommendations before we can consider your manuscript for acceptance. Your revisions should address the specific points made by each reviewer and we encourage you to respond to particular issues Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

In addition, when you are ready to resubmit, please be prepared to provide the following:

(1) A detailed list of your responses to the review comments and the changes you have made in the manuscript. We require a file of this nature before your manuscript is passed back to the editors.

(2) A copy of your manuscript with the changes highlighted (encouraged). We encourage authors, if possible to show clearly where changes have been made to their manuscript e.g. by highlighting text.

(3) A striking still image to accompany your article (optional). If the image is judged to be suitable by the editors, it may be featured on our website and might be chosen as the issue image for that month. These square, high-quality images should be accompanied by a short caption. Please note as well that there should be no copyright restrictions on the use of the image, so that it can be published under the Open-Access license and be subject only to appropriate attribution.

Before you resubmit your manuscript, please consult our Submission Checklist to ensure your manuscript is formatted correctly for PLOS Computational Biology: http://www.ploscompbiol.org/static/checklist.action. Some key points to remember are:

- Figures uploaded separately as TIFF or EPS files (if you wish, your figures may remain in your main manuscript file in addition).

- Supporting Information uploaded as separate files, titled 'Dataset', 'Figure', 'Table', 'Text', 'Protocol', 'Audio', or 'Video'.

- Funding information in the 'Financial Disclosure' box in the online system.

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com  PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

We hope to receive your revised manuscript within the next 30 days. If you anticipate any delay in its return, we ask that you let us know the expected resubmission date by email at ploscompbiol@plos.org.

If you have any questions or concerns while you make these revisions, please let us know.

Sincerely,

Stacey Finley, Ph.D.

Associate Editor

PLOS Computational Biology

Daniel Beard

Deputy Editor

PLOS Computational Biology

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: I am satisfied with the author responses

Reviewer #2: The revised version of this manuscript is significantly improved and I commend the authors for their effort. I only have the following point for clarification.

Figure 1B and Figure 2: I still don't understand what it means to describe a linear gradient of receptors when they are located on the arc length of the membrane. I appreciate the authors clarification from the figure but since this is one of the important features in their model, it should be stated very clearly what this means. I am very confused by this point and I think clarity on this will not only help eliminate confusion but increase accessibility of spatial models.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

3 Feb 2020

Dear Prof. Haugh,

We are pleased to inform you that your manuscript 'Mechanistic models of PLC/PKC signaling implicate phosphatidic acid as a key amplifier of chemotactic gradient sensing' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch within two working days with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Stacey Finley, Ph.D.

Associate Editor

PLOS Computational Biology

Daniel Beard

Deputy Editor

PLOS Computational Biology

***********************************************************

31 Mar 2020

PCOMPBIOL-D-19-01266R2

Mechanistic models of PLC/PKC signaling implicate phosphatidic acid as a key amplifier of chemotactic gradient sensing

Dear Dr Haugh,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Laura Mallard

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol