3D 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 stimulations 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 hypothesis.


Introduction
How cells sense and interpret external chemoattractant cues and use this information to 17 direct cell movement is one of the most fundamental processes of biology. Unicellular 18 organisms rely on this mechanism to seek nutrients and survive. In multicellular 19 organisms, it is a fundamental process during embryonic development [1,2] as well as 20 responsible for the proper operation of the mammalian immune system [3,4]. Perversely, 21 it is through the development of directed cell migration that cancer cells become 22 metastatic [5][6][7]. 23 Mathematical models have been fundamental in elucidating the mechanisms that 24 cells use to direct the cell migration [8][9][10]. There is a broad consensus that cells such 25 as the social amoeba Dictyostelium discoideum and mammalian neutrophils sense the 26 chemoattractant gradient through a local excitation, global inhibition (LEGI) 27 mechanism based on an incoherent feedforward loop motif that was originally proposed 28 to explain perfect adaptation [11][12][13]. By incorporating different diffusion properties on 29 the signal components, the mechanism senses static spatial gradients without 30 movement [12,13]. 31 In response to a spatially uniform stimulus, cells display an initial transient response 32 in which Ras and downstream PI(3,4,5)P 3 and F-actin activities increase and decrease 33 several times, before eventually returning close to the pre-stimulus basal states, resulting 34 in "near-perfect" adaptation. Signaling motifs responsible for perfect adaptation fall 35 into one of two broad classes: a negative feedback (NFB) loop with a buffering node, or 36 an incoherent feedforward (IFF) loop with a proportioner node [14]. In Dictyostelium, 37 experimental evidence favors the presence of IFF-based adaptation [12,[15][16][17]. The 38 LEGI mechanism, a form of IFF, assumes the existence of a fast local excitor and a slow 39 globally diffusive inhibitor. The productions of the excitor and the inhibitor are 40 independently driven by the receptor occupancy [13]. A local response regulator, which 41 is activated by the excitor and inhibited by the inhibitor, drives downstream signals.
feature of the stochastic nature of the biochemical reactions and depends on the state of 66 the dynamical system [31,32]. To our knowledge, however, all existing computational 67 models use partial differential equations and generate these fluctuations by injecting 68 noise as an additional input into these differential equations. 69 To overcome the aforementioned limitations, here we present a new model of the 70 biological signaling mechanism that regulates motility. Specifically, to account for the 71 noise accurately we eschew the partial-differential equation approach and instead 72 incorporate the reaction-diffusion processes into the URDME (Unstructured Reaction 73 Diffusion Master Equation) software [33]. This methodology does not make a priori 74 assumptions on the size of the stochastic perturbations; instead, the fluctuations are 75 inherent in the underlying chemical master equation. Hence, the noise is controlled by 76 number of the molecules of the reacting species. Moreover, we consider a realistic 77 geometry consisting of a three-dimensional cell with membrane and cortex elements. 78 Finally, we use detailed biochemical models of the receptor dynamics, LEGI and 79 excitable modules and verify the sufficiency of these three interconnected networks in 80 explaining the spatiotemporal dynamics of chemotactic signaling. 81

82
Overall system architecture 83 For the present study we have divided the signaling pathways that drive the observed 84 excitable Ras dynamics in Dictyostelium cells into three signaling subsystems (Fig. 1A): 85 1) The G-protein coupled receptor (GPCR) subsystem is responsible for sensing 86 chemoattractants.
2) The receptor occupancy information from GPCR is passed down 87 to the Local Excitation and Global Inhibition (LEGI) mechanism which is responsible all-or-nothing responses, refractory periods and wave propagation, and provides a 92 characteristic response both in the presence or absence of chemoattractant signals. We 93 have excluded downstream networks that directly influence the actin dynamics. The 94 output of our system could readily be coupled to the cytoskeletal signaling network but 95 that would likely necessitate numerous additional entities, including elements with a 96 large number of molecules (e.g., ∼10 8 molecules of actin [34]) leading to a presently 97 unmanageable computational burden for these stochastic simulations. 98 Geometry 99 We used a stationary hemispherical structure of radius 5 µm to capture the general 100 shape of adherent amoeboid cells devoid of any cytoskeletal components (as if treated 101 with Latrunculin) (Fig. 1B). Because most of the species known to be involved in 102 generating the spatial patterns are either membrane-bound, or in the cytoskeletal 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). but consider this a sink from some molecules can move in and out of. We assumed that 105 the cortex is 200 nm thick and discretized a shell of this thickness into nodes resulting 106 from an unstructured tetrahedral mesh. The allowable minimum and maximum 107 distances between the nodes were set at 100 and 300 nm, respectively. Overall, this gave 108 rise to 5,500 nodes and 16,469 tetrahedral elements with volume distribution 109 26.8 ± 4.6 × 10 −4 µm (mean±std. dev.). Because a coarser mesh affects the smoothness 110 of diffusion, a finer mesh would be preferred, but this increases the simulation cost 111 considerably and would become a major constraining factor as the number of 112 biochemical elements in the model grew (S1 FigA-D). For this reason, we compromised 113 by choosing a medium size mesh so that, in an element of average size, a single molecule 114 corresponds to a concentration of approximately 60 nM. The nodal volumes obtained 115 from the dual of the tetrahedral mesh has a leptokurtic distribution with parameters: 116 80.1 ± 16.3 × 10 −4 µm (mean±std. dev.). The nodes on the surface are denoted as the 117 membrane nodes whereas the others are assigned as cortex nodes (Fig. 1C). The volume 118 distributions of cortex and membrane nodes were similar, with means of 8.4 × 10 −3 µm 119 and 7.7 × 10 −3 µm, respectively. The mean subtracted distributions failed to reject the 120 null hypothesis using a Mann-Whitney-Wilcoxon test, indicating that the two size 121 distributions are not significantly different. The distributions of membranes nodes at 122 the basal and apical surfaces were also similar, with mean volumes of 7.3 × 10 −3 µm and 123 7.9 × 10 −3 µm, respectively.

125
GPCR signaling module 126 The initiation of chemotaxis in Dictyostelium involves binding of chemoattractant 127 molecules (ligand) to surface-bound G-protein coupled receptor (GPCR) molecules. In 128 the present study, we focused on cAR1 as the major receptor for cAMP. The 129 heterogeneous 3',5'-cyclic adenosine monophosphate (cAMP) binding in Dictyostelium 130 has been modeled using three states indicating different levels of affinities [35]. For 131 simplicity, we ignored the sequential binding dynamics and followed the classification of 132 GPCRs on the basis of affinity only. The three states of unoccupied GPCRs have high 133 (H), low (L) and slow (S) binding affinity with respect to extracellular cAMP [35]. The 134 corresponding occupied receptor states are labeled H:C, L:C and S:C, respectively.

199
At the same time, near perfect adaptation is observed in the activated Ras response [16]. 200 This suggests that adaptation happens upstream of Ras and, as the response of G βγ is 201 persistent during cAMP stimulation, we treated this as the excitation process in the 202 LEGI module (Fig. 3A). The parameters for the G βγ dynamics were chosen ( Table 2) to 203 match experimentally measured half-times of dissociation (on application of a saturating 204 stimulus) and reassociation (on removal of stimulus) of the G α2 and G βγ subunits [17]. 205 Table 2. Parameters for G-Protein Dynamics. The parameter values are estimated to match the half times of dissociation (t 1/2,diss. ) and reassociation (t 1/2,reass. ) from [17] ntot,G-Protein ntot,GPCR  We simulated the system and looked at the dissociated subunits and observed a 206 highly nonlinear relationship between the response curves of receptor occupancy and the 207 dissociated G-protein (S3 FigA). The fraction of dissociated G-proteins decreased as the 208 total number of G-proteins increased (S3 FigA); however, as a function of the maximal 209 response, it was independent of the total number of G-proteins reaching 50% 210 dissociation at ∼23% R.O. (S3 FigB) , indicating the existence of "spare receptors" in 211 the system [48]. Additionally, we observed that fluctuations (mean coefficient of as for the higher values there was not difference in the relative fluctuation level.

215
Unlike the excitation process, inhibition is G-protein independent [17] in 216 Dictyostelium. To account for this, we assumed an inhibitor existing in both active (I ) 217 and inactive (I) forms. The diffusion constants for the inhibitor states were assumed 218 high (a number typical of cytosolic entities [47]) to satisfy LEGI requirements that the 219 global inhibitor have higher diffusivity than the local G α2βγ and G βγ . To account for a 220 possible multistep translocation to and from the membrane, we made the inhibitor 221 dynamics slow compared to those of G βγ .

222
The interaction of the excitor G βγ and the activated inhibitor jointly regulate a 223 response regulator, RR, whose action can be modeled using either a difference or a ratio 224 scheme, depending on how the activation and inhibition processes regulate it (Fig. 3B between G βγ and I (Suppl. S3 FigC, S1 File). An alternative realization, the ratio 229 scheme, involves G βγ -dependent activation and I -dependent inactivation of RR along 230 with independent basal activation/inactivation and leads to a steady-state of RR that is 231 proportional to the ratio of affine functions of the G βγ to I ) (S3 FigD, S1 File). As 232 there are a number of potential biochemical elements that could serve as the response 233 regulator, we considered both schemes. We implemented these schemes explicitly, using 234 the corresponding reactions, as well as implicitly, using a quasi-steady-state 235 approximation of the reactions. By ignoring the transient dynamics, the implicit 236 formulations served to reduce the computational burden of explicitly simulating these 237 systems (S1 File).

238
One of the possible limitations of the adaptive networks described above is the 239 sensitivity of the adaptation to noise, resulting in variances that do not adapt, but 240 depend on the stimulus level [49]. As a third alternative, we adopted the Antithetic

241
Integral Feedback (AIF) model [50] by making some of the terms local and global, and 242 then compared its performance with the two schemes described above (Fig. 3B) Fig S4A,B). In all variants considered, the initial responses depended on the stimulus 254 level; i.e., the peak amplitude (red) increased and the peak time (orange) decreased 255 with %R.O. (Fig. 3F). In all three networks, the mean nodal response regulator activity 256 returned to prestimulus levels (blue in Fig. 3F). However, the steady-state variance 257 showed a stimulus-level dependency that increased for the difference scheme, decreased 258 slightly for the ratio scheme, and was smallest and fairly constant for the AIF scheme peak amplitude of the initial response for the ratio scheme was higher than for the 265 difference scheme. The lower coefficient of variation in the peak response profile of the 266 difference scheme provides more certainty of response to change of stimulus than the 267 ratio scheme. Upon removal of the stimulus, the difference scheme returned to the basal 268 level more quickly than to the ratio scheme (S4 FigC,D). For the rest of study, we used 269 the implicit difference scheme of LEGI to avoid unnecessary repetition. 270 We next simulated the response to a gradient generated by releasing 271 chemoattractant from a micropipette. To this end, we imposed a stationary 3D

272
Gaussian profile centered around an edge point at the basal plane and considered both 273 the apical and basal (Fig. 4A) and perimeter responses ( Fig. 4B; S1 File). Following 274 imposition of the gradient, free G βγ rose quickly at the front where %R.O. was highest; 275 at the rear the rise was slower (Fig. 4C). In contrast, the inhibitor rose more slowly and 276 was fairly uniform around the perimeter of the cell, owing to the relatively high 277 diffusion (Fig. 4C,D). The resultant response regulator rose sharply at the front and 278 dropped below the basal level at the rear (Fig. 4C,D). Note that the steady-state level 279 of the response regulator at the front was lower than the peak, as the level of the 280 inhibitor increased, but still displayed levels above basal. To examine the role of 281 diffusion in these patterns, we varied the inhibitor diffusion over a wide range (S5 Fig). 282 Increasing the ratio of inhibitor-to-G βγ diffusion resulted in greater differences in the 283 response regulator between front and back. The mean difference between RR molecules 284 between front and back doubled for a 10-fold difference in diffusion coefficients.

285
At steady state, the linear profiles of local entities such as R.O. and G βγ displayed 286 gradients that were more diffused on the top surface than on the basal surface 287 (Fig. 4E-G). In contrast, the inhibitor profile was relatively flat along the length of the 288 cell, though random variations were apparent. The resultant response regulator profile 289 showed a gradient and was higher/lower at the front/rear relative to the basal level. We 290 repeated the simulations for a shallower gradient and observed similar behavior, with a 291 smaller degree of localization and slower response in RR (S6 Fig).   . 12. cytoskeleton-based network, and a slower signaling network that drives the cytoskeletal 295 system [24,25,44]. As we are modeling cells lacking an intact cytoskeleton, we focused 296 on the latter. Our proposed STEN consists of five species: RasGDP, RasGTP, activated 297 and inactivated protein kinase B substrates (PKB*s and PKBs, respectively), and 298 membrane phosphatidylinositol bisphosphate (PIP 2 ), which represents contributions 299 from both PI(3,4)P 2 and PI(4,5)P 2 (Fig. 5A). In our scheme, RasGTP acts as the activator of the excitable system and serves as a 301 front marker of the cell. Recent experiments have demonstrated that lowering PI(4,5)P 2 302 results in increased Ras activity [43]. Similarly, lowering PI(3,4)P 2 increased Ras 303 activity through the regulation of RasGAP2 and RapGAP3 [51]. Thus, we incorporated 304 the PIP 2 -mediated hydrolysis of RasGTP to RasGDP into the model. This closes a 305 positive feedback loop that is formed through mutual inhibition between RasGTP and 306 PIP 2 .

307
A slower, negative feedback loop is achieved through the RasGTP-mediated 308 activation of PKBs. This activation is achieved partly by having PH-domain containing 309 PKBA translocate to the membrane to bind to PI(3,4,5)P 3 , as well as TorC2-mediated 310 phosphorylation of PKBR1 [52]. This negative feedback loop is closed as activated  Finally, we coupled this excitable network to the LEGI in two ways. First, through 321 an RR-dependent term that converts RasGDP to RasGTP, consistent with the 322 possibility that the RR is a RasGEF or activates a RasGEF. Second, we also included a 323 term to account for the activation of Phospholipase C (PLC) by G α2 which leads to 324 PI(4,5)P 2 hydrolysis [55]. Table 4 lists the detailed reaction terms and parameter values. 325  Combining the GPCR, LEGI, and STEN modules 326 We first simulated our model in the absence of any stimulus. We observed characteristic 327 excitable behavior (Fig. 5B,C). When viewing a single node, RasGTP and PIP 2 showed 328 mutually exclusive behavior, in which the number of molecules at any one time of one 329 species dominated the other. For example, when PIP 2 dominated, there was 330 approximated 40 molecules of PIP 2 compared to approximately one molecule of 331 RasGTP (Fig. 5B, t = 30 s). The transition to a state in which RasGTP dominated 332 (∼25 molecules of RasGTP vs. ∼0-2 molecules of PIP 2 ) was rapid. In contrast, the 333 number of PKB*s molecules varied considerably less, ranging from ∼40 to ∼60 with 334 slow increases happening after the node transition to a high RasGTP state (Fig. 5B, 335 t = 60 s). When viewed across the surface of the cells, we observed wave activity 336 (Fig. 5C,S7 FigA,B,S1 Video). The cell was typically in a back state in which high PIP 2 337 levels dominated. However, when stochastic perturbations led to a spot of high RasGTP 338 activity, a wave of activation swept across the cell, moving between the basal and apical 339 surfaces. This eventually extinguished and the cell returned to its basal back state. following short or long pulses of chemoattractant stimulation, consistent with the notion 360 that the underlying signaling network is excitable [7,24]. To test this in our model, we 361 repeatedly applied short (2 s) and long (30 s) global stimuli to models with varying 362 parameter values and compared the respective RR and RasGTP responses (Fig. 5G).

363
The response profiles were nearly indistinguishable, with both peaking about 7-8 s after 364 the application of the stimulus, followed by a return to the pre-stimulus level. being longer than the duration of a pulse from the excitable system; it has been seen 374 experimentally both with signaling [56] and cytoskeletal biosensors [57].

375
Cells displays refractory periods following stimulation during which further 376 excitation fails to trigger a response, or diminished in intensity [7,24]. We applied a 377 series of double pulses, each of duration 2 s with variable delays (10-90 s) and compared 378 the peak of the second response to the first response (Fig. 5H, left panel). For short 379 delays (< 10 s) there was no distinct second response. However, as we increased the 380 time delay, we observed an increase in the peak amplitude of the second response. After 381 a sufficiently long delay (80 s), the second peak almost matched the first peak with a 382 50% recovery occurring at ∼50 s. (Fig. 5H, right panel).

383
When stimulated by a gradient, cells display persistent high levels of activity at the 384 side of the cell facing the chemoattractant source. We simulated these experiments by 385 introducing a gradient of receptor occupancy across the cell. RasGTP showed a 386 localized persistent patch of high activity facing the side with highest receptor 387 occupancy (Fig. 6A, S5 Video). We also simulated two experiments in which the spatial 388 gradient was combined with a temporal stimulus. In the first, we introduced a gradient 389 and waited until the cell displayed a spatial response; we then removed it for a variable 390 time, before finally reintroducing it [58]. The crescent disappeared following removal of 391 the stimulus but did not return to its full strength until the delay was ∼60 s (Fig. 6B,

392
S6 Video), matching our previous observations in the double pulse experiment. In the 393 second simulation, we applied a large global stimulus following the establishment of a 394 crescent in response to a gradient, and then removed all cAMP. In this case, the global 395 stimulus elicited a response everywhere (Fig. 6C, S7 Video). These simulations show the 396 complex interactions of spatial and temporal components of coupled LEGI and STEN 397 systems.  [7,44]. In this case we saw elevated levels of activity and the cell commenced 405 periodic whole-cell increases in activity similar to what have been observed 406 experimentally (Fig. 7A). Similar, though less acute increases in activity were seen in 407 simulations in which PKB*s levels in the cell were lowered through a reduction in the 408 RasGTP-mediated activation rate (Fig. 7B). In this case, bursts of wave initiations were 409 observed, but the resulting waves did not cover the whole cell surface.

410
Experiments have also demonstrated that mechanical contacts can trigger excitable 411 behavior [59][60][61]. Finally, we considered the possibility that the threshold was 412 differentially regulated by mechanical contact with the substrate, with a lower threshold 413 at the basal surface than at the apical surface. To take this into account, we assumed 414 that the PKBs-mediated inhibition of RasGTP was different between the basal and 415 apical surfaces, with lower inhibition at the former (Fig. 7C). In these simulations we 416 Fig 7. Effect of threshold on STEN dynamics. (A,B) The kymographs show the effect of lowering the STEN threshold by inhibiting PIP 2 (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. observed frequent wave initiations at the basal surface, but these waves were quickly 417 extinguished when they reached the apical surface (Fig. 7D). At the apical surface we 418 rarely saw any de novo waves, and when they did appear, they did not spread like the 419 basal counterpart and were extinguished rapidly (Fig. 7E). observations, such as receptor-ligand interactions [35], adaptation [13,16,63], 426 amplification [64], wave propagation [65][66][67], excitability [25,47,68], or polarity [69][70][71], 427 by concentrating on a specific functional block sometimes ignoring its connection to the 428 overall network. Part of this problem is a lack of experimental data as most studies consider molecules and species that exist in relatively small numbers [31,32].

443
Our simulations here demonstrate, however, that even in signaling system that consist 444 of biochemical species with high numbers, noise can play an important role in regulating 445 cell physiology. The stochastic framework that we used to simulate our model does not 446 require an artificial injection of noise. Instead, noise is a natural consequence of the 447 stochastic description of the reaction-diffusion equations. In this stochastic setting, the 448 LEGI mechanism adapts perfectly as in the deterministic setting, as long as we focus on 449 the mean level of activity (Fig 3D, S4A,B). However, we showed that the variance in the 450 response increases as a function of the stimulus strength (Fig. 3F) and that the size of 451 this increase depends on the particular way that the LEGI is implemented. The increase 452 in the variance observed could account for chemokinesis, the increased speed of random 453 migration seen following a uniform chemokine stimulus in neutrophils [73]. Moreover,  [74]. Simulations of cells stimulated by gradients of chemoattractant show 462 that the RR activity remains high at the cell front and is suppressed elsewhere. One 463 would expect that the resultant spatially localized lowered threshold would trigger 464 activity that would travel away from their point of origin. Further activity would be 465 possible, after the delay imposed by the refractory period. To explain this, we modified 466 the STEN dynamics such that, following stimulation, the system undergoes an 467 excitable-to-bistable transition (S7 FigG,H). In this case, a persistent, elevated level of 468 RR leads to a new "high" state of activity in the STEN, and persistent crescents were 469 observed. Though our implementation is similar to the wave pinning scheme used to 470 explain persistent polarization [75], it differs due to the fact that our system is only in 471 the bistable region in the presence of a persistent high stimulus.  Experimentally, several pharmacological/genetic perturbations have been used to study 482 this effect of alteration in Dictyostelium cells [43,44,76]. Recent publications suggested 483 of possible existence of a number of mechanical feedbacks on the leading edge 484 protrusion [59,77] and hence on the system threshold. Cao et al. reported that the 485 basal surface waves are mostly restricted at the bottom surface unless the threshold of 486 the system is lowered [60]. In that case, on reaching the bottom boundary the basal 487 surface waves starts travelling upwards long the apical surface.

488
Finally, while our study has centered on chemoattractant signaling, we note that 489 various components, such as GPCR signaling, the presence of feedforward adaptive 490 mechanisms, and biochemical excitability are concepts that extend well beyond the 491 realm of cell migration. For example, a recent report suggests that stochastic effects 492 may play an important role in rhodopsin signaling [78]. Thus, our study provides a 493 means for analyzing the effect of noise in these systems and could help to understand 494 more complex interactions in these other settings.   Table 3).   Table 3).