Spatial models of pattern formation during phagocytosis

Phagocytosis, the biological process in which cells ingest large particles such as bacteria, is a key component of the innate immune response. Fcγ receptor (FcγR)-mediated phagocytosis is initiated when these receptors are activated after binding immunoglobulin G (IgG). Receptor activation initiates a signaling cascade that leads to the formation of the phagocytic cup and culminates with ingestion of the foreign particle. In the experimental system termed “frustrated phagocytosis”, cells attempt to internalize micropatterned disks of IgG. Cells that engage in frustrated phagocytosis form “rosettes” of actin-enriched structures called podosomes around the IgG disk. The mechanism that generates the rosette pattern is unknown. We present data that supports the involvement of Cdc42, a member of the Rho family of GTPases, in pattern formation. Cdc42 acts downstream of receptor activation, upstream of actin polymerization, and is known to play a role in polarity establishment. Reaction-diffusion models for GTPase spatiotemporal dynamics exist. We demonstrate how the addition of negative feedback and minor changes to these models can generate the experimentally observed rosette pattern of podosomes. We show that this pattern formation can occur through two general mechanisms. In the first mechanism, an intermediate species forms a ring of high activity around the IgG disk, which then promotes rosette organization. The second mechanism does not require initial ring formation but relies on spatial gradients of intermediate chemical species that are selectively activated over the IgG patch. Finally, we analyze the models to suggest experiments to test their validity.


Introduction
All cells must be able to respond to changes in their environment, and often the proper response requires cells to adopt a new morphology. For example, cell shape changes occur during migration, division, and phagocytosis. Typically, these changes are initiated when receptors on the cell surface are activated by an external cue [1]. Receptor activation initiates a signaling cascade that results in spatiotemporal regulation of the actin cytoskeleton. The Rho family of GTPases are a class of signaling molecules that play key roles in this process [2][3][4][5][6]. These proteins act as molecular switches. They are in an inactive state when bound with GDP and become active when GDP is exchanged for GTP. Once active, Rho GTPases interact with effector molecules including those that regulate the actin cytoskeleton. Due to the nonlinear nature of the signaling pathways that regulate GTPase activity, understanding the molecular mechanisms that generate cell shape changes has proven challenging [1]. Therefore, many recent studies have turned to mathematical modeling to explore mechanisms capable of generating complex molecular structures [7][8][9][10][11].
Here we focus on Fcγ Receptor (FcγR)-mediated phagocytosis because of its biological importance in the innate immune response [12,13] and because phagocytosis provides an ideal system for studying how Rho GTPases organize the cytoskeleton into well-defined structures. Phagocytosis is initiated by the binding of the antibody immunoglobulin G (IgG) to FcγR. Upon FcγR clustering, receptor cross-linking leads to phosphorylation of activation motif domains, enabling downstream signaling [12][13][14]. To study the events that initiate phagocytosis under well-controlled conditions, IgG is micropatterned in small disks on a glass coverslip ( Fig 1A). Because the antibody is attached to the coverslip it cannot be internalized, and the experimental system is therefore referred to as "frustrated" phagocytosis [15]. Following receptor activation, actin-enriched, adhesion-like structures termed podosomes [13,16,17] form in a circle around the IgG disk (Fig 1B and 1C). Podosomes classically have been considered rod or cone-like structures of dense actin [16,17]. However, we recently demonstrated an hourglass-like shape [18]. Podosomes recruit many additional molecules and are thought to coordinate interactions between the actin cytoskeleton and the extracellular matrix [16,17,19]. They also form the leading edge of the phagocytic cup [20,21], and have been referred to as "teeth" coordinating the "jaw" during phagocytosis [21]. The mechanisms responsible for podosome formation and patterning are not known. Therefore, we turned to mathematical modeling to establish sufficient conditions for pattern formation during frustrated phagocytosis.
Beginning with Turing's seminal paper [22] and continuing with developments by Gierer and Meinhardt [23] and Meinhardt [24], reaction-diffusion models have been used to investigate pattern formation in biological systems. These models rely on positive feedback to amplify local fluctuations in signaling activity and some form of global inhibition or substrate depletion to keep regions of high activity localized [23]. Another key requirement of these models is that at least one of the chemical species in the system diffuses at a different rate from the others [22,23]. The hydrolysis cycle of GTPases satisfies the requirements for spontaneous polarization [7][8][9][10]25,26]. GTPases cycle between an active state when GTP-bound and an inactive state when GDP-bound. Their activation is catalyzed by guanine nucleotide exchange factors (GEFs), which promote the exchange of GDP to GTP. This exchange typically occurs at the cell membrane where diffusion is slow as compared to the cytosol [7,8,10,26,27]. When in the active state, some GTPases have been shown to recruit their own GEFs forming a positive feedback loop [25,26,[28][29][30][31]. GTPase inactivation is accelerated by GTPase-activating proteins (GAPs) [3][4][5][6]. When inactive, GTPases are sequestered in the cytosol by guanine nucleotide dissociation inhibitors (GDIs) and diffuse rapidly [3][4][5][6].
There are now many reaction-diffusion models that describe how GTPases can generate cell polarity and patterning in various systems [8][9][10][11]. One of the best characterized cases is in yeast (Saccharomyces cerevisiae) budding, in which the GTPase Cdc42 generates a single, active site to determine the location of a bud site or mating projection. In yeast, autocatalysis is welldefined: active Cdc42 binds to the scaffold protein Bem1, which also binds to the GEF Cdc24 that locally activates more GTPase [25,26,28,29]. Other examples of Rho GTPase-driven pattern formation include single-cell wound healing, in which Rho and Cdc42 form distinct rings through regulation by the dual GAP-GEF Abr [32,33], actin waves observed during cytokinesis that are driven by RhoA and its GAP RGA-3/4 and GEF Ect2 [34,35], RhoA driven pulsed contractions observed during embryonic development in C. elegans [34], and tip growth in pollen tubes and fungal hyphae [10].
Here, we expanded upon existing reaction-diffusion models for GTPase activity to demonstrate how these systems can generate the "rosette" pattern of podosomes observed during frustrated phagocytosis. We explore the behavior of a recent model for GTPase activity that includes a negative feedback loop formed through the activation of a GAP [9,10]. Depending on the choice of parameter values, this model generates a range of patterns including spots, mazes, and inverse spots [10]. We use the model to identify two distinct mechanisms for generating a rosette pattern. In the first scenario, an intermediate species forms a ring of activity that promotes the formation of active GTPase spots in the ring. Next, we use a parameterization approach involving an evolutionary algorithm followed by Markov chain Monte Carlo to evolve systems that do not rely on initial ring formation to generate the rosette pattern. A common theme that emerges from this analysis is that rosette formation requires the activation of both a positive and negative regulator of GTPase activity over the IgG disk. This creates spatial gradients of these regulators, which in turn are sufficient to drive the formation of the rosette pattern. Finally, we analyzed the behavior of the models to suggest experiments to test our proposed mechanisms.

Experimental observations suggest Cdc42, but not myosin, is required for rosette patterning
Macrophages (RAW 264.7 cells) were observed during frustrated Fcγ receptor IIa (FcγR) mediated phagocytosis, where cells attempt to phagocytose fixed, micropatterned disks of immunoglobulin G (IgG). Actin, a major downstream effector during FcγR-mediated phagocytic signaling, formed in rings of small puncta, just outside of the IgG disks (Fig 1A-1C). These puncta were podosomes: actin-rich, adhesion-like structures observed during phagocytosis but more commonly known for their roles in motility and extracellular matrix interactions [16,19]. This superstructural organization of podosomes in a circular arrangement has previously been termed a podosome "rosette" [36][37][38]. Due to the dynamic nature of phagocytosis, actomyosin contractility is known to play an integral role during the engulfment process [12,13,21,39] and myosin II has been observed to localize to phagocytic podosomes and podosome rosettes [17,21]. Therefore, we wondered whether actomyosin contractility was important for podosome rosette formation. To test this possibility, we treated cells with the Rho kinase inhibitor Y27632. Inhibition of Rho kinase during frustrated phagocytosis led to the complete disassembly of myosin II filaments, demonstrating that myosin II contractility was inhibited (S1A and S1B Fig). However, podosome rosettes still formed (S1A and S1B Fig), which suggested that the formation and maintenance of podosomes during phagocytosis is independent of actomyosin contractility and that a biochemical mechanism may underlie rosette formation.
Rho family GTPases, including Cdc42, are known to be activated during FcγR-mediated phagocytic signaling [2,12,13,40,41]. Cdc42 is a regulator of the actin cytoskeleton, so we next examined its localization during frustrated phagocytosis. We recently developed tools for observing and analyzing single molecule conformational changes in living cells [42]. We made use of these techniques to visualize Cdc42 during frustrated phagocytosis (Figs 1D-1G and S1C). Cdc42 localized to the podosome rosette, with individual molecule tracks clustered near podosomes (Fig 1D-1G).
Taken together these results suggest podosome rosette organization involves localized Cdc42 activity but does not require active myosin-mediated force generation. The involvement of Cdc42 in rosette formation is also supported by other studies; Cdc42 levels are reduced when actin is reduced at the phagocytic site [43], and Cdc42 is recruited to the tips of pseudopodia early in phagocytosis [40,44]. Therefore, we decided to investigate if a mechanism involving Cdc42 might underlie formation of the podosome rosette.

Forming coexisting clusters of active GTPase
The core components of mathematical models for polarity establishment include an inactive form of a GTPase that is cytosolic and diffuses rapidly, an active form that is membrane bound and diffuses slowly, and positive feedback through autoactivation [7][8][9][10][11]. Additionally, these models often assume that mass is conserved and, therefore, do not include protein synthesis and degradation. In their simplest form, these models typically form a single polarity site [7][8][9][10]25]. Recent investigations have focused on establishing mechanisms that generate coexisting active GTPase clusters. Chiou et al. investigated a minimal model of a GTPase polarity circuit that assumed mass conservation and did not include any form of negative feedback [8]. They demonstrated how, in this model, local depletion increases the competition time between clusters so that coexistence is maintained over biologically relevant time scales. To generate stable coexisting sites requires the inclusion of protein synthesis and degradation [10,45,46], adding negative feedback to limit the growth of active cluster [10,46], or including an intermediate state of the Rho GTPase, so it cannot be immediately reactivated following deactivation [9]. We note that any of these mechanisms might underlie rosette pattern formation. However, the GAPs ARHGAP12, ARHGAP25, and SH3BP1 have been shown to play a role in phagocytosis [47]. Therefore, we focused on the Wave-Pinning GAP (WPGAP) model studied by Jacobs et al. in which negative feedback occurs through activation of a GAP [10].
The WPGAP model (Fig 2A, [10]) is described mathematically by the following set of reaction-diffusion equations: where u is the concentration of active GTPase, v is the concentration of inactive GTPase, G is the concentration of active GAP, and g is the concentration of inactive GAP. The basal GTPase activation rate is b, the maximum self-positive feedback rate is γ, K is the concentration of active GTPase when the feedback is at the half-maximal response, the basal GTPase inactivation rate is σ, the GAP-mediated negative feedback rate is e, the GAP activation is c, and the GAP inactivation rate is d. The total mass of both species is conserved: where the integrals are over the volume of the system. A requirement for polarization is that the membrane-bound active form of GTPase diffuses slowly in comparison to the cytosolic inactive form: Both the active and inactive forms of the GAP are treated as cytosolic species that diffuse rapidly compared to the membrane-bound active GTPase. It is tempting for simplicity and computational efficiency to perform initial investigations on a 1D ring to find regions of parameter space that support coexisting spots. However, the existence of coexisting clusters on a 1D ring does not guarantee a localized ring of clusters will form when the system is extended to 2D. One important difference between 1D and 2D is that  (Table 1). Red arrows show most parameters can transition from a low intensity homogenous regime (circle cap) to a spot patterning regime (arrow cap). C) Spatial profile of parameters modulated by species M. M either activates (first 3 panels) or inhibits (last 3 panels) parameters. D) Active GTPase concentrations using the parameter profiles shown in C (above each panel, respectively). Consistent with parameter sweeps shown in B, each parameter is capable of being modulated to form a rosette pattern except for the basal activation rate b.
https://doi.org/10.1371/journal.pcbi.1010092.g002 in 2D spots have curvature, which can contribute to their stability [8]. Furthermore, our system requires low activity without patterning away from the spots, which could not be captured in 1D. Therefore, we performed WPGAP simulations in 2D.
For the WPGAP model, Jacobs et al. [10] explored how changes in the total amount of GTPase (T) impacted bistability and the types of patterns that formed. They found the system could form spots (localized regions of high GTPase activity), mazes, and negative spots (localized regions of low GTPase activity). In this study, we are interested in coexisting spots, and thus, guided by the work of Jacobs et al. [10], we began with parameter values that placed the system in this regime. From single particle tracking during frustrated phagocytosis (Fig 1D-1G), we estimated the diffusivity for the membrane-bound form of Cdc42 as~0.04 μm 2 s -1 (S1C Fig). Using this information, we fixed D u = 0.04 μm 2 s -1 and we let D v = D g = D G = 100 D u . Radii of podosomes have been observed to be anywhere from~0.15-0.6 microns [17,18,48,49]. Therefore, we tuned the rate constants (Table 1) so that the model produced spots with radius 0.31 ± 0.02 μm. Because most of the rate constants in the model have not been measured experimentally, our choices for the initial parameter values are somewhat arbitrary. However, these rate constants result in coexisting spots of appropriate size when using the estimated membrane-bound diffusion coefficient for Cdc42.
To explore how GTPase cluster size depends on the relative time scales of diffusion and chemical kinetics, we fixed the rate constants and varied the diffusion coefficient of the active GTPase D u , keeping the other diffusion coefficients 100-fold greater than D u (S2A Fig). Note that mathematically this is equivalent to fixing the diffusion coefficients and varying the time scale for the chemical kinetics, because solutions to the model equations only depend on the ratio of these two time scales. We found that increasing the diffusivity increased the spot size (S2A and S2B Fig), but did not impact the patterns formed by the system (quantified by spot eccentricity, S2C Fig).
We set the average Cdc42 concentration to 808 molecules per μm 2 ( Table 1). The number of Cdc42 molecules per podosome is not known. However, particle-based simulations of Cdc42 polarity in yeast typically predict around 3000 molecules in a polarity site of roughly 1 μm radius, resulting in a surface concentration of~950 molecules per μm 2 [50]. Using a GTPase spot radius of 0.31 μm, the area of a podosome is approximately 0.30 μm 2 , and there are approximately 240 Cdc42 molecules per podosome. This produces a surface concentration of~900 molecules per μm 2 at active spots, which is consistent with particle-based models for the polarity site in yeast. Next, we explored how varying individual parameters changed pattern formation. When sweeping the GTPase activation rate b, lower values resulted in spots, but higher values resulted in a spatially homogenous steady state with an intermediate level of GTPase activity ( Fig 2B). This suggested that a finite basal GTPase activation rate b was not required or had to be quite small to facilitate patterning. Interestingly, each of the other parameter sweeps resulted in changes in the observed patterning types, including spots, mazes, and holes ( Fig  2B). For the GTPase inactivation rate σ, a high value resulted in a single low concentration throughout the domain, and decreasing this value led to spots, then mazes. However, minimizing this value did not cause the entire domain to be at a single, high steady state, due to negative feedback from GAPs. In contrast, the self-positive feedback rate γ, the GAP inactivation rate c, the GAP activation rate d, and the GAP-mediated negative feedback rate e were capable of all patterning types, from a single low state to spots, mazes, holes, and a single high state ( Fig 2B). Overall, these observations suggested that it would be possible to spatially modulate these parameters to go from the low, homogenous state to the spot forming state.

A two-step model for rosette formation
We next sought to determine if the WPGAP model could be modified to enable rosette formation. One possible explanation for how a rosette could form is if two distinct steps occur: 1) an initial ring of high or low concentration of some species (M) forms and 2) this species modulates a key parameter in the pattern forming, WPGAP model. To test this model, we assumed a preexisting ring in the concertation profile of M. We note that rings of high GTPase activity have been observed and modeled in other contexts, such as in wound healing, in which a chemical gradient and modulation of a bistable GTPase resulted in distinct rings of activity [32,33].
For our initial investigations, we assumed that a modulator M affected a rate in the WPGAP model through the functional form: where ω 1 is the basal rate, ω 2 models the effect of M on ω ± and r measures the radial distance from the center of an IgG disk. In our simulations, we consider a single IgG disk and use polar coordinates with the origin located at the center of the disk. The computational domain consists of a disk of radius of R with reflective boundaries at r = 0 to r = R (see Methods). Unless otherwise noted R = 4 μm. M(r) was modeled as a Gaussian-shaped function centered at r = 2 μm with variable variance. This form of ω ± allowed us to tune model parameters so that spot formation was only promoted within the ring.
For parameters that increase GTPase activity (GTPase activation b, GAP inactivation d, and the maximum self-positive feedback rate γ), the WPGAP model was coupled to a ring of high M concentration ω + (Fig 2C and 2D, three leftmost columns). For parameters that decrease GTPase activity (GAP activation c, GTPase inactivation σ, and GAP-mediated GTPase inactivation e), the WPGAP model was coupled to an inverted ring of M, ω - (Fig 2C and 2D, three rightmost columns). For each model parameter, ω 1 and ω 2 were varied to determine if the system could generate rosette organization. As an initial guess, the parameter values were chosen based on the results from the parameter sweeps ( Fig 2B).
As expected from the parameter sweep results, modulating the basal GTPase activation rate b did not appear sufficient to form a rosette pattern, because this produced spot formation throughout the entire domain. However, modulating the other rate constants, such as the positive feedback rate γ, all resulted in a rosette forming (Fig 2C and 2D). Interestingly, when we modulated the rates for GTPase inactivation σ and the GAP-mediated negative feedback e, we found that the rates required to form rosettes were higher than expected (Fig 2C and 2D). For example, to form a rosette, the rate required for the GTPase inactivation σ within the ring was 15 s -1 , which resulted in no patterning when used as the global rate in the isolated WPGAP model (Fig 2B-2D).

Gradient establishment by a simple reaction-diffusion model
The analysis presented above demonstrated that the rosette pattern can form following the establishment of a ring of activity. Therefore, we next wanted to determine if rosette formation could occur in the absence of such an initial ring. Our first goal was to establish the concentration profile of a signaling molecule that is activated by IgG-bound receptors over the disk and deactivated throughout the domain. This profile was later used to model the spatially-dependent rate constants in the WPGAP model. Let X represent the active form of this molecule. To model the spatiotemporal dynamics of X, we used the following equation: where k x (r) is the spatially-dependent activation rate, δ x is the deactivation rate and D x is the diffusion coefficient for X ( Fig 3A). Note that if k x is independent of r, then at steady state X(r) = k x /δ x . To model the IgG disk, we treat k x (r) as a step-function: where k disk is the IgG-induced activation rate, k basal is the basal activation rate, and μ F is the radius of the disk. For fixed values of k disk , k basal , and δ x , this system was simulated for varying diffusion coefficients D x (Fig 3B). The range of diffusion coefficients was chosen to approximate the orders of magnitude observed in cellular diffusion rates [51,52], from a slow membrane-bound rate (10 −3 μm 2 s -1 ) to a fast cytosolic rate (10.0 μm 2 s -1 ). For a domain with a radius R = 4 μm, changing the diffusivity of X resulted in changes in both the gradient steepness and the difference between the maximum and minimum values of X ( Fig 3C). For small values of D x , the distribution of X was switch-like, and X approached the expected steady state, k x (r)/δ x . However, for larger values of D x , the gradient in X was shallower, and deviated significantly from k x (r)/δ x with a lower total amplitude (Fig 3B and 3C). We wondered if there was a functional form that approximated the solution to this equation, and found that our simulation results could be well approximated using a logistic function with the form: where β and α+β are the minimum and maximum values of f(r), respectively, k m is the logistic decay rate, and f(r 0 ) = β+α/2 ( Fig 3D). Therefore, for simplicity and computational speed we used f(r) when performing parameter searches for the WGAP model. However, after performing parameter searches, we then verified that directly simulating the diffusing species rather than using f(r) also generated rosette formation (see S1 Note).

Rosette formation through gradients
To determine if rosette formation is possible in the WPGAP model without the formation of an initial ring around the IgG disk, we treated the positive feedback rate γ and the GAP activation rate c as spatially-dependent with profiles given by f(r): Note that these functions asymptotically approach their maximum values of γ max = γ α + γ β or c max = c α + c β as r decreases. Unlike the two-step model, it was difficult to empirically determine parameter values that form a rosette. Thus, we used a two-step approach to search parameter space. We first used an evolutionary algorithm (EA) [53] to perform a global search and subsequently performed a more local sampling of parameter space using a Delayed Rejection Adaptive Metropolis Markov chain Monte Carlo (DRAM-MCMC, see Methods) [54,55].
To implement the two-step approach requires a score function that provides a quantitative measure for how close a simulated result is to the desired rosette pattern. For the desired pattern, a GTPase rosette formed by the two-step model was used (Fig 2D). From this, we measured the radial average and radial standard deviation for the active GTPase u (Fig 4A-4C).
For new simulations, we measured the radial average of u (Fig 4A and 4C). We then divided the system into octants and measured the radial standard deviation within each octant (Fig 4B  and 4C). The difference between the means of the desired output and simulation result and  the differences between the standard deviation of the desired output and standard deviations in each octant were calculated. These nine measurements (Fig 4A-4C) were then averaged to produce a single score. This score function was accurate, but flexible enough to allow for various numbers of spots and spot locations.
Simulations were initialized with a small amount of random noise and seeded with an initial concentration of active GTPase in the shape of a rosette (see Methods). We seeded simulations with a rosette to decrease the time for pattern formation, because parameter estimation requires a significant number of simulations. Also, from observations of initial parameterization attempts, the GTPase activation rate b, the GTPase inactivation rate σ, the minimum selfpositive feedback rate γ α , and the minimum GAP activation rate c α were typically quite small and were thus fixed at 2e-3, 0.4, 5e-3, and 5e-3 s -1 , respectively. For 99 individual EA runs (100 individuals, 100 generations), most runs were able to discover parameters capable of rosette organization (top 80 appeared successful, S3A and S3B Fig). The best parameter set found by the EAs was then used to initialize DRAM-MCMC simulations. DRAM-MCMCs were simulated until they appeared to converge, with all but the final 5,000 iterations removed as a "burn-in" period (S3C Fig, see Methods).
The parameter distributions generated by MCMC sampling appeared Gaussian (Fig 4D, on the diagonal). We took the mean values of the individual parameter distributions as our representative parameter set (Figs 4D-4F and S3D and Table 2). To check how well the MCMC performed, we also simulated the worst scoring parameter set, and these parameters also resulted in rosette organization (S3D and S3E Fig and S1 Table). Our simulations predicted that GTPase rosettes mostly form by 20 s and appear stable by 40 s (S4 Fig). We also verified that simulating the diffusing species (see the section Gradient Establishment by a Simple Reaction-Diffusion Model) rather than using the logistic function approximation f(r) would generate rosette formation (see S1 Note, S5 Fig).
Inspection of the spatially-dependent rates revealed how the system was capable of rosette patterning (Fig 4E). When the ratio between the positive to negative feedback (γ(r) /c(r)) is plotted as a function of r using the identified parameter sets, in all cases the ratio is maximized just beyond r = 2 μm, near where the spots formed (Fig 4E and 4F). Relative to the self-positive feedback rate γ(r) (Fig 4G, magenta), the GAP activation rate c(r) (Fig 4E, green) is high over the disk and away from the disk. However, c transitions more rapidly than γ between its elevated level over the disk and its basal level away from the disk (Fig 4E). Thus, while the negative feedback is relatively high over the disk and away from it, there is a zone near the edge of the disk where the negative feedback is low relative to the positive feedback. It is in this region that rosette formation occurs.
To gain further insight into the model's behavior we looked for pairwise correlations between model parameters. Several parameters demonstrated strong correlations (Fig 4D). There was a strong anti-correlation between c max , the maximum GAP activation rate, and e, the rate constant for GAP-mediated GTPase inactivation. This likely indicates a sensitivity of the model to the total amount of GAP activity. The other correlations were not as intuitively apparent, so to further explore parameter-dependent model behavior, we performed individual parameter sweeps using the representative parameter set (Table 2 and Fig 5A-5F). Parameters typically moved from no patterning to rosette organization to ring formation (c km , γ max , Individual parameter distributions are shown on the diagonal. Lower triangular plots show kernel density estimates for parameter pairs. Circles in the upper triangle represent the Spearman correlation coefficient between parameters. E) Radial profiles for the non-constant parameters from the parameter estimation. The positive to negative feedback ratio is also shown. For visualization purposes, the GAP activation rate and ratio are scaled to be similar orders of magnitude. The solid lines are the results for the mean parameter values from D (on the diagonal, Table 2) and the shaded regions indicate one standard deviation. F) Active GTPase concentration using the representative parameter set ( and d, Fig 5B, 5C and 5E), or vice versa (c max , γ km, and e, Fig 5A, 5D and 5F). Thus, the anticorrelations between c km and γ max as well as γ km and e likely result from a balancing of the effects produced by varying the individual parameters. However, the reason for the positive correlation between c max and γ km is not readily apparent but may result from the logistic function having a lower maximum value for shallower gradients (i.e., low decay rates, Fig 3). Also, note that this model can generate rings of activity (Fig 5A-5F), as is required for the two-step model discussed above.   Finally, we performed a two parameter sweep for γ and c. We restricted the sweeps to the region of parameter space where spots formed (Fig 4E and 4F). Individual simulations were performed using a constant value for γ and c (Fig 5G). For high values of γ as compared to c, the model was in a high, homogenous regime (Fig 5G, upper left). For high values of c as compared to γ, the model was in a low, homogenous regime (Fig 5G, lower right). For intermediate values of c and γ, various types of patterning occurred, from spots to mazes to holes (Fig 5G,  bottom left to upper right). We next plotted c(r) vs. γ(r) within this region using the representative parameter set given in Table 2 (Fig 5G, red curve). The curve is typically in the low, homogenous regime but passes through the patterning area, demonstrating why spots can form only within a certain spatial zone.

Experimentally testable predictions
Finally, we analyzed the model with the goal of motivating experimental investigations. First note that if the GAP responsible for deactivating Cdc42 is knocked down, the model predicts that the rosette pattern should change to a ring of activity, whereas overexpression of the GAP should result in a loss of pattern formation (Fig 5A and 5F). Next, we simulated the model using varying disk sizes (Fig 6A and 6B). We observed a linear relationship between the disk radius and the number of spots, indicating that the distance between spots remains constant as the disk size increases (Fig 6A). To compare our simulation results to experimental results, we counted the number of podosomes per site for IgG disks of radius 1. The observed number of spots qualitatively compared well with the numbers predicted by the model (Fig 6B). For simulations using disk sizes of radius greater than 0.24 μm, 3 or more spots of active GTPase formed, while smaller disks formed 2 or fewer spots (Figs 6A, 6B and S6C).

PLOS COMPUTATIONAL BIOLOGY
Similarly, we explored simulations with holes lacking IgG of various sizes (Table 2, negative c km and γ km , Figs 6C and S6D). Rosettes always formed for a hole of radius 2.1 μm, and never formed for a hole of radius 1.6 μm. However, holes with radii between 1.7 and 2.0 μm appeared to be capable of forming rosettes, but rosettes sometimes failed to organize properly depending on the initial conditions, suggesting that the system is bistable within this regime (Fig 6C).

Discussion
We used mathematical modeling to demonstrate how a Wave-Pinning GAP system with spatially-dependent parameters can generate the rosette pattern of podosomes observed during phagocytosis. Other studies have investigated the effects of spatially-varying parameters on Turing systems. Most of these investigations were motivated by pattern formation during embryonic development or morphogenesis. They revealed that spatial gradients can induce transitions between spots and stripes [56,57], demonstrated that spatial gradients can generate "Turing-like" patterns in regions of parameter space where the spatially homogenous system does not admit a Turing instability [57], and revealed that more complex patterns than spots and stripes, such as traveling waves, are possible when spatially-dependent parameters are considered [58]. Spatial gradients have also been studied at the single cell level. For example, Payne and Grierson studied how the position of the GTPase ROP polarity site is determine by a gradient of the plant hormone auxin [59]. Their model was able to explain many of the observed phenotypes of various mutants. In a follow up study, a detailed bifurcation analysis of the model was performed [60]. This analysis revealed that transitions from a boundary patch to an internal patch to the coexistence of multiple patches are controlled by the auxin gradient.
Our study was motivated by changes in cell morphology that occur in many different physiological contexts, including cell migration, division, and differentiation. In eukaryotes, cell shape changes are driven by forces generated by the actin cytoskeleton. The Rho family of GTPases are primary regulators of the actin cytoskeleton. Therefore, understanding how these signaling molecules generate spatiotemporal patterns is fundamental to understanding cellular morphodynamics. An emerging theme in Rho GTPase signaling is that pattern formation occurs through a combination of positive autoregulation and differences in diffusivity between active and inactive GTPase states. It is now well understood how these elements can generate cell polarity (i.e., determining a cell front and back). Recent studies have investigated how GTPase-driven systems can generate more complicated spatial patterns [8][9][10][11]26]. Therefore, we turned to mathematical modeling to determine whether these systems can generate the rosette of podosomes observed during frustrated phagocytosis. Experimental evidence for such a biochemical mechanism comes from our observations that Cdc42 colocalizes to rosette structures and actomyosin contractility does not appear to be needed for rosette formation.
The starting point of our analysis was the WPGAP model that was previously shown to be capable of forming co-existing clusters of high GTPase activity [10]. We used the model to investigate two potential mechanisms for rosette formation. In the first scenario, a ring of high or low concentration of a regulator of GTPase activity initially forms. Our analysis revealed that rosette formation was possible if the species forming the ring regulated rates associated with positive feedback, GTPase inactivation, GAP-mediated negative feedback, GAP activation, or GAP inactivation. However, modulating the basal GTPase activation rate did not generate rosette formation, but instead spots of high GTPase activity formed throughout the entire domain.
We next used the model to demonstrate that rosette patterning could occur in the absence of initial ring formation. This scenario required the following conditions to be met: 1) the positive and negative feedback strengths increased over the IgG disk, 2) negative feedback dominated over positive feedback over the disk and far from the disk and 3) the negative feedback transitioned from its elevated level over the disk to its basal level more rapidly than that of the positive feedback. These three features generated a small region outside of the disk where the positive to negative feedback ratio is sufficiently high to enable spots of active GTPase. This scenario occurs if the positive regulator of GTPase activity diffuses rapidly, whereas the negative regulator of GTPase activity diffuses slowly in comparison.
Further investigations are required to rule out other mechanisms for rosette formation and determine the specific molecular species involved in this process. A starting point for these investigations is the GEFs and GAPs known to a play role in phagocytosis. In particular, the GAPs ARHGAP12, ARHGAP25, and SH3BP1 were found to be important for the successful engulfment of larger targets (beads of diameter = 8.3 μm) but not smaller ones (diameter = 1.6 μm) [47]. However, Schlam et al. suggest that these GAPs are required for the replenishment of a limiting species in actin polymerization, which is not important for small targets [47]. The main particle size we observed (diameter = 3.5 μm) is between their "large" and "small" bead sizes, and cells engaging in frustrated phagocytosis will (by definition) not successfully engulf a target. Thus, these GAPs are good initial targets for further studies, but careful consideration will be required to elucidate their target size-dependent and time-dependent roles.
Another important next step is to obtain estimates for model parameters, as many rates in the model have not been measured in vivo. Our labs recently demonstrated how single molecule measurements can be used to infer the chemical kinetics of conformationally dependent signaling molecules and the number of molecules within clusters in vivo [42]. These approaches should be valuable tools for estimating model parameters and determining important features of Rho GTPases and other podosome-related molecules that the model should capture. Evidence that supports our current choice of parameters is that the model matches the size and spatial organization of podosomes, uses a membrane-bound diffusion rate for Cdc42 that is estimated from experimental data, and uses a reasonable estimate for the number of Cdc42 molecules expected within podosomes. Finally, we note that the time scale predicted by the model for Cdc42 to form rosettes (~20-40s, S4 Fig) appears to be faster than our experimental observation for podosome rosette formation (~200 s) [18]. However, we believe the model's prediction is consistent with that of podosome formation because the recruitment of actin and other podosome-related molecules would occur following GTPase activation.
To investigate the model's behavior, we performed simulations using varying sizes of either IgG disks or holes. Experimental results for the number of spots formed using different IgG disk sizes were consistent with our simulation results. We also noticed that for simulations on disks of radius R = 0.24 μm, 3 distinct spots of GTPase activity were produced, whereas disks of smaller radii produced 1 or 2 sites. This result suggests a threshold for the minimum size of a particle that can be internalized via phagocytosis. If we assume three or more podosomes are required to engulf a target, then the model is consistent with the reported value of 0.5 μm as the minimum size for phagocytic targets [61,62], although it is unclear whether GTPase patterning is the limiting factor in this process. Simulations on holes lacking IgG of different sizes revealed that holes of radii 1.6 μm or less do not form rosettes, while holes of 2.1 μm are capable of rosette formation. Interestingly, for holes between 1.7 μm and 2.0 μm rosette formation depended on initial conditions, suggesting the system is bistable in this regime. Thus, follow up studies exploring frustrated phagocytosis on holes and other micro-patterned shapes of IgG may prove insightful.
In summary, our analysis revealed how relatively minor additions to a Rho GTPase and GAP circuit were sufficient to generate the rosette of podosomes observed during frustrated phagocytosis. It is likely that this same polarity circuit can generate more complex patterns when additional regulatory elements are added, and the computational framework we developed may prove useful for future investigations into the spatial regulation of GTPase activity. Our analysis focused on static cytoskeletal structures. However, the actin cytoskeleton is a dynamic system, and phagocytosis requires exact spatiotemporal control of cellular morphodynamics during engulfment. Therefore, in future studies it will be important to consider in more detail the time-dependent behavior of Rho GTPase signaling during phagocytosis and extend the model to three-dimensions.

Cell culture and transfection
RAW 264.7 macrophages were obtained from the American Type Culture Collection (ATCC) and maintained in culture medium: RPMI 1640 medium GlutaMAX Supplement (Thermo-Fisher Scientific, 61870127) containing 10% heat-inactivated FBS (HI-FBS, GEMINI Bio, 100-106) in a 5% CO 2 humidified incubator at 37˚C. To detach RAW 264.7 cells from the Falcon tissue culture dish (Fisher Scientific, 08-772E), the cells were treated with Accutase (Thermo-Fisher Scientific, A1110501) at 37˚C for 5 min before gentle scraping (CytoOne, CC7600-0220). The plasmids FTractin-tdTomato and myosin regulatory light chain (MRLC)-EGFP were described previously [63,64]. RAW 264.7 cells were electroporated with the Neon Transfection System (ThermoFisher Scientific) following the manufacturer's protocol. In brief, 5x10 6 cells were electroporated with 1 μg plasmid in R buffer at a setting of 1680 V, 20 ms, and 1 pulse using 10 microliter Neon pipette tip. The cells were transferred into a well of 12-well plate, with each well containing 1 ml of culture medium. After 12 hours of incubation, the transfected macrophages were ready for the frustrated phagocytosis experiments.
Bone marrow cells were isolated from 6 to 12 weeks C57BL/6 mice and differentiated into macrophages for 5-7 days in RPMI 1640 medium containing 10% heat inactivated FBS and 10% M-CSF (L929 conditioned medium) described elsewhere [65,66]. These macrophages were detached from the flask using Accutase and gentle scraping.

Microcontact printing
The IgG patterns on glass coverslips were made using the microcontact printing of Polydimethylsiloxane (PDMS) as previously described [67]. The silicon master with an array of 3.5 μm holes spaced 8 μm apart or 10 μm holes spaced 20 μm apart was made using photoresist lithography, and PDMS stamping on glass coverslips was carried out as described previously [18].

Inhibition treatment and immunofluorescence staining
To inhibit actomyosin contractility and disassemble myosin II filaments, RAW 264.7 macrophages were plated on patterned IgG coverslips in Ham's F12 medium (Caisson Labs, UT) supplemented with 2% HI-FBS and 20 μM Rho kinase inhibitor Y-27632 (Hello Bio, HB2297) for 25 min of inhibition during frustrated phagocytosis. For frustrated phagocytosis against 10 μm IgG spots, bone marrow-derived macrophages were plated on patterned IgG in the above medium without inhibitor and incubated at 37˚C for 15 min before staining.
The cells were fixed with 4% paraformaldehyde at 37˚C for 10-15 min and permeabilized using 0.1% Triton-X-100 (Sigma-Aldrich) in PBS for 5 min. Cells were then thoroughly washed with PBS and fixative quenched with 0.1 M glycine for 20 min followed by incubation with 2% BSA fraction V (Thermo, 15260037) in PBS for 30 min. Actin was stained with Alexa-Fluor 568 phalloidin (dilution 1:500, ThermoFisher Scientific A12380) diluted in 2% BSA in PBS at room temperature for 20 min followed by one wash with 1xPBS/0.05% Tween for 10 min, and two washes with 1x PBS for 15 min.

Imaging of podosome structures during frustrated phagocytosis
Total internal reflection fluorescence structured illumination microscopy (TIRF-SIM) was used to image podosomes in F-tractin-tdTomato transfected live RAW 264.7 macrophages. Fluorescence emission was recorded using an sCMOS camera (Hamamatsu, Orca Flash 4.0 v2 sCMOS). Lasers with wavelengths 560 and 647 nm and an Olympus UApo N 100x oil NA 1.49 objective were used, and fluorescence emission was recorded using an sCMOS camera (Hamamatsu, Orca Flash 4.0 v2 sCMOS). A Nikon SIM microscope was used to image podosomes in fixed RAW 264.7 macrophages after Y-27632 inhibition, using 488 and 561 nm lasers. A 100x oil immersion objective (1.49 NA, Nikon CFI Apochromat TIRF 100x) and EMCCD camera (Andor DU-897) were used. To image podosomes in bone marrow-derived macrophages, a Zeiss confocal microscope LSM880 built around AxioObserver 7 with a 63x 1.4 NA oil objective (Zeiss) was used.

Single particle tracking
Single particle tracking was performed using a home-built total internal reflection microscope based on an Olympus IX81. The microscope was equipped with four solid state lasers (Coherent OBIS 405 nm, 488 nm, 561 nm, and 647 nm), a 100X TIRF objective (Olympus, UPLA-PO100XOHR) and an sCMOS camera (Photometrics Prime 95B) for fluorescence collection. Raw cells were co-transfected with mScarlet-F-tractin (Excitation, 561 nm; Emission, Semrock, FF01-600/52) and Cdc42-HaloTag [68]. Note that the HaloTag was attached to the N terminus of Cdc42 to ensure that membrane interactions were not impaired. Cells were incubated with 100 pM dye JF646-Halo (Emission, Semrock, FF01-698/70) for 30 minutes and washed with culture medium three times before imaging. Super-resolved F-tractin images were acquired at 100 Hz for 5 seconds and subjected to Super-Resolution Radial Fluctuations analysis [69]. For single particle tracking of Cdc42, we streamed for 40 seconds at 50 Hz (2000 frames).
Single molecule diffusion analysis was done as before [42]. Briefly, individual molecules were identified by a wavelet decomposition based approach [70] and precise centroids were obtained by fitting with a 2D Gaussian function. Single molecule trajectories were built through a well-established linking algorithm [71] and the mean-square-displacement was then calculated [72,73] to color encode the tracks.

Numerical simulations
Ordinary differential equation (ODE) simulations were performed by using the Python package odeint from Scipy [74]. Reaction-diffusion equations were solved using the spectral differential equation solver Python package Dedalus [75]. For simulations using Cartesian coordinates, the system was spatially discretized using a Fourier basis in x and a Chebyshev basis in y with the recommended dealiasing factor of 1.5, as done before [10]. The system had periodic boundary conditions in x and Neumann (reflective) boundary conditions in y. Similarly, for simulations using polar coordinates, the system was spatially discretized using a Fourier basis in ϕ, and a Chebyshev basis in r (dealiasing factor of 1.5) with periodic boundary conditions in ϕ, and Neumann (no flux) boundary conditions in r. The typical grid size used for simulations was 256 x 128 (ϕ, r respectively) which was informed by mesh grid refinement (i.e., larger grid sizes resulted in the same outcome). However, a grid size of 64 x 64 was used for parameterization steps to decrease simulation time. Simulations were typically performed using a time step dt = 0.01 s or 0.025 s. Reaction steps were solved using 4 th order Runge-Kutta, although 2 nd order Runge-Kutta was used for parameterization steps.
Homogeneous steady states were determined by running the ODE system (without diffusion) for t = 1000 s using odeint. For initial conditions of the reaction-diffusion equations, each species was set to its steady state value throughout the domain and subsequently noise was added by converting a small fraction of inactive species to the active form. The fraction of concentration converted was determined by each simulation but was typically generated by uniform sampling between 0 and 0.2v ss (where v ss is the steady state concentration for the inactive species). For seeded simulations, the same random noise (but between 0 and 0.1v ss ) was converted and additionally the normalized seed (i.e., a rosette) was scaled by 0.1v ss and converted to active GTPase. For simulations with non-constant coefficients, the initial steady states were determined by using the basal values for the spatially-dependent rates.
To fit the logistic equation to the simple reaction-diffusion profiles, we used the Python package minimize (using method = "SLSQP") from Scipy [74].
Note that for the coupled model, where the simple reaction diffusion model was simulated in place of the logistic function, we used the same MCMC pipeline for sampling to discover a working coupled model. For proof of concept, we simply ran this pipeline for 1,000 steps and took the best scoring parameter set.

Spot size determination
Simulations were performed as described above using cartesian coordinates (t final = 40 s). Each system was interpolated to a uniform grid with the same grid size (128 x 128). A mask was generated by thresholding at the mean of the maximum and minimum concentrations within the system. Using this mask, features were quantified using the Python package scikit-image [78]. The effective radius was defined as the average between the major and minor axis lengths. After an initial pass at fitting Radius eff based on log(D u ), simulations were rerun using different grid sizes to ensure the number of spots counted for each simulation were similar (N~100).

Counting podosomes per site
For experimental results, the number of podosomes per site were calculated using a pipeline we developed previously [18] (https://github.com/elstonlab/PodosomeImageAnalysis). In essence, this pipeline uses persistent homology, a type of topological data analysis, to identify significantly persistent features (connected components, holes) within images followed by post-processing.
For simulated results, a mask was created by thresholding at the average between the maximum and minimum intensity within a simulation. From this mask, the number of features was counted using the Python package ndimage.label from Scipy [74].
Supporting information S1 Fig. Additional experimental evidence for GTPases in frustrated phagocytosis. A) Control RAW 264.7 macrophages were marked for actin (phalloidin staining) and myosin II (RLC-eGFP) during frustrated phagocytosis. B) RAW 264.7 macrophages were marked for actin (phalloidin staining) and myosin II (RLC-eGFP) when treated with 20mM Rho kinase inhibitor Y-27632 for 25 min during frustrated phagocytosis. Proper actin rosette formation despite Rho kinase inhibition suggests that actomyosin contractility is not necessary for rosette formation. C) Diffusion coefficient estimates for Cdc42 during frustrated phagocytosis. Quantifications from the mean squared displacement analysis of the single particle tracking performed in Fig 1D-1G Table 2). The worst scoring parameter set shown by magenta diamonds (S1 Table). E) Active GTPase concentration for the worst scoring parameter set (S1 Table). from a coupled model, where the diffusing species are simulated (dashed lines, see S1 Note) and the logistic function approximations using the representative parameter set ( Table 2). B) Active GTPase concentration for the coupled model shown in A. C) Spatial concentration of the diffusing species modulating the GAP activation rate, c, in the coupled model in A,B. D) Spatial concentration of the diffusing species modulating the self-positive feedback rate, γ, the coupled model in A,B. E) The positive to negative feedback ratio γ/c forms a ring. For visualization purposes, the GAP activation rate, c, is multiplied by 100 to be a similar order of magnitude to the self-positive feedback rate γ. (TIF)

S6 Fig. Additional information for experimental and simulated results on varying disk and hole sizes. A)
Representative experimental results for disks of radius 1.75 μm. Podosomes are indicated with red circles. B) Same as A but using disks of radius 5 μm. C) Simulated number of GTPase spots versus disk radius for small disks with radius less than 0.25 μm. D) Simulations for the representative parameter set (Table 2, negative c km and γ km ) when changing the hole size. (TIF) S1 Table. Worst scoring parameter set after running the MCMC. Fixed values same as in Table 2. (DOCX) S1 Note. Verification of a fully coupled model. Details for using reaction-diffusion equations to determine spatially-dependent rates rather than the logistic function approximation.