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

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.


Specifying reactions in URDME
Here we demonstrate the process of specifying reactions in particular to our system. Additional details can be found in the Github page for URDME (https://github.com/URDME/urdme/blob/ master/doc/manual.pdf). We used a 3D spherical shell to mimic the cell shape where we defined the nodes on the outer surface as part of the "membrane" subdomain whereas the interior as the "cortex" subdomain. To restrict the diffusion of the molecules of exclusively membrane bound species, for example, receptors to the cortex, the entries in the diffusion matrix correspond to the cortex elements were manually made zero.
As an example of how to translate the reactions of Tables 1-4 into the URDME code, we consider reaction 1 from Table 1 In this case, PKB * s (PKBsa) does not appear in the left side as in the previous example because it acts as an enzyme in the conversion reaction from RasGTP to RasGDP. When considering the effect of a spatially-differentiated thresholds, we used different values for the reaction constant a 2 (a 20 , a 21 ) in en2 for different subdomains. In this context, we divided our simulation domain further into subdomains "basal" (sd = 0) and "apical" (sd = 1). In this case, the reactions were entered as follows.
In the basal domain: Finally, an example of a generation reaction is no. 5 in Table 4. The reaction is with propensity function b 1 . In this case, the reaction is written as: The reaction is similar to a birth process where the symbol "@" implies the production is de novo and b 1 is in µMs −1 ; the extra "vol" term is multiplied.

Implicit and explicit LEGI schemes
Though our study relies exclusively on stochastic descriptions of the reactions (as illustrated in the previous section), to explain the difference(s) between the "implicit" and "explicit" formalisms of the LEGI scheme we will use deterministic description based on PDEs.
The explicit dynamics (Fig. S3E,F) in the deterministic setup can be represented as follows.
First, for the difference scheme: Assuming a quasi-steady-state and a spatially-homogeneous solution, we obtain: Replacing from 2a into 2b yields: and solving for [RR] ss we obtain: Thus, we directly used the steady-state expression as the coupling term between LEGI and STEN as follows: where BG = G βγ and Ia = I * .

Kymograph generation
To create the cell perimeter kymographs, the nodes lying between the z sections between heights 1 and 2 µm (Fig. 4B), were projected on a plane, resulting in an annular section. This was then divided into 90 sectors each with a 4 • span. The maximum value for the specific species over all the voxels in each sector was used as the representative one. Thus, at each time point we generated an array of 90 representative values. We continued this process over the intended period of time to generate the respective kymographs (Fig. 4D, S6 Fig.D, Fig. 6A-C, Fig. 7A,B, S7 Fig.F).

Parameter selection and parameter variation
Whenever the direct experimentally measured rate constants [1,2] were available, we used the values in our model (Table 1). In other cases, we relied on some indirect approaches. For example, the experimental half times of G α2βγ association and dissociation were available [3], so we chose reaction rates, k E0 , k E , k −E ( Table 2) to achieve the similar time constants from simulated response. The total number of G-protein (G α2βγ +G βγ ) and the LEGI-inhibitor (I+I * ) were assumed to have the 210,000 molecules on average. The rest of the parameters of LEGI schemes were chosen to achieve similar mean level of steady-state values and adaptation time. The adaptation time was computed as the time taken by the mean nodal profile of RR to reach within 2% of the final steady-state value. The peak amplitude was computed as the maximum value of the smoothened mean nodal profile of RR and the peak time was the time taken to reach that after the stimulus was applied. Unlike the LEGI implicit difference and ratio scheme, in LEGI-AIF, the species X, Y and RR were explicitly modeled (reactions no. 7-15 in Table 3). When choosing the reaction rate constants for the LEGI-AIF scheme, we tried to suppress the oscillations in the mean profile of RR. For saturating dose of cAMP, the maximum value of the total number of intermediate X and Y were approximately 30,000 and 45,000, respectively.
For the selection of the parameter set for the excitable system, we relied on the experimentallydetermined wave speed, qualitative and quantitative behaviors of different entities [3][4][5][6][7]. Most of the reactions in STEN module are single-step reactions except for the autocatalytic reaction (reaction 3 in Table 4). The complex propensity function of reaction 3, involves several steps: inhibition by PIP2, basal activation as well as RR-dependent activation. To reduce the computation burden of simulating these first/ second-order reactions as well as estimation of more number of unknown reaction parameters, we simplified this multistep reaction process using quasi-steady-state assumption. PKBA, a component of PKB*s shuttle between membrane and cytosol, we assumed that inactivated PKBs is in the cytosol which on activation comes to the membrane as PKB * s. We did not model the activation/inactivation of PKBs explicitly, rather we assumed it as a birth/death process (reactions no. 9&10 in Table 4).
In the present study, we investigated the effect of the intrinsic noise of the system as well as the effect of the cell-to-cell variations. In the latter, we varied all the parameters (both reaction rates and the diffusion constants) following a Gaussian distribution with nominal values and standard deviations. All the nominal values and the standard deviations (typically 10%) are listed in the respective tables. We assumed homogeneity in the cell but heterogeneity across the cells in terms of the parameters.

Computational and Statistical analysis
All simulations were run using URDME package and custom codes written in MATLAB 2019a (MathWorks, Natick, MA, USA) on an macOS (version 10 and 11). The statistical analysis were performed with custom codes written in MATLAB 2019a (MathWorks, Natick, MA, USA) and RStudio 3.5.1 (RStudio, Boston, MA, USA). The statistical tests used, as indicated in this manuscript, were as follows: Mann-Whitney-Wilcoxon (MWW) test, Kolmogorov-Smirnov (K-S) test. For all tests, the significance level was set to less than 0.05.