Three-dimensional stochastic simulation of chemoattractant-mediated excitability in cells

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.

Fig 1. Model schematic and simulation domain.

(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).

Fig 2. GPCR signaling.

(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.

Fig 3. Response of the LEGI mechanism to global stimulation.

(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).

Fig 4. Steep 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 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).

Fig 5. Response of STEN to global stimulus.

(A) Schematic of STEN showing the entities: RasGTP, RasGDP, PIP₂ 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), PIP₂ (red) and PKB*s (blue) at a random node that has fired spontaneously. (C–E) Spatiotemporal profiles of RasGTP (green) and PIP₂ (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 PIP₂ responses (left), and same at a single random node (center) and the spatiotemporal profile (RasGTP, PIP₂) 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. PIP₂ peaks at approximately 1,000 molecules/μm² similar to that reported in [58]. (G) Basal subtracted normalized RR ($\overline{\text{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.

Fig 6. Response to combinations of temporal and gradient stimuli.

(A–C) Shown are the temporal profiles of RasGTP (green), PIP₂ (red), PKBs*(blue), G_α2 (teal), G_βγ(orange) and basal subtracted normalized RR ($\overline{\text{RR}}$, cyan) at single nodes at the front and back of the cell. The kymographs (right) show RasGTP (green) and PIP₂ (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).

Fig 7. Effect of threshold on STEN dynamics.

(A,B) The kymographs show the effect of lowering the STEN threshold by inhibiting PIP₂(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.