https://orcid.org/0000-0003-4141-4034
Figures
Abstract
During the last decade, a consensus has emerged that the stochastic triggering of an excitable system drives pseudopod formation and subsequent migration of amoeboid cells. The presence of chemoattractant stimuli alters the threshold for triggering this activity and can bias the direction of migration. Though noise plays an important role in these behaviors, mathematical models have typically ignored its origin and merely introduced it as an external signal into a series of reaction-diffusion equations. Here we consider a more realistic description based on a reaction-diffusion master equation formalism to implement these networks. In this scheme, noise arises naturally from a stochastic description of the various reaction and diffusion terms. Working on a three-dimensional geometry in which separate compartments are divided into a tetrahedral mesh, we implement a modular description of the system, consisting of G-protein coupled receptor signaling (GPCR), a local excitation-global inhibition mechanism (LEGI), and signal transduction excitable network (STEN). Our models implement detailed biochemical descriptions whenever this information is available, such as in the GPCR and G-protein interactions. In contrast, where the biochemical entities are less certain, such as the LEGI mechanism, we consider various possible schemes and highlight the differences between them. Our simulations show that even when the LEGI mechanism displays perfect adaptation in terms of the mean level of proteins, the variance shows a dose-dependence. This differs between the various models considered, suggesting a possible means for determining experimentally among the various potential networks. Overall, our simulations recreate temporal and spatial patterns observed experimentally in both wild-type and perturbed cells, providing further evidence for the excitable system paradigm. Moreover, because of the overall importance and ubiquity of the modules we consider, including GPCR signaling and adaptation, our results will be of interest beyond the field of directed migration.
Author summary
Though the term noise usually carries negative connotations, it can also contribute positively to the characteristic dynamics of a system. In biological systems, where noise arises from the stochastic interactions between molecules, its study is usually confined to genetic regulatory systems in which copy numbers are small and fluctuations large. However, noise can have important roles when the number of signaling molecules is large. The extension of pseudopods and the subsequent motion of amoeboid cells arises from the noise-induced trigger of an excitable system. Chemoattractant signals bias this triggering thereby directing cell motion. To date, this paradigm has not been tested by mathematical models that account accurately for the noise that arises in the corresponding reactions. In this study, we employ a reaction-diffusion master equation approach to investigate the effects of noise. Using a modular approach and a three-dimensional cell model with specific subdomains attributed to the cell membrane and cortex, we explore the spatiotemporal dynamics of the system. Our simulations recreate many experimentally-observed cell behaviors thereby supporting the biased-excitable network hypothesis.
Citation: Biswas D, Devreotes PN, Iglesias PA (2021) Three-dimensional stochastic simulation of chemoattractant-mediated excitability in cells. PLoS Comput Biol 17(7): e1008803. https://doi.org/10.1371/journal.pcbi.1008803
Editor: Martin Meier-Schellersheim, National Institutes of Health, UNITED STATES
Received: January 29, 2021; Accepted: June 8, 2021; Published: July 14, 2021
Copyright: © 2021 Biswas et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the manuscript and its Supporting information files.
Funding: This work was funded in part by DARPA, under contract number HR0011-16-C-0139 (P.A.I.) and the National Institutes of Health grant NIH/NIGMS R35 GM118177 (P.N.D.). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
How cells sense and interpret external chemoattractant cues and use this information to direct cell movement is one of the most fundamental processes of biology. Unicellular organisms rely on this mechanism to seek nutrients and survive. In multicellular organisms, it is a fundamental process during embryonic development [1, 2] as well as responsible for the proper operation of the mammalian immune system [3, 4]. Perversely, it is through the development of directed cell migration that cancer cells become metastatic [5–7].
Mathematical models have been fundamental in elucidating the mechanisms that cells use to direct the cell migration [8–10]. There is a broad consensus that cells such as the social amoeba Dictyostelium discoideum and mammalian neutrophils sense the chemoattractant gradient through a local excitation, global inhibition (LEGI) mechanism based on an incoherent feedforward loop motif that was originally proposed to explain perfect adaptation [11–13]. By incorporating different diffusion properties on the signal components, the mechanism senses static spatial gradients without movement [12, 13].
In response to a spatially uniform stimulus, cells display an initial transient response in which Ras and downstream PI(3,4,5)P3 and F-actin activities increase and decrease several times, before eventually returning close to the pre-stimulus basal states, resulting in “near-perfect” adaptation. Signaling motifs responsible for perfect adaptation fall into one of two broad classes: a negative feedback (NFB) loop with a buffering node, or an incoherent feedforward (IFF) loop with a proportioner node [14]. In Dictyostelium, experimental evidence favors the presence of IFF-based adaptation [12, 15–17]. The LEGI mechanism, a form of IFF, assumes the existence of a fast local excitor and a slow globally diffusive inhibitor. The productions of the excitor and the inhibitor are independently driven by the receptor occupancy [13]. A local response regulator, which is activated by the excitor and inhibited by the inhibitor, drives downstream signals.
LEGI mechanisms explain gradient sensing but do not account for several aspects of the chemotactic cells, including the ability to move in the absence of external cues. Excitable systems recreate many of the observed properties of randomly migrating cells [18–23] including the stereotypic nature of pseudopods during migration, as well as the spatial pattern of activities exhibited by both signaling and cytoskeletal elements in cells [24, 25]. Recently, LEGI mechanisms have been coupled to an excitable system to explain how a cell’s ability to migrate randomly can be steered in the direction of the external gradient [17, 22]. When combined with a memory-like ability to polarize the chemotactic machinery, these models account for nearly all the observed behavior of chemotaxing cells [22, 26]. Most mathematical models of chemoattractant signaling have adopted a variant of the FitzHugh-Nagumo (FHN) model of neuronal excitability [27, 28]. These models, however, suffer from limitations because of the phenomenological aspect of the model. For example, using reaction terms in polynomial form as in the FHN model does not permit direct biological interpretability. Moreover, system states, which would typically represent concentrations, are not constrained to be nonnegative [19–22, 29].
Though the LEGI mechanism successfully explains the temporal and spatial responses to chemoattractant in Dictyostelium [15–17, 30], little attention has been paid to the effect of noise on the LEGI mechanism. Moreover, in the excitable paradigm, the ability to generate random protrusions, as seen in unstimulated migrating cells, relies on stochastic fluctuations triggering the excitable system. A proper account of these fluctuations is vital, however, it is their relative size that determines whether the external signal directs movement properly. In practice, noise arises as an intrinsic feature of the stochastic nature of the biochemical reactions and depends on the state of the dynamical system [31, 32]. To our knowledge, however, all existing computational models use partial differential equations and generate these fluctuations by injecting noise as an additional input into these differential equations.
To overcome the aforementioned limitations, here we present a new model of the biological signaling mechanism that regulates motility. Specifically, to account for the noise accurately we eschew the partial-differential equation approach and instead incorporate the reaction-diffusion processes into the URDME (Unstructured Reaction Diffusion Master Equation) software [33]. This methodology does not make a priori assumptions on the size of the stochastic perturbations; instead, the fluctuations are inherent in the underlying chemical master equation. Hence, the noise is controlled by number of the molecules of the reacting species. Moreover, we consider a realistic geometry consisting of a three-dimensional cell with membrane and cortex elements. Finally, we use detailed biochemical models of the receptor dynamics, LEGI and excitable modules and verify the sufficiency of these three interconnected networks in explaining the spatiotemporal dynamics of chemotactic signaling.
Methods
Stochastic simulations
The Unstructured Reaction Diffusion Master Equation (URDME) approach [33] adopted in the present study uses the Next Sub-volume Method (NSM) [34], which is a variant of the Next Reaction Method (NRM) [35]. The URDME approach is computationally efficient and easily applicable to unstructured meshes allowing it to consider complex geometries [33, 36, 37]. In the URDME framework, the time when the next reaction/diffusion event takes place is computed using Gillespie’s direct method [38]. Spatial heterogeneity is achieved by discretizing the simulation domain into tetrahedral voxels. For determining in which voxel an event occurs, NRM uses an event queue. If the event is a chemical reaction, then only the values of the respective species in that voxel are modified. Alternatively, if the event represents diffusion, then values in both the “sending” and “receiving” voxels are updated. In this formulation, the simulation time is proportional to the logarithm of the number of reactions [35]. Specific examples of how to include different types of reactions in the URDME framework are given in S1 File.
Overall system architecture
For the present study we have divided the signaling pathways that drive the observed excitable dynamics in Dictyostelium cells into three signaling subsystems (Fig 1A): 1) The G-protein coupled receptor (GPCR) subsystem is responsible for sensing chemoattractants. 2) The receptor occupancy information from GPCR is passed down to the Local Excitation and Global Inhibition (LEGI) mechanism which is responsible for adaptation in the presence of a global stimulus as well as interpreting directional cues when cells are in a chemoattractant gradient. 3) The Signal Transduction Excitable Network (STEN), accounts for the features of the excitable behavior, including all-or-nothing responses, refractory periods and wave propagation, and provides a characteristic response both in the presence or absence of chemoattractant signals. We have excluded downstream networks that directly influence the actin dynamics. The output of our system could readily be coupled to the cytoskeletal signaling network but that would likely necessitate numerous additional entities, including elements with a large number of molecules (e.g., ∼108 molecules of actin [39]) leading to a presently unmanageable computational burden for these stochastic simulations.
(A) Reaction scheme adopted in the present study involves a receptor module describing GPCR (G-Protein Coupled Receptor) dynamics, a LEGI (Local Excitation and Global Inhibition) module that provides adaptation and directional sensing, and a STEN (Signal Transduction Excitable Network) that describes the excitable behavior of the cell. (B) Isometric (left) and cross-sectional side view (right) of the hemispherical simulation domain of radius 5 μm and thickness of 200 nm. The outer surface is the membrane and the interior of the shell is the cortex.(C) Frequency distribution of all the nodal volumes (black), membrane nodes (red) and cortex nodes (gray).
Geometry
We used a stationary hemispherical structure of radius 5 μm to capture the general shape of adherent amoeboid cells devoid of any cytoskeletal components (as if treated with Latrunculin) (Fig 1B). Because most of the species known to be involved in generating the spatial patterns are either membrane-bound, or in the cytoskeletal cortex, we focus on these areas of the cell. We do not incorporate nodes for the cytosol but consider this a sink from which some molecules, such as PKBA, a component of PKB*s, can shuttle to the membrane [40, 41]. We assumed that the cortex is 200 nm thick, and discretized a shell of this thickness into nodes resulting from an unstructured tetrahedral mesh. The allowable minimum and maximum distances between the nodes were set at 100 and 300 nm, respectively. Overall, this gave rise to 5,500 nodes and 16,469 tetrahedral elements with volume distribution 26.8 ± 4.6 × 10−4 μm3 (mean±std. dev.). Because a coarser mesh affects the smoothness of diffusion, a finer mesh would be preferred, but this increases the simulation cost considerably and would become a major constraining factor as the number of biochemical elements in the model grew (S1 Fig). For this reason, we compromised by choosing a medium size mesh so that, in an element of average size, a single molecule corresponds to a concentration of approximately 60 nM. The nodal volumes obtained from the dual of the tetrahedral mesh have a leptokurtic distribution with parameters: 80.1 ± 16.3 × 10−4 μm3 (mean±std. dev.). The membrane nodes are those close to the surface (0–20 nm); the others are cortex nodes (Fig 1C). The volume distributions of cortex and membrane nodes were similar, with means of 8.4 × 10−3 μm3 and 7.7 × 10−3 μm3, respectively. The mean subtracted distributions failed to reject the null hypothesis using a Mann-Whitney-Wilcoxon test, indicating that the two size distributions are not significantly different. The distributions of membranes nodes at the basal and apical surfaces were also similar, with mean volumes of 7.3 × 10−3 μm3 and 7.9 × 10−3 μm3, respectively.
Results
GPCR signaling module
The initiation of chemotaxis in Dictyostelium involves binding of chemoattractant molecules (ligand) to surface-bound G-protein coupled receptor (GPCR) molecules. In the present study, we focused on cAR1 as the major receptor for cAMP. The heterogeneous 3’,5’-cyclic adenosine monophosphate (cAMP) binding in Dictyostelium has been modeled using three states indicating different levels of affinities [42]. For simplicity, we ignored the sequential binding dynamics and followed the classification of GPCRs on the basis of affinity only. The three states of unoccupied GPCRs have high (H), low (L) and slow (S) binding affinity with respect to extracellular cAMP [42]. The corresponding occupied receptor states are labeled H:C, L:C and S:C, respectively. Interaction of cAR1 with cAMP also results in the desensitization of receptors due to phosphorylation [43, 44]. Despite any notable loss in the binding sites, a 3–5-fold decrease in the binding affinity of the low affinity class (L) has been reported [43]. This motivated us to consider additional phosphorylated receptor states: PH, PL, and respective occupied states: PH:C, PL:C which resulted due to interaction among receptors (H,L,H:C,L:C) and inorganic phosphate, P. Our complete receptor model consists of 10 reacting species and 26 reactions in total (Fig 2A). The detailed reactions and parameter values in mesoscopic form are listed in Table 1.
(A) Detailed schematic of the different states of the G-protein coupled receptor (GPCR) and cAMP binding. Unoccupied receptors exist in high (H) and low (L) affinities, and a third slow (S) binding state. Occupied receptors are denoted H:C, L:C and S:C. Phosphorylated states are denoted by a superscript P. (B) Dose response curve. The circled numbers denote different concentration levels of cAMP corresponding to (1) low (4%), (2) mid (50%) and (3) high (100%) levels of R.O. (C) Steady-state R.O. in response to different concentrations of cAMP. (D) Distribution of nodes based on R.O. for mid (light red) and saturating (red) cAMP doses. (E,F) Temporal profile of number of total occupied receptors (H:C+L:C+S:C+PH:C+PL:C) at a single random node (E) and in the cell (F) for low (black) and saturating (red) doses of cAMP. (G) Temporal profile of total free (black) and occupied (red) receptors in a cell in response to application and removal of the high dose of cAMP. The shaded regions denote the respective standard deviations (n = 10 independent simulations in which the parameter values were varied according to the distributions from Table 1 to account the cell-to-cell variation.
Parameter values for reaction no. 1–10 are from [42]; for reaction no. 11–26, the values were estimated to match the experimental observations from [43, 44]; for reaction no. 27–28, the values were obtained by minimizing the difference between the responses of full- and reduced-order system. The diffusion constants for all the GPCRs are assumed to be 2.7 × 10−2 μm2s−1 [45]. The total number of receptor molecules followed a Gaussian distribution with mean and standard deviation of 70,000 and 5,000 molecules, respectively.
To study the system response to different levels of cAMP, we mapped the relationship between the dose of the cAMP (in terms of the total number of molecules of cAMP experienced by the cell) and receptor occupancy (R.O.) (Fig 2B). We considered three different cAMP doses corresponding to low (4%), mid (50%) and high (100%, saturating dose) levels of R.O. The respective profiles at the basal surface show little variation at the two cAMP concentration extremes, but considerably more heterogeneity (skewness = 0.04) at the mid-point (Fig 2C and 2D). In this case, the individual nodes displayed an approximately normal distribution, but some nodes had as few as 10% (n = 4) or as high as 85% R.O. (n = 2). With higher cAMP concentrations, the distribution became more skewed (skewness = −1.77), but there was a small number of nodes (n = 18) with as little as 90% R.O. (Fig 2D).
In the presence of cAMP, the free GPCR states get converted to the respective occupied states, among which the occupied low affinity receptors (L:C) form the greater population (S2(A) and S2(B) Fig). The phosphorylated states showed slower dynamics (t1/2 = 198 s) compared to the unphosphorylated states (t1/2 = 10 s). The average Fano factor was 0.95 (S2(C) Fig). Whereas the absolute value of the temporal fluctuations at individual nodes increased with higher cAMP concentrations (Fig 2E), the relative values decreased when expressed in terms of coefficient of variation (COV = 0.69 at 4% R.O., 0.19 at 100% R.O.). While observing the global response of total occupied receptors in a cell (summed across all the nodes), the temporal fluctuations become almost negligible (Fig 2F). Furthermore, we looked into the heterogeneity in terms of the mean occupied receptor states among nodes for a range of doses of cAMP. We observed that the relative internodal noise in the system decreased with the increasing dose of cAMP and eventually approached a Poisson-like distribution (average Fano factor = 0.85) for the saturating dose of cAMP from the initial sub-Poissonian distribution (S2(D) and S2(E) Fig).
So far, we have studied the intrinsic noise of the system. In a population of cells, the number of receptors and other signaling molecules differ from cell to cell. To study the effect of this cell-to-cell heterogeneity, we varied reaction rates, diffusion constants and the total number of receptors, following a Gaussian distribution with nominal values and standard deviations. In terms of the parameters, we assumed homogeneity in the cell but heterogeneity across the cells. All the nominal values and the standard deviations are listed in Table 1 with further details in S1 File. In this combined study of intrinsic and extrinsic fluctuations, we only observed small differences (std. dev. ∼2.5% at the saturating dose, and ∼10% during initial decay phase following removal of stimulus, n = 10 simulations) indicating a high degree of robustness (Fig 2G).
Lastly, because of the computational burden of dealing with the large number of states and reactions, we considered the possibility that GPCR binding could be represented by a reduced-order model that could capture the essential spatiotemporal dynamics as well as most of the noise characteristics of the full model. To this end, we used a model consisting of one free (R) and a single occupied (R:C) state (S2(F) Fig) and varied the parameters of the reduced model so as to minimize the error between the temporal response of the full and reduced systems. The responses matched closely with only small differences during the initial decay phase following the removal of the stimulus (S2(G) and S2(H) Fig). Importantly, the noise characteristics were also similar over a wide range of cAMP doses; for example, at a saturating dose of cAMP, the standard deviation of the full order was 4.58 molecules compared to 4.59 for the reduced-order model (S2(C) and S2(E) Fig). This finding suggests that there is little error in using the more computationally efficient reduced-order model in lieu of the detailed one.
LEGI module
The occupied receptors: H:C, L:C, S:C, PH:C, PL:C (collectively referred to as RL) pass down the receptor occupancy information to the immediate down-stream signaling network which, follows the LEGI mechanism [13]. Though LEGI successfully explains the temporal and spatial responses to chemoattractant in Dictyostelium [15–17, 30], little attention has been paid to the effect of noise on the LEGI mechanism. Moreover, there are a number of possible ways of implementing LEGI, several of which we considered. In Dictyostelium, chemoattractant signaling depends on the presence of G-proteins. G-proteins form a heterotrimer, consisting of α (we focused on Gα2, which mediates cAMP signaling downstream of the cAR1 receptor [46, 47], β and γ subunits. Whereas these subunits are together in an unoccupied receptor, the α2 and βγ subunits dissociate upon stimulation at which time the latter are free to signal to down-stream elements. This dissociation is persistent [47]. The Gβγ works upstream of Ras and is crucial in chemotaxis [48, 49]. The Ras response, as well as downstream signals, display properties of excitable systems, including a refractory period [24], an all-or-nothing threshold [50], and wave propagation [51–54]. At the same time, near perfect adaptation is observed in the activated Ras response [16]. This suggests that adaptation happens upstream of Ras and, as the response of Gβγ is persistent during cAMP stimulation, we treated this as the excitation process in the LEGI module (Fig 3A). The parameters for the Gβγ dynamics were chosen (Table 2) to match experimentally measured half-times of dissociation (on application of a saturating stimulus) and reassociation (on removal of stimulus) of the Gα2 and Gβγ subunits [17].
(A) Output of receptor module, RL = H:C + L:C + S:C + PH:C + PL:C, drives the LEGI module. LEGI scheme involves a local activator (Gβγ) and global inhibitor (I⋆). Their interaction creates a response regulator (RR) which positively affects the conversion of RasGDP to RasGTP. (B) The regulation of RR can be realized either through a Difference scheme (top) or through a Ratio scheme (bottom). Both involve basal and Gβγ-dependent production of RR. Whereas inactivation is mediated by I⋆ through an intermediate X in the difference scheme, in the ratio scheme it depends on both basal and I⋆-dependent terms. (C) LEGI with Antithetic Integral Feedback (LEGI-AIF). In this scheme, Gβγ and I⋆ create intermediates (X and Y, respectively) that annihilate each other. The RR is created by Y and catalyzes the production of X. (D) Temporal dynamics of components of the different LEGI schemes. The top panel shows nodal average profiles of Gβγ (green) and I⋆ (red) in response to a staircase profile of cAMP stimulus: 0–2 min: 0% R.O.; 2–8 min: 4% R.O.; 8–12 min: 100% R.O. The bottom panels show the corresponding RR profile for the different schemes. The shaded regions denote standard deviations among all nodes from a single simulation. (E) Basal surface profile of Gβγ (green), I⋆ (red) and RR (blue) from different schemes at the time points indicated. (F) Effect of concentrations of cAMP (in terms of % R.O.) on peak amplitude (red), steady-state amplitude (blue), peak time (orange) and adaptation time (green) of the nodal average profile of RR for different schemes. The solid lines and the shaded regions show the respective mean and standard deviations (std. dev. here is the measure of the inter-nodal variation).
Here, we varied the ratio of total number of G-protein (Gα2βγ+Gβγ) to the total number of receptors; for each case, the reaction rate constants of G-protein dynamics were estimated to match the half times of dissociation (t1/2,diss.) and reassociation (t1/2,reass.) from [17]. In reaction no. 2, RL is the total occupied receptors (H:C + L:C + S:C + PH:C + PL:C).
We simulated the system and looked at the dissociated subunits and observed a highly nonlinear relationship between the response curves of receptor occupancy and the dissociated G-protein (S3(A) Fig). The fraction of dissociated G-proteins decreased as the total number of G-proteins increased (S3(A) Fig); however, as a function of the maximal response, it was independent of the total number of G-proteins reaching 50% dissociation at ∼23% R.O. (S3(B) Fig), indicating the existence of “spare receptors” in the system [55]. Additionally, we observed that fluctuations (mean coefficient of variation) in the G-protein response decreased with higher total number of G-proteins (S3(C) Fig). We chose the total G-proteins to be three times the total number of receptors as for the higher values there was no difference in the relative fluctuation level.
Unlike the excitation process, inhibition is G-protein independent in Dictyostelium [17]. To account for this, we assumed an inhibitor existing in both active (I⋆) and inactive (I) forms. The diffusion constants for the inhibitor states were assumed high (a number typical of cytosolic entities [54]) to satisfy LEGI requirements that the global inhibitor have higher diffusivity than the local Gα2βγ and Gβγ. To account for a possible multistep translocation to and from the membrane, we made the inhibitor dynamics slow compared to those of Gβγ.
The interaction of the excitor Gβγ and the activated inhibitor jointly regulate a response regulator, RR, whose action can be modeled using either a difference or a ratio scheme, depending on how the activation and inhibition processes regulate it (Fig 3B and S3(D)–S3(F) Fig). In the difference scheme, basal and Gβγ-dependent activation together with I⋆-mediated deactivation of RR through an intermediate lead to a steady-state RR concentration that is an affine (linear plus a constant) function of the difference between Gβγ and I⋆ (S3(E) Fig and S1 File). An alternative realization, the ratio scheme, involves Gβγ-dependent activation and I⋆-dependent inactivation of RR along with independent basal activation/inactivation and leads to a steady-state of RR that is proportional to the ratio of affine functions of the Gβγ to I⋆ (S3(F) Fig and S1 File). As there are a number of potential biochemical elements that could serve as the response regulator, we considered both schemes. We implemented these schemes explicitly, using the corresponding reactions, as well as implicitly, using a quasi-steady-state approximation of the reactions. By ignoring the transient dynamics, the implicit formulations served to reduce the computational burden of explicitly simulating these systems (S1 File).
One of the possible limitations of the adaptive networks described above is the sensitivity of the adaptation to noise, resulting in variances that do not adapt, but depend on the stimulus level [56]. As a third alternative, we adopted the Antithetic Integral Feedback (AIF) model [57] by making some of the terms local and global, and then compared its performance with the two schemes described above (Fig 3B). Here, Gβγ and I⋆ create intermediates that annihilate each other. One intermediate activates RR which, in turn, catalyzes the production of the other. The detailed reactions and parameter values of all the three LEGI schemes are listed in Table 3.
RL is the total occupied receptors (H:C + L:C + S:C + PH:C + PL:C). The reactions no. 7–15 are for only AIF scheme of LEGI mechanism. The nominal values for the diffusion constants are: , DI = DI* = 20,
. The standard deviations in the diffusion constants were chosen to be 20% of the respective nominal values. Further details are in S1 File.
Simulating the combined GPCR and LEGI modules
To analyze the adaptation characteristics of proposed LEGI networks, we applied spatially uniform increasing concentrations of cAMP starting from no stimulus, (0% R.O.; 0–2 min) rising to a low dose (4% R.O.; 2–8 min) and finally reaching a saturating dose (100% R.O., 8–12 min). As expected, both Gβγ and I⋆ rose with the latter showing a slower response (Fig 3D and 3E and S5(A) and S5(B) Fig). In all variants considered, the initial responses depended on the stimulus level; i.e., the peak amplitude (RRpeak, red) increased and the peak time (Tpeak, orange) decreased with %R.O. (Fig 3F). In all three networks, the mean nodal response regulator activity returned to prestimulus levels (RRss, blue in Fig 3F) but the variance shows stimulus level dependency.
Though the parameter values of the different schemes were chosen to ensure similar mean level of steady-state characteristics (steady-state value, RRss and adaptation time, Tadaptation, S4(A) Fig), we observed differences among the schemes. Across a single cell, the variance over the node increased for the difference scheme, decreased slightly for the ratio scheme, and was smallest and fairly constant for the AIF scheme (blue in Fig 3F). The adaptation times in all three schemes decreased monotonically as a function of %R.O., in agreement with published experimental data [16], with the AIF scheme showing the least sensitivity (green in Fig 3F). Unlike the difference and ratio schemes, the temporal profile of the LEGI-AIF showed some oscillatory behavior during the adaptation at high levels of R.O. The variance in the peak amplitude of the initial response for the ratio scheme was higher than for the difference scheme. The lower coefficient of variation in the peak response profile of the difference scheme provides more certainty of response to change of stimulus than the ratio scheme. Upon removal of the stimulus, the difference scheme returned to the basal level more quickly compared to the ratio scheme (S5(C) and S5(D) Fig).
We next allowed for varying internal concentrations and parameter values so as to capture cellular heterogeneity. We observed the steady-state characteristics (RRss, Tadaptation) were more robust than peak characteristics (RRpeak, Tpeak, S4(B) and S4(C) Fig) to changes in the parameter values. This was true for all the three schemes. For the rest of study, we used the implicit difference scheme of LEGI to avoid unnecessary repetition.
We next simulated the response to a gradient generated by releasing chemoattractant from a micropipette. To this end, we imposed a stationary 3D Gaussian profile centered around an edge point at the basal plane and considered both the apical and basal (Fig 4A) and perimeter responses (Fig 4B and S1 File). Following imposition of the gradient, free Gβγ rose quickly at the front where %R.O. was highest; at the rear the rise was slower (Fig 4C). In contrast, the inhibitor rose more slowly and was fairly uniform around the perimeter of the cell, owing to the relatively high diffusion (Fig 4C and 4D). The resultant response regulator rose sharply at the front and dropped below the basal level at the rear (Fig 4C and 4D). Note that the steady-state level of the response regulator at the front was lower than the peak, as the level of the inhibitor increased, but still displayed levels above basal. To examine the role of diffusion in these patterns, we varied the inhibitor diffusion over a wide range (S6 Fig). Increasing the ratio of inhibitor-to-Gβγ diffusion resulted in greater differences in the response regulator between front and back. The mean difference between RR molecules between front and back doubled for a 10-fold difference in diffusion coefficients.
(A) Schematic showing the diameter at the bottom surface (red) and a semicircular arc on the curved surface (blue) connecting front and back of the cell with respect to the needle position (light blue dot). (B) Circular cross sections of the hemispherical domain at z = 1 (a) and z = 2 μm (b). The front (closest to needle) and back of the cell is marked as 0 and 180°, respectively. (C) Temporal profiles of Gβγ (green), I⋆ (red) and RR (blue) at cell front (0°, darker shade) and back (180°, lighter shade). (D–G) Spatial response of the system for receptor occupancy (R.O.), Gβγ, I⋆ and response regulator (RR). The kymographs (D) are based on the maximal projection of the hemispherical domain (nodes between a and b) of panel B. The white dashed line indicates the time instant when cAMP gradient was applied. Panel E shows the spatial profiles at the basal and apical surfaces at t = 3 min. Panel G shows the spatial profiles along the lines marked in pane A. Lines denote mean and the shaded regions standard deviations (n = 5 independent simulations with parameter fluctuations as in Fig 2G).
At steady state, the linear profiles of local entities such as R.O. and Gβγ displayed gradients that were more diffused on the top surface than on the basal surface (Fig 4E–4G). In contrast, the inhibitor profile was relatively flat along the length of the cell, though random variations were apparent. The resultant response regulator profile showed a gradient and was higher/lower at the front/rear relative to the basal level. We repeated the simulations for a shallower gradient and observed similar behavior, with a smaller degree of localization and slower response in RR (S7 Fig). Additionally, we considered a number of different profiles varying the receptor occupancy at the back (from 12–57%) and front (80–100%). In all cases, the LEGI mechanism was able to sense gradients (S8 Fig).
Signal transduction excitable network (STEN) module
Dictyostelium cells have two excitable systems that work in tandem: a fast cytoskeleton-based network, and a slower signaling network that drives the cytoskeletal system [24, 25, 51]. As we are modeling cells lacking an intact cytoskeleton, we focused on the latter. Our proposed STEN consists of five species: RasGDP, RasGTP, activated and inactivated protein kinase B substrates (PKB*s and PKBs, respectively), and membrane phosphatidylinositol bisphosphate (PIP2), which represents contributions from both PI(3,4)P2 and PI(4,5)P2 (Fig 5A). To reduce the simulation burden, we omitted inactivated protein kinase B substrates (PKBs) from explicit modeling.
(A) Schematic of STEN showing the entities: RasGTP, RasGDP, PIP2 PKBs, PKB*s and their interactions. Green arrows denote positive feedback whereas, orange arrows complete the negative feedback on RasGTP. (B) Temporal profiles of RasGTP (green), PIP2 (red) and PKB*s (blue) at a random node that has fired spontaneously. (C–E) Spatiotemporal profiles of RasGTP (green) and PIP2 (red) of the membrane showing wave-traveling (C), splitting (D) and annihilation (E). (F) Response to a spatially uniform dose of cAMP. Shown are the global RasGTP, PKB*s and PIP2 responses (left), and same at a single random node (center) and the spatiotemporal profile (RasGTP, PIP2) at the basal surface of the cell (right). The yellow arrowheads denote wave initiation sites. Colors are as in panels B–E. The shaded region in the left panel is the standard deviation as in Fig 4D (n = 10). The peaks of the RasGTP and PKB*s profiles correspond to approximately 165,000 and 275,000 molecules. PIP2 peaks at approximately 1,000 molecules/μm2 similar to that reported in [58]. (G) Basal subtracted normalized RR (, left) and RasGTP (right) responses to short (2 s, blue) and long (30 s, red) stimuli. The solid lines and the shaded regions are as in panel F (n = 10). (H) RasGTP response to the two short pulses (2 s) of spatially uniform stimuli with variable delays. Left: temporal mean nodal profile of RasGTP (n = 10). Right: plot of normalized peak of the second response to the first versus the delay between the stimuli.
In our scheme, RasGTP acts as the activator of the excitable system and serves as a front marker of the cell. Recent experiments have demonstrated that lowering PI(4,5)P2 results in increased Ras activity [50]. Similarly, lowering PI(3,4)P2 increased Ras activity through the regulation of RasGAP2 and RapGAP3 [59]. Thus, we incorporated the PIP2-mediated hydrolysis of RasGTP to RasGDP into the model. This closes a positive feedback loop that is formed through mutual inhibition between RasGTP and PIP2.
A slower, negative feedback loop is achieved through the RasGTP-mediated activation of PKBs. This activation is achieved partly by having PH-domain containing PKBA translocate to the membrane to bind to PI(3,4,5)P3, as well as TorC2-mediated phosphorylation of PKBR1 [40, 41]. This negative feedback loop is closed as activated PKBs (PKB*s) negatively regulate RasGTP. There are two proposed mechanisms of negative feedback: 1) controlling the localization of the Sca1/RasGEF/PP2A complex on the plasma membrane through phosphorylation of Sca1; 2) phosphorylation and activation of PI5K which increases PIP2 [41, 60]. As we do not specifically model PI(3,4,5)P3, we implemented this loop with PKBs being activated directly by RasGTP, with a slower time scale to account for the omitted intermediate steps (e.g., PI3K activation, PI(3,4,5)P3 formation and PKB translocation and subsequent phosphorylation) and that some reactions involve cytosolic species, compared to the mutually inhibitory positive feedback loop.
Finally, we coupled this excitable network to the LEGI in two ways. First, through an RR-dependent term that converts RasGDP to RasGTP, consistent with the possibility that the RR is a RasGEF or activates a RasGEF. Second, we also included a term to account for the activation of Phospholipase C (PLC) by Gα2 which leads to PI(4,5)P2 hydrolysis [61]. Table 4 lists the detailed reaction terms and parameter values.
The nominal values for the diffusion constants are: DRasGDP = DRasGTP = DPIP2 = 0.05, . The standard deviations in the diffusion constants were chosen to be 20% of the respective nominal values. The peak values are: RasGTP = 165,000 molec; PKB*s = 275,000 molec; PIP2 = 1,000 molec/μm2 [58]. Further details are in S1 File.
Combining the GPCR, LEGI, and STEN modules
We first simulated our model in the absence of any stimulus. We observed several characteristics of excitable behaviors (Fig 5B–5E), such as traveling waves which were annihilated on collision [62] and occasionally split into smaller waves. When viewing a single node, RasGTP and PIP2 showed mutually exclusive behavior, in which the number of molecules at any one time of one species dominated the other. For example, when PIP2 dominated, there was approximated 40 molecules of PIP2 compared to approximately one molecule of RasGTP (Fig 5B, t = 30 s). The transition to a state in which RasGTP dominated (∼25 molecules of RasGTP vs. ∼0–2 molecules of PIP2) was rapid. In contrast, the number of PKB*s molecules varied considerably less, ranging from ∼40–60 with slow increases happening after the node transitions to a high RasGTP state (Fig 5B, t = 60 s). When viewed across the surface of the cells, we observed wave activity (Fig 5C and S9(A) and S9(B) Fig and S1 Video). The cell was typically in a back state in which high PIP2 levels dominated. However, when stochastic perturbations led to a spot of high RasGTP activity, a wave of activation swept across the cell, moving between the basal and apical surfaces. This eventually extinguished and the cell returned to its basal back state.
Interestingly, after the RasGTP wave went through a region, the subsequent PIP2 recovery was actually higher than before the wave, indicating an overshoot of the basal level (c.f. the PIP2 intensity between the first and last panels in Fig 5C and S9(A) Fig (bottom)). The elevated PIP2 regions played an important role in the steering of waves. Whenever a traveling wave encountered a region of supra-basal PIP2 on its path, it moved around this region which often resulted in the splitting of the wave (Fig 5D and S2 Video). Generally, wave splits were rare because of the small size of the cell relative to the wave. This is consistent with experimental findings which have prompted experimentalists studying waves to consider giant fused cells [51, 53]. Consistent with properties of excitable waves, when two wave fronts met, they annihilated (Fig 5E and S3 Video). The annihilation time for wave collisions depended on the time scale and relative diffusivity of RasGTP and PKB*s.
Response to stimulation
We simulated the response of our model to a spatially uniform stimulation using saturating dose of cAMP (100% R.O., Fig 5F and S4 Video). After an initial delay of 2–3 seconds, we observed multiple wave initiations (arrowheads in Fig 5F) and wave spreading. The waves eventually spanned the whole cell at which time (∼8–9 s) total Ras activity peaked. The activity died down as the RR level returned to the pre-stimulus condition. The experimentally observed response of cells is quite similar following short or long pulses of chemoattractant stimulation, consistent with the notion that the underlying signaling network is excitable [7, 24]. To test this in our model, we repeatedly applied short (2 s) and long (30 s) global stimuli to models with varying parameter values and compared the respective RR and RasGTP responses (Fig 5G). The response profiles were nearly indistinguishable, with both peaking about 7–8 s after the application of the stimulus, followed by a return to the pre-stimulus level. They only differed in the final phase in which the response to the short stimulus reached its steady-state faster (∼20 vs. 30 s). Interestingly, the RasGTP response was only noticeable after ∼2 s at which time the short stimulus had already been removed. This is consistent with excitable systems which reach a point of no return following stimulation. When the stimulus remained present beyond the time taken for the excitable system to return to its basal state (> 60 s), we observed a smaller and more patchy second wave of activity (S9(C)–S9(F) Fig). This second peak of the response is due to the partially adapted state of the LEGI module which is due to the adaptation time being longer than the duration of a pulse from the excitable system; it has been seen experimentally both with signaling [63] and cytoskeletal biosensors [64].
Cells display refractory periods following stimulation during which further excitation fails to trigger a response, or diminished in intensity [7, 24]. We applied a series of double pulses, each of duration 2 s with variable delays (10–90 s) and compared the peak of the second response to the first response (Fig 5H, left panel). For short delays (<10 s) there was no distinct second response. However, as we increased the time delay, we observed an increase in the peak amplitude of the second response. After a sufficiently long delay (80 s), the second peak almost matched the first peak with a 50% recovery occurring at ∼50 s. (Fig 5H, right panel).
When stimulated by a gradient, cells display persistent high levels of activity at the side of the cell facing the chemoattractant source. We simulated these experiments by introducing a gradient of receptor occupancy across the cell. RasGTP showed a localized persistent patch of high activity facing the side with highest receptor occupancy (Fig 6A and S5 Video). We also simulated two experiments in which the spatial gradient was combined with a temporal stimulus. In the first, we introduced a gradient and waited until the cell displayed a spatial response; we then removed it for a variable time, before finally reintroducing it [65]. The crescent disappeared following removal of the stimulus but did not return to its full strength until the delay was ∼60 s (Fig 6B and S6 Video), matching our previous observations in the double pulse experiment. In the second simulation, we applied a large global stimulus following the establishment of a crescent in response to a gradient, and then removed all cAMP. In this case, the global stimulus elicited a response everywhere (Fig 6C and S7 Video). These simulations show the complex interactions of spatial and temporal components of coupled LEGI and STEN systems.
(A–C) Shown are the temporal profiles of RasGTP (green), PIP2 (red), PKBs*(blue), Gα2 (teal), Gβγ(orange) and basal subtracted normalized RR (, cyan) at single nodes at the front and back of the cell. The kymographs (right) show RasGTP (green) and PIP2 (red). Solid white line denotes when the gradient was applied, and the dashed line shows the needle position (front). Whereas Panel A shows the response to a single gradient, B and C show the response to two gradient stimuli with a delay of 60 s (B) and to a gradient stimulation followed by a global one (C).
Effect of lowering threshold
Recent experiments in which the activity and motility of the cell were altered suggest that these came about through changes in the threshold of the excitable signaling system. Our model allows to test some of these perturbations. We first considered lowering PIP2 levels by adding an extra degradation term, so as to recreate experiments in which the phosphatase Inp54p is brought synthetically to the cell membrane [7, 51]. In this case, we saw elevated levels of activity and the cell commenced periodic whole-cell increases in activity similar to what have been observed experimentally (Fig 7A). Similar, though less acute increases in activity were seen in simulations in which PKB*s levels in the cell were lowered through a reduction in the RasGTP-mediated activation rate (Fig 7B). In this case, bursts of wave initiations were observed, but the resulting waves did not cover the whole cell surface.
(A,B) The kymographs show the effect of lowering the STEN threshold by inhibiting PIP2(A) and PKB*s (B) as denoted in the schematics. The white lines indicate the time at which the respective species were lowered. (C) Schematic for incorporating the differential threshold between top and bottom surface of the cell through altering the PKB*s mediated inhibition on RasGTP. (D) Increased threshold restricts the wave activities at the basal surface and the waves were not allowed to travel to the apical surface. (E) Higher threshold on the apical surface made small waves with fewer number of wave initiations.
Experiments have also demonstrated that mechanical contacts can trigger excitable behavior [66–68]. Finally, we considered the possibility that the threshold was differentially regulated by mechanical contact with the substrate, with a lower threshold at the basal surface than at the apical surface. To take this into account, we assumed that the PKBs-mediated inhibition of RasGTP was different between the basal and apical surfaces, with lower inhibition at the former (Fig 7C). In these simulations, we observed frequent wave initiations at the basal surface, but these waves were quickly extinguished when they reached the apical surface (Fig 7D). We rarely saw any de novo waves at the apical surface, and when they did appear, they did not spread like the basal counterpart, and were extinguished rapidly (Fig 7E).
Discussion
Mathematical models have been a popular means of understanding how cells move in response to chemoattractant stimuli. These models have taken a modular view of the overall signaling network, breaking the overall system down into simpler functional blocks [69]. In doing so, most models have focused on a specific set of experimental observations, such as receptor-ligand interactions [42], adaptation [13, 16, 70], amplification [71], wave propagation [72–75], excitability [25, 54, 76], or polarity [77–79], by concentrating on a specific functional block sometimes ignoring its connection to the overall network. Part of this problem is a lack of experimental data as most experimental papers focus on specific aspects of the chemotactic behavior. Here, we have proposed an integrated model that takes the receptor-ligand binding to PKBs, which have been shown to be a link to the cytoskeletal network. Where the experiments do not offer single out a unique alternative, we have presented various possibilities and simulated scenarios that could be used to distinguish among them, for example, we presented different schemes of LEGI.
Because of the importance of noise-induced transitions in triggering the signal transduction excitable network, thereby allowing unstimulated cells to move randomly, our modeling approach has emphasized the role of noise. While the role of stochastic fluctuations in patterning cellular responses is now firmly appreciated, most of these studies consider molecules and species that exist in relatively small numbers [31, 32]. As our simulations here demonstrate, even in signaling systems that consist of biochemical species with high copy numbers, noise can play an important role in regulating cell physiology.
There are a number of approaches for the stochastic simulation of biological reactions [80, 81]. The chemical Langevin equation (CLE) is one of the most widely used as it is easy to obtain from the corresponding deterministic ODE model. However, there are some limitations. The moments of CLE and those of the Chemical Fokker-Planck Equation (CFPE) are only equivalent for unimolecular reaction systems [82] but are not the same for highly nonlinear systems like those considered here [83]. Furthermore, the Gaussian noise approximation inherent in the CLE does not restrict the values of the states to be non-negative which can lead to the physically implausible results [84, 85]. This problem becomes more acute when dealing with low copy numbers. There are simulation methods developed to avoid the limitations of CLE, but most often they are found to introduce some additional errors [86]. In presence of diffusion, the stochastic ordinary differential equations of the CLE framework become partial differential equations (PDEs). Solving these stochastic PDEs, particularly on unstructured meshes is difficult.
The URDME approach used here has several advantages over the stochastic PDE approach. The stochastic framework that we used to simulate our model does not require an artificial injection of noise. Instead, noise is a natural consequence of the stochastic description of the reaction-diffusion equations. Moreover, highly complex geometries can be also be considered. Nevertheless, there are some drawbacks. For example, quasi-steady-state assumptions that are usually made (e.g. reaction 3 in Table 4) can introduce errors in stochastic simulations even in cases when the approximation would be appropriate in a deterministic setting [87, 88]. The URDME approach is also computationally expensive, particularly as the copy numbers of the species increases. We should point out that other approaches exist all with their advantages and disadvantages; see, for example, the description in [89, 90].
In this stochastic setting, the LEGI mechanism adapts perfectly as in the deterministic setting, as long as we focus on the mean level of activity (Fig 3D and S5(A) and S5(B) Fig). However, we showed that the variance in the response increases as a function of the stimulus strength (Fig 3F) and that the size of this increase depends on the particular way that the LEGI is implemented. The increase in the variance observed could account for chemokinesis, the increased speed of random migration seen following a uniform chemokine stimulus in neutrophils [91]. Moreover, because the different implementations of the LEGI mechanism have distinctive noise patterns, our results suggest a way of elucidating the precise nature of the adaptation mechanism experimentally.
Our models define the essential characteristics of the LEGI inhibitor and facilitate its identification. This has been difficult as there are inhibitors in each of the modules. For example, STEN possesses an inhibitor that accounts for fast shutoff of several biosensors (e.g. RBD). This inhibition is mediated by PKBs [51, 92]. However, we have shown experimentally that this STEN-inhibition is separate from the LEGI inhibitor [17]. There is also inhibition that contributes to the polarization of cells and comes from the cytoskeleton [51] which we did not model here. We and others [93–95] have suggested that this may be related to mechanical properties of the cell. However, while both these processes provide inhibition, they are not the I that we propose in this paper, which is specific to the LEGI mechanism. Importantly, this LEGI inhibitor would be easy to miss by solely looking at cell migration, as cells lacking it can chemotax, albeit it with lower efficiency than cells with it [22]. In fact, certain cells like fibroblasts do not have robust adaptation, but can chemotax though they respond only to a narrow range of chemoattractant concentrations [96]. Nevertheless, there are a several facts that we can state about the inhibitor we seek. First, experimental data of Gβγ-null cells in presence of cAMP substantiates the existence of an inhibitor in LEGI. In Gβγ-null cells, basal activity (PH or LimE) decreases over one or two minutes and remains down as long as the stimulus is present [17]. On removal of the stimulus, the activity returns within one or two minutes. This demonstrates that there is a cAMP-mediated, but Gβγ-independent persistent inhibition (consistent with the role that we assign to I* in this paper). We think that this is the best assay to identify the inhibitor. One other aspect of this inhibitor is that it is not a feedback inhibitor. Comparing model predictions of feedback and feedforward inhibition with experiments in which multiple stimuli are applied in series [17] demonstrated that the inhibition does not use feedback. A similar conclusion was also reached by considering the response to multiple chemoattractant doses [16].
The biased-excitable network hypothesis suggests that chemoattractant stimulation leads to a signal that lowers the threshold of activation of the STEN, and various models and experiments support this view. However, how cells display a persistent crescent of front-markers towards the source of a chemoattractant gradient remains an open question [97]. Simulations of cells stimulated by gradients of chemoattractant show that the RR activity remains high at the cell front and is suppressed elsewhere. One would expect that the resultant spatially localized lowered threshold would trigger activity that would travel away from their point of origin. Further activity would be possible, after the delay imposed by the refractory period. To explain this, we modified the STEN dynamics such that, following stimulation, the system undergoes an excitable-to-bistable transition (S9(G) and S9(H) Fig). In this case, a persistent, elevated level of RR leads to a new “high” state of activity in the STEN, and persistent crescents were observed. Though our implementation is similar to the wave pinning scheme used to explain persistent polarization [98], it differs due to the fact that our system is only in the bistable region in the presence of a persistent high stimulus.
Our model recreates most of the observed responses, including adaptation to persistent global stimuli (Fig 5C–5G), presence of secondary peaks (S9(C)–S9(F) Fig, [99]), spatially localized responses to gradients (Fig 5H–5I) along with the typical excitable behaviors such as wave annihilation and refractory period. Moreover, our simulations show the power of having an excitable system at the heart of the chemotactic signaling system. This threshold becomes the most important parameter shaping the cellular response, dictating overall activity and eventually, migratory modes [50]. Increased number of wave initiations as well as more wave spreading is often associated with lowering of the threshold, whereas increase in the threshold corresponds to opposite effect. The threshold in the chemotactic signaling network can be altered in several ways and not all the alterations have the same effect on the system behavior. Theoretically weakening all the inhibitory pathways or strengthening the catalytic pathways acting on RasGTP can result in lowering of the threshold. Experimentally, several pharmacological/genetic perturbations have been used to study this effect of alteration in Dictyostelium cells [50, 51, 100]. Recent publications suggested of possible existence of a number of mechanical feedbacks on the leading edge protrusion [66, 101] and hence on the system threshold. Cao et al. reported that the basal surface waves are mostly restricted at the bottom surface unless the threshold of the system is lowered [67]. In that case, on reaching the bottom boundary the basal surface waves starts traveling upwards long the apical surface.
In our study, we only have considered the mechanical input to the signaling network through mechanical contacts, as suggested in [67]. Others factors, like cell geometry, membrane curvature and cortical tension are also likely to influence signaling waves. In particular, wave pinning could be caused by the localization of signaling molecules in cell structures that represent significantly different surface-to-volume characteristics, like filopodia, dendritic spines, and primary cillia [102]. The resultant changes in cortical tension could also affect polarization in cells [37, 94, 95]. Change in the membrane curvature also recruit curvature-sensing proteins such as Bar-domain proteins which are known to interact with the signaling network at various levels [103, 104]. In this present study, we did not include any cytoskeletal proteins. Moreover, our solution domain is fixed and of regular geometry; we leave these research avenues for the future study. We note that the URDME approach was recently combined with a moving boundary framework to describe yeast polarization [105].
Finally, while our study has centered on chemoattractant signaling, we note that various components, such as GPCR signaling, the presence of feedforward adaptive mechanisms, and biochemical excitability are concepts that extend well beyond the realm of cell migration. For example, a recent report suggests that stochastic effects may play an important role in rhodopsin signaling [106]. Thus, our study provides a means for analyzing the effect of noise in these systems and could help to understand more complex interactions in these other settings.
Supporting information
S1 Fig. Effect of mesh refinement.
(A) Total number of nodes (in logarithmic scale) for different values of Hmax (maximum allowable distance between nodes). Hmin (minimum allowable distance between nodes) was chosen to be 1/3 of the respective Hmax value. The red dot shows the mesh size adopted for the present study whereas the gray dots are used for the comparison in panel (C). (B) Simulation time (in logarithmic scale) for different values of Hmax corresponds to a one second simulation of a diffusion process involving a single entity. Colors are as in (A). (C) Comparison of the simulation output (basal surface profile) of a diffusion process involving single entity at t = 1 s using fine, medium and coarse meshes as indicated by the Hmax values. (D) Simulation time (in logarithmic scale) for different values of Hmax corresponding to a one second simulation of GPCR and LEGI modules combined, in the absence of cAMP. Colors are same as in (A).
https://doi.org/10.1371/journal.pcbi.1008803.s001
(EPS)
S2 Fig. Temporal response of different full- and reduced-order receptor states.
(A) Temporal responses of free (H,L,S), occupied (H:C,L:C,S:C), free phosphorylated (PH,PL) and occupied phosphorylated (PH:C,PL:C) states of the receptors, in the absence and presence of a saturating dose of cAMP. The gray boxes denote the time segments used to compute steady-state concentration of different entities in (B). (B) Relative steady-state amount of different receptor states in absence (top) and presence (bottom) of cAMP. (C) Mean Fano factor of the temporal fluctuations among nodes for full-(blue solid) and reduced-order (red dashed) modules in response to varying cAMP levels (% R.O.). (D) Plot of total number of occupied receptors at nodes. Solid line and shaded region show the mean and the standard deviation respectively. (E) Comparison between internode noise characteristics of full- (solid) and reduced-order (dashed) modules in terms of the coefficient of variation (blue) and Fano factor (red). (F) Schematic of a two-state reduced order model of receptors containing single unoccupied (R) and occupied (R:C) states. (G,H) Comparison between global responses of full- (blue solid) and reduced-order (red dashed) receptor modules in presence and absence of cAMP. The region where the two responses differ is denoted by the gray box in (G), which is zoomed into in (H).
https://doi.org/10.1371/journal.pcbi.1008803.s002
(EPS)
S3 Fig. G-Protein response and LEGI signaling.
(A) Dissociated G-protein response (%) to varying degree of receptor occupancy (% R.O.) with different total number of G-protein molecules (represented by different colors and denoted as a ratio to the total receptors molecules). (B) G-protein responses of (A) normalized to respective maximum values. The inset shows a zoomed-in version of a smaller section of the plot. Colors are same as in (A). (C) Coefficient of variation as a function of %R.O. Colors are same as in (A). (D) Schematic of implicit LEGI mechanisms. (E,F) Comparison between implicit (left) and explicit (right) schemes of difference (C) and ratio (bottom) mechanism. Whereas the explicit mechanism uses detailed reactions involving Gβγ, I* and RR, the implicit schemes use the steady-state expression of RR directly in the model ignoring the transient dynamics. (G) Effect of concentrations of cAMP (in terms of % R.O.) on steady-state amplitude (blue), peak amplitude (red), adaptation time (green) and peak time (orange) of the nodal average profile of RR in the ratio mechanism for different basal activation/inactivation rate. The solid lines and the shaded regions show the respective mean and standard deviations (n = 10 independent simulations in which the parameter values were varied according to the distributions of Table 3).
https://doi.org/10.1371/journal.pcbi.1008803.s003
(EPS)
S4 Fig. Comparison among noise levels in LEGI schemes.
(A) Effect of concentrations of cAMP (in terms of % R.O.) on steady-state amplitude (top panel; blue solid: difference, light blue dashed: ratio), peak amplitude (top panel; red solid: difference, magenta dashed: ratio), coefficient of variation (COV) of steady-state (middle panel; blue solid: difference, light blue dashed: ratio) and peak amplitude (middle panel; red solid: difference, magenta dashed: ratio), adaptation time (bottom panel; green solid: difference, light green dashed: ratio) and peak time (bottom panel; orange solid: difference, orange dashed: ratio) of the nodal average profile of RR. Nominal values of the parameters were used. (B,C) Effect of parameter fluctuations on the profiles in (A) where parameters were varied according to the distributions of Table 3. The solid/dashed lines and the shaded regions show the respective mean and standard deviations (n = 10). The colors are same as in (A).
https://doi.org/10.1371/journal.pcbi.1008803.s004
(EPS)
S5 Fig. Temporal response of LEGI mechanisms.
(A,B) Temporal nodal profiles of Gβγ(green), I⋆ (red) and RR (blue) in response to a staircase profile of cAMP for the difference (A) and ratio (B) mechanisms. (C,D) Temporal average (C) and individual (D) nodal profiles of Gβγ (green), I⋆ (red) and RR (blue) in response to the application and withdrawal of stimulus for difference (center) and ratio (bottom) mechanisms. The shaded regions in (C) denote standard deviations among all nodes from a single simulation.
https://doi.org/10.1371/journal.pcbi.1008803.s005
(EPS)
S6 Fig. Effect of diffusion of LEGI inhibitor during gradient stimulation.
(A,B,C) Temporal difference (between front and back of the cell) profiles of Gβγ (A), I⋆ (B) and response regulator (C, scaled) for different relative diffusive strengths of LEGI inhibitor, I⋆. The gray boxes denote the time segments used to generate the average values of the individual profiles for (D). (D) Plot showing steady-state temporal average values of the difference profiles Gβγ (green), I⋆ (red) and response regulator (blue) for different relative diffusive strengths of I⋆ (in logarithmic scale). Mean±SEM is shown (n = 5 independent simulations in which the parameter values were varied according to the distributions of Table 3).
https://doi.org/10.1371/journal.pcbi.1008803.s006
(EPS)
S7 Fig. Shallow gradient sensing by LEGI difference scheme.
(A) Schematic showing the diameter at the bottom surface (red) and a semicircular arc on the curved surface (blue) connecting front and back of the cell with respect to the needle position (light blue dot). (B) Circular cross sections of the hemispherical domain at z = 1 (a) and z = 2 μm (b). The front (closest to needle) and back of the cell is marked as 0 and 180°, respectively. (C) Temporal profile of Gβγ (green), I⋆ (red) and RR (blue) at cell front (0°, darker shade) and back (180°, lighter shade). (D-G) Spatial response of the system for receptor occupancy (R.O.), Gβγ, I⋆ and response regulator (RR). The kymographs (D) are based on the maximal projection of the hemispherical domain (nodes between a and b) of panel B. The white dashed line indicates the time instant when cAMP gradient was applied. Panel E shows the spatial profiles at the basal and apical surfaces at t = 3 min. Panel G shows the spatial profiles along the lines marked in panel A. Lines denote mean and the shaded regions standard deviations (n = 5 independent simulations in which the parameter values were varied according to the distributions of Table 3).
https://doi.org/10.1371/journal.pcbi.1008803.s007
(EPS)
S8 Fig. Gradient sensing by LEGI difference scheme for varying receptor occupancy at the front and the back.
(A) Box plot showing the steady-state temporal distribution of the receptor occupancy at the front (0°, black) and the back (180°, gray) of the cell. (B,C) Temporal profile of Gβγ (green), I⋆ (red) and RR (blue) at cell front (0°, darker shade) and back (180°, lighter shade). The solid lines and the shaded regions show the respective mean and standard deviations (n = 10 independent simulations in which the parameter values were allowed to vary according to the distributions of Table 3). (D) Pre-(cyan) and post-stimulus (blue) profiles of RR. The profiles are obtained by temporal averaging over a span of 5 s as shown by the rectangular boxes of respective colors in (C). These simulations used nominal parameter values in the LEGI difference model.
https://doi.org/10.1371/journal.pcbi.1008803.s008
(EPS)
S9 Fig. STEN dynamics.
(A) Spatiotemporal profiles of RasGTP-PKB*s (green-blue, top) and PIP2-PKB*s (red-blue, bottom) on the basal surface showing wave-traveling. The regions bounded by white dashed lines show the membrane regions with high RasGTP. (B) The Color-coded overlays show the progression of waves as a function of time and computation of the wave speed. (C) The nodal responses of RasGTP (green), PIP2 (red) and PKB*s (blue) at three random membrane nodes in presence of a global cAMP stimulus. (D) Basal surface profile of RasGTP (green) and PIP2 (red) in presence of a global stimulus at the time points indicated. (E,F) Normalized global response of RasGTP (D) and the corresponding kymograph in response to a sustained, spatially uniform cAMP dose. (G) Illustration showing how an excitable-to-bistable transition could occur. Shown are hypothetical nullclines for RasGTP (green) and PKB*s (blue) for a deterministic, two-state model for the excitable system. The darker and the lighter green shade represent the corresponding RasGTP nullcline in the excitable and bistable regimes, respectively. (H) Illustration showing the effect of a gradient stimulus. The dashed line represents the RasGTP-nullcline in the absence of a stimulus. The darker and the lighter green shade represent the RasGTP nullclines at the front and back of the cell, respectively, when subjected to cAMP gradient. The scale bars in (A,B,D) represents 2 μm.
https://doi.org/10.1371/journal.pcbi.1008803.s009
(EPS)
S2 File. MATLAB implementation.
Matlab scripts used for simulating the various models.
https://doi.org/10.1371/journal.pcbi.1008803.s011
(ZIP)
S1 Video. Wave traveling.
Video showing various views of the wave traveling. Corresponds to Fig 5C.
https://doi.org/10.1371/journal.pcbi.1008803.s012
(MP4)
S2 Video. Wave splitting.
Video showing various views of the wave splitting. Corresponds to Fig 5D.
https://doi.org/10.1371/journal.pcbi.1008803.s013
(MP4)
S3 Video. Wave annihilation.
Video showing various views of the wave annihilation. Corresponds to Fig 5E.
https://doi.org/10.1371/journal.pcbi.1008803.s014
(MP4)
S4 Video. Global stimulus.
Video showing various views of the response to a global chemoattractant stimulus. Corresponds to Fig 5F.
https://doi.org/10.1371/journal.pcbi.1008803.s015
(MP4)
S5 Video. Gradient stimulus.
Video showing various views of the response to a gradient chemoattractant stimulus. Corresponds to Fig 6A.
https://doi.org/10.1371/journal.pcbi.1008803.s016
(MP4)
S6 Video. Response to a successive application of gradient stimulii.
Corresponds to Fig 6B.
https://doi.org/10.1371/journal.pcbi.1008803.s017
(MP4)
S7 Video. Response to a gradient stimulation followed by a global stimulation.
Corresponds to Fig 6C.
https://doi.org/10.1371/journal.pcbi.1008803.s018
(MP4)
Acknowledgments
The authors would like to thank all members of the P.A. Iglesias (Johns Hopkins University), P.N. Devreotes and D.N. Robinson laboratories (Johns Hopkins University School of Medicine) for their valuable contributions in various points of the research process.
References
- 1. Richardson BE, Lehmann R. Mechanisms guiding primordial germ cell migration: strategies from different organisms. Nat Rev Mol Cell Biol. 2010;11(1):37–49.
- 2. Norden C, Lecaudey V. Collective cell migration: General themes and new paradigms. Curr Opin Genet Dev. 2019;57:54–60.
- 3. Nourshargh S, Alon R. Leukocyte migration into inflamed tissues. Immunity. 2014;41(5):694–707.
- 4. Michael M, Vermeren S. A neutrophil-centric view of chemotaxis. Essays Biochem. 2019;63(5):607–618.
- 5. Friedl P, Alexander S. Cancer invasion and the microenvironment: plasticity and reciprocity. Cell. 2011;147(5):992–1009.
- 6. Bravo-Cordero JJ, Hodgson L, Condeelis J. Directed cell invasion and migration during metastasis. Curr Opin Cell Biol. 2012;24(2):277–83.
- 7. Zhan H, Bhattacharya S, Cai H, Iglesias PA, Huang CH, Devreotes PN. An excitable Ras/PI3K/ERK signaling network controls migration and oncogenic transformation in epithelial cells. Dev Cell. 2020;54(5):608–623.e5.
- 8. Iglesias PA, Devreotes PN. Navigating through models of chemotaxis. Curr Opin Cell Biol. 2008;20(1):35–40.
- 9. Holmes WR, Carlsson AE, Edelstein-Keshet L. Regimes of wave type patterning driven by refractory actin feedback: Transition from static polarization to dynamic wave behaviour. Phys Biol. 2012;9(4):046005.
- 10. Rappel WJ, Edelstein-Keshet L. Mechanisms of cell polarization. Curr Opin Syst Biol. 2017;3:43–53.
- 11. Koshland DE Jr. A response regulator model in a simple sensory system. Science. 1977;196(4294):1055–63.
- 12. Parent CA, Devreotes PN. A cell’s sense of direction. Science. 1999;284(5415):765–70.
- 13. Levchenko A, Iglesias PA. Models of eukaryotic gradient sensing: Application to chemotaxis of amoebae and neutrophils. Biophys J. 2002;82(1 Pt 1):50–63.
- 14. Ma W, Trusina A, El-Samad H, Lim WA, Tang C. Defining network topologies that can achieve biochemical adaptation. Cell. 2009;138(4):760–73.
- 15. Janetopoulos C, Ma L, Devreotes PN, Iglesias PA. Chemoattractant-induced phosphatidylinositol 3,4,5-trisphosphate accumulation is spatially amplified and adapts, independent of the actin cytoskeleton. Proc Natl Acad Sci U S A. 2004;101(24):8951–6.
- 16. Takeda K, Shao D, Adler M, Charest PG, Loomis WF, Levine H, et al. Incoherent feedforward control governs adaptation of activated Ras in a eukaryotic chemotaxis pathway. Sci Signal. 2012;5(205):ra2. pmid:22215733
- 17. Tang M, Wang M, Shi C, Iglesias PA, Devreotes PN, Huang CH. Evolutionarily conserved coupling of adaptive and excitable networks mediates eukaryotic chemotaxis. Nat Commun. 2014;5:5175.
- 18. Meinhardt H, de Boer PA. Pattern formation in Escherichia coli: A model for the pole-to-pole oscillations of Min proteins and the localization of the division site. Proc Natl Acad Sci U S A. 2001;98(25):14202–14207.
- 19. Xiong Y, Huang CH, Iglesias PA, Devreotes PN. Cells navigate with a local-excitation, global-inhibition-biased excitable network. Proc Natl Acad Sci U S A. 2010;107(40):17079–17086.
- 20. Neilson MP, Veltman DM, van Haastert PJM, Webb SD, Mackenzie JA, Insall RH. Chemotaxis: A feedback-based computational model robustly predicts multiple aspects of real cell behaviour. PLoS Biol. 2011;9(5):e1000618.
- 21. Shao D, Levine H, Rappel WJ. Coupling actin flow, adhesion, and morphology in a computational cell motility model. Proc Natl Acad Sci U S A. 2012;109(18):6851–6.
- 22. Shi C, Huang CH, Devreotes PN, Iglesias PA. Interaction of motility, directional sensing, and polarity modules recreates the behaviors of chemotaxing cells. PLoS Comput Biol. 2013;9(7):e1003122.
- 23. Alonso S, Stange M, Beta C. Modeling random crawling, membrane deformation and intracellular polarity of motile amoeboid cells. PLoS One. 2018;13(8):e0201977.
- 24. Huang CH, Tang M, Shi C, Iglesias PA, Devreotes PN. An excitable signal integrator couples to an idling cytoskeletal oscillator to drive cell migration. Nat Cell Biol. 2013;15(11):1307–1316.
- 25. van Haastert PJM, Keizer-Gunnink I, Kortholt A. Coupled excitable Ras and F-actin activation mediates spontaneous pseudopod formation and directed cell movement. Mol Biol Cell. 2017;28(7):922–934.
- 26. Cooper RM, Wingreen NS, Cox EC. An excitable cortex and memory model successfully predicts new pseudopod dynamics. PLoS One. 2012;7(3):e33528.
- 27. FitzHugh R. Impulses and physiological states in theoretical models of nerve membrane. Biophys J. 1961;1(6):445–466.
- 28. Hecht I, Kessler DA, Levine H. Transient localized patterns in noise-driven reaction-diffusion systems. Phys Rev Lett. 2010;104(15):158301.
- 29. Flemming S, Font F, Alonso S, Beta C. How cortical waves drive fission of motile cells. Proc Natl Acad Sci U S A. 2020;117(12):6330–6338.
- 30. Ma L, Janetopoulos C, Yang L, Devreotes PN, Iglesias PA. Two complementary, local excitation, global inhibition mechanisms acting in parallel can explain the chemoattractant-induced regulation of PI(3,4,5)P3 response in dictyostelium cells. Biophys J. 2004;87(6):3764–74.
- 31. Rao CV, Wolf DM, Arkin AP. Control, exploitation and tolerance of intracellular noise. Nature. 2002;420(6912):231–7.
- 32. Balázsi G, van Oudenaarden A, Collins JJ. Cellular decision making and biological noise: from microbes to mammals. Cell. 2011;144(6):910–25.
- 33. Drawert B, Engblom S, Hellander A. URDME: A modular framework for stochastic simulation of reaction-transport processes in complex geometries. BMC Syst Biol. 2012;6(1):76.
- 34. Fange D, Elf J. Noise-induced Min phenotypes in E. coli. PLoS Comput Biol. 2006;2(6):e80.
- 35. Gibson MA, Bruck J. Efficient exact stochastic simulation of chemical systems with many species and many channels. J Phys Chem A. 2000;104(9):1876–1889.
- 36. Drawert B, Hellander S, Trogdon M, Yi TM, Petzold L. A framework for discrete stochastic simulation on 3D moving boundary domains. J Chem Phys. 2016;145(18):184113.
- 37. Trogdon M, Drawert B, Gomez C, Banavar SP, Yi TM, Campàs O, et al. The effect of cell geometry on polarization in budding yeast. PLoS Comp Biol. 2018;14(6):e1006241. pmid:29889845
- 38. Gillespie DT. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J Comp Phys. 1976;22(4):403–434.
- 39. Ditlev JA, Vacanti NM, Novak IL, Loew LM. An open model of actin dendritic nucleation. Biophys J. 2009;96(9):3529–42.
- 40. Kamimura Y, Xiong Y, Iglesias PA, Hoeller O, Bolourani P, Devreotes PN. PIP3-independent activation of TorC2 and PKB at the cell’s leading edge mediates chemotaxis. Curr Biol. 2008;18(14):1034–1043.
- 41. Charest PG, Shen Z, Lakoduk A, Sasaki AT, Briggs SP, Firtel RA. A Ras signaling complex controls the RasC-TORC2 pathway and directed cell migration. Dev Cell. 2010;18(5):737–749.
- 42. Van Haastert PJ, De Wit RJ. Demonstration of receptor heterogeneity and affinity modulation by nonequilibrium binding experiments. The cell surface cAMP receptor of Dictyostelium discoideum. J Biol Chem. 1984;259(21):13321–8.
- 43. Caterina MJ, Devreotes PN, Borleis J, Hereld D. Agonist-induced loss of ligand binding is correlated with phosphorylation of cAR1, a G protein-coupled chemoattractant receptor from Dictyostelium. J Biol Chem. 1995;270(15):8667–72.
- 44. Caterina MJ, Hereld D, Devreotes PN. Occupancy of the Dictyostelium cAMP receptor, cAR1, induces a reduction in affinity which depends upon COOH-terminal serine residues. J Biol Chem. 1995;270(9):4418–23.
- 45. Ueda M, Sako Y, Tanaka T, Devreotes P, Yanagida T. Single-molecule analysis of chemotactic signaling in Dictyostelium cells. Science. 2001;294(5543):864–7.
- 46. Okaichi K, Cubitt AB, Pitt GS, Firtel RA. Amino acid substitutions in the Dictyostelium Gα subunit Gα2 produce dominant negative phenotypes and inhibit the activation of adenylyl cyclase, guanylyl cyclase, and phospholipase C. Mol Biol Cell. 1992;3(7):735–47.
- 47. Janetopoulos C, Jin T, Devreotes P. Receptor-mediated activation of heterotrimeric G-proteins in living cells. Science. 2001;291(5512):2408–11.
- 48. Jin T, Zhang N, Long Y, Parent CA, Devreotes PN. Localization of the G protein βγ complex in living cells during chemotaxis. Science. 2000;287(5455):1034–6.
- 49. Lehmann DM, Seneviratne AMPB, Smrcka AV. Small molecule disruption of G protein βγ subunit signaling inhibits neutrophil chemotaxis and inflammation. Mol Pharmacol. 2008;73(2):410–8.
- 50. Miao Y, Bhattacharya S, Edwards M, Cai H, Inoue T, Iglesias PA, et al. Altering the threshold of an excitable signal transduction network changes cell migratory modes. Nat Cell Biol. 2017;19(4):329–340. pmid:28346441
- 51. Miao Y, Bhattacharya S, Banerjee T, Abubaker-Sharif B, Long Y, Inoue T, et al. Wave patterns organize cellular protrusions and control cortical dynamics. Mol Syst Biol. 2019;15(3):e8585. pmid:30858181
- 52. Bretschneider T, Anderson K, Ecke M, Müller-Taubenberger A, Schroth-Diez B, Ishikawa-Ankerhold HC, et al. The three-dimensional dynamics of actin waves, a model of cytoskeletal self-organization. Biophys J. 2009;96(7):2888–2900. pmid:19348770
- 53. Gerhardt M, Walz M, Beta C. Signaling in chemotactic amoebae remains spatially confined to stimulated membrane regions. J Cell Sci. 2014;127(23):5115–25.
- 54. Fukushima S, Matsuoka S, Ueda M. Excitable dynamics of Ras triggers spontaneous symmetry breaking of PIP3 signaling in motile cells. J Cell Sci. 2019;132(5). pmid:30745337
- 55. Stephenson RP. A modification of receptor theory. Br J Pharmacol Chemother. 1956;11(4):379–93.
- 56. Iglesias PA, Shi C. Comparison of adaptation motifs: Temporal, stochastic and spatial responses. IET Syst Biol. 2014;8(6):268–281.
- 57. Briat C, Gupta A, Khammash M. Antithetic integral feedback ensures robust perfect adaptation in noisy biomolecular networks. Cell Syst. 2016;2(2):133.
- 58. Loew LM. Where does all the PIP2 come from? J Physiol. 2007;582(3):945–951.
- 59. Li X, Edwards M, Swaney KF, Singh N, Bhattacharya S, Borleis J, et al. Mutually inhibitory Ras-PI(3,4)P2 feedback loops mediate cell migration. Proc Natl Acad Sci U S A. 2018;115(39):E9125–E9134. pmid:30194235
- 60. Fets L, Nichols JME, Kay RR. A PIP5 kinase essential for efficient chemotactic signaling. Curr Biol. 2014;24(4):415–21.
- 61. Kortholt A, King JS, Keizer-Gunnink I, Harwood AJ, Van Haastert PJM. Phospholipase C regulation of phosphatidylinositol 3,4,5-trisphosphate-mediated chemotaxis. Mol Biol Cell. 2007;18(12):4772–9.
- 62. Keener JP. A geometrical theory for spiral waves in excitable media. SIAM J Appl Math. 1986;46(6):1039–1056.
- 63. Postma M, Roelofs J, Goedhart J, Loovers HM, Visser AJWG, Van Haastert PJM. Sensitization of Dictyostelium chemotaxis by phosphoinositide-3-kinase-mediated self-organizing signalling patches. J Cell Sci. 2004;117(Pt 14):2925–35.
- 64. Chen L, Janetopoulos C, Huang YE, Iijima M, Borleis J, Devreotes PN. Two phases of actin polymerization display different dependencies on PI(3,4,5)P3 accumulation and have unique roles during chemotaxis. Mol Biol Cell. 2003;14(12):5028–37.
- 65. Xu X, Meier-Schellersheim M, Yan J, Jin T. Locally controlled inhibitory mechanisms are involved in eukaryotic GPCR-mediated chemosensing. J Cell Biol. 2007;178(1):141–53.
- 66. Artemenko Y, Axiotakis L Jr, Borleis J, Iglesias PA, Devreotes PN. Chemical and mechanical stimuli act on common signal transduction and cytoskeletal networks. Proc Natl Acad Sci U S A. 2016;113(47):E7500–E7509.
- 67. Cao Y, Ghabache E, Miao Y, Niman C, Hakozaki H, Reck-Peterson SL, et al. A minimal computational model for three-dimensional cell migration. J R Soc Interface. 2019;16(161):20190619. pmid:31847757
- 68. Belotti Y, McGloin D, Weijer CJ. Analysis of barotactic and chemotactic guidance cues on directional decision-making of Dictyostelium discoideum cells in confined environments. Proc Natl Acad Sci U S A. 2020;117(41):25553–25559.
- 69. Iglesias PA, Levchenko A. Modeling the cell’s guidance system. Sci STKE. 2002;2002(148):re12.
- 70. Levine H, Kessler DA, Rappel WJ. Directional sensing in eukaryotic chemotaxis: a balanced inactivation model. Proc Natl Acad Sci U S A. 2006;103(26):9761–6.
- 71. Wang CJ, Bergmann A, Lin B, Kim K, Levchenko A. Diverse sensitivity thresholds in dynamic signaling responses by social amoebae. Sci Signal. 2012;5(213):ra17.
- 72. Driscoll MK, McCann C, Kopace R, Homan T, Fourkas JT, Parent C, et al. Cell shape dynamics: from waves to migration. PLoS Comput Biol. 2012;8(3):e1002392. pmid:22438794
- 73. Khamviwath V, Hu J, Othmer HG. A continuum model of actin waves in Dictyostelium discoideum. PLoS One. 2013;8(5):e64272.
- 74. Bhowmik A, Rappel WJ, Levine H. Excitable waves and direction-sensing in Dictyostelium discoideum: Steps towards a chemotaxis model. Phys Biol. 2016;13(1):016002.
- 75. Cheng Y, Othmer H. A model for direction sensing in Dictyostelium discoideum: Ras activity and symmetry breaking driven by a Gβγ-Mediated, Gα2-Ric8—dependent signal transduction network. PLoS Comput Biol. 2016;12(5):e1004900.
- 76. Nishikawa M, Hörning M, Ueda M, Shibata T. Excitable signal transduction induces both spontaneous and directional cell asymmetries in the phosphatidylinositol lipid signaling system for eukaryotic chemotaxis. Biophys J. 2014;106(3):723–734.
- 77. Onsum M, Rao CV. A mathematical model for neutrophil gradient sensing and polarization. PLoS Comput Biol. 2007;3(3):e36.
- 78. Marée AFM, Grieneisen VA, Edelstein-Keshet L. How cells integrate complex stimuli: the effect of feedback from phosphoinositides and cell shape on cell polarization and motility. PLoS Comput Biol. 2012;8(3):e1002402.
- 79. Ghose D, Lew D. Mechanistic insights into actin-driven polarity site movement in yeast. Mol Biol Cell. 2020;31(10):1085–1102.
- 80. Higham DJ. Modeling and simulating chemical reactions. SIAM Rev. 2008;50(2):347–368.
- 81. Muñoz-Cobo JL, Berna C. Chemical kinetics roots and methods to obtain the probability distribution function evolution of reactants and products in chemical networks governed by a master equation. Entropy. 2019;21(2):181.
- 82. Swain S. Handbook of stochastic methods for physics, chemistry and the natural sciences. Optica Acta. 1984;31(9):977–978.
- 83. Grima R, Thomas P, Straube AV. How accurate are the nonlinear chemical Fokker-Planck and chemical Langevin equations? J Chem Phys. 2011;135(8):084103.
- 84. Higham DJ. Stochastic ordinary differential equations in applied and computational mathematics. IMA J Appl Math. 2011;76(3):449–474.
- 85. Bostani N, Kessler DA, Shnerb NM, Rappel WJ, Levine H. Noise effects in nonlinear biochemical signaling. Phys Rev E. 2012;85(1):011901.
- 86. Schnoerr D, Sanguinetti G, Grima R. The complex chemical Langevin equation. J Chem Phys. 2014;141(2):07B606_1.
- 87. Kim JK, Josić K, Bennett MR. The relationship between stochastic and deterministic quasi-steady state approximations. BMC Syst Biol. 2015;9:87.
- 88. Lawson MJ, Petzold L, Hellander A. Accuracy of the Michaelis-Menten approximation when analysing effects of molecular noise. J R Soc Interface. 2015;12(106). pmid:25833240
- 89. Pablo M, Ramirez SA, Elston TC. Particle-based simulations of polarity establishment reveal stochastic promotion of Turing pattern formation. PLoS Comput Biol. 2018;14(3):e1006016.
- 90. Ramirez SA, Pablo M, Burk S, Lew DJ, Elston TC. A novel stochastic simulation approach enables exploration of mechanisms for regulating polarity site movement. bioRxiv. 2020;
- 91. Worbs T, Mempel TR, Bölter J, von Andrian UH, Förster R. CCR7 ligands stimulate the intranodal motility of T lymphocytes in vivo. J Exp Med. 2007;204(3):489–95.
- 92. Charest PG, Firtel RA. “TORCing” neutrophil chemotaxis. Dev Cell. 2010;19(6):795–6.
- 93. Kabacoff C, Xiong Y, Musib R, Reichl EM, Kim J, Iglesias PA, et al. Dynacortin facilitates polarization of chemotaxing cells. BMC Biol. 2007;5:53. pmid:18039371
- 94. Houk AR, Jilkine A, Mejean CO, Boltyanskiy R, Dufresne ER, Angenent SB, et al. Membrane tension maintains cell polarity by confining signals to the leading edge during neutrophil migration. Cell. 2012;148(1-2):175–88. pmid:22265410
- 95. Zmurchok C, Collette J, Rajagopal V, Holmes WR. Membrane tension can enhance adaptation to maintain polarity of migrating cells. Biophys J. 2020;119(8):1617–1629.
- 96. Schneider IC, Haugh JM. Quantitative elucidation of a distinct spatial gradient-sensing mechanism in fibroblasts. J Cell Biol. 2005;171(5):883–892.
- 97. Matsuoka S, Ueda M. Mutual inhibition between PTEN and PIP3 generates bistability for polarity in motile cells. Nat Commun. 2018;9(1):4481.
- 98. Mori Y, Jilkine A, Edelstein-Keshet L. Wave-pinning and cell polarity from a bistable reaction-diffusion system. Biophys J. 2008;94(9):3684–3697.
- 99. Postma M, Roelofs J, Goedhart J, Gadella TWJ, Visser AJWG, Van Haastert PJM. Uniform cAMP stimulation of Dictyostelium cells induces localized patches of signal transduction and pseudopodia. Mol Biol Cell. 2003;14(12):5019–27.
- 100. Edwards M, Cai H, Abubaker-Sharif B, Long Y, Lampert TJ, Devreotes PN. Insight from the maximal activation of the signal transduction excitable network in Dictyostelium discoideum. Proc Natl Acad Sci U S A. 2018;115(16):E3722–E3730.
- 101. Lee S, Shen Z, Robinson DN, Briggs S, Firtel RA. Involvement of the cytoskeleton in controlling leading-edge function during chemotaxis. Mol Biol Cell. 2010;21(11):1810–24.
- 102. Ramirez SA, Raghavachari S, Lew DJ. Dendritic spine geometry can localize GTPase signaling in neurons. Mol Biol Cell. 2015;26(22):4171–4181.
- 103. Wu Z, Su M, Tong C, Wu M, Liu J. Membrane shape-mediated wave propagation of cortical protein dynamics. Nat Commun. 2018;9(1):1–12.
- 104. Graziano BR, Town JP, Sitarska E, Nagy TL, Fošnarič M, Penič S, et al. Cell confinement reveals a branched-actin independent circuit for neutrophil polarity. PLoS Biol. 2019;17(10):e3000457. pmid:31600188
- 105. Banavar SP, Trogdon M, Drawert B, Yi TM, Petzold LR, Campàs O. Coordinating cell polarization and morphogenesis through mechanical feedback. PLoS Comput Biol. 2021;17(1):e1007971.
- 106. Hayashi F, Saito N, Tanimoto Y, Okada K, Morigaki K, Seno K, et al. Raftophilic rhodopsin-clusters offer stochastic platforms for G protein signalling in retinal discs. Commun Biol. 2019;2:209. pmid:31240247