Diffusion-limited phase separation in eukaryotic chemotaxis

Results and Discussion

We simulated the kinetics of the network of chemical reactions that represents the ubiquitous biochemical backbone of the directional sensing module and contains the following: (i) binding of PI3K to activated membrane receptors, (ii) binding of PTEN to PIP₂, (iii) catalytic activity of PI3K and PTEN, and (iv) phosphoinositide diffusion within the plasma membrane. Because the chemical system is characterized by extremely low concentrations of chemical factors [0–50 nM for enzymes (10) and 0.1–1 μM for phosphoinositides (11)] and evolution takes place out of equilibrium, we used a stochastic approach (8, 9). Indeed, rare, large fluctuations are likely to be relevant for kinetics in the presence of unstable or metastable states. This assumption is also consistent with the observation that cell response to chemotactic stimuli has a stochastic character (12). Simulated reactions and diffusion processes are summarized as follows:

1. PI3K(cytosol) + Rec(i) ↔ PI3K·Rec(i)

2. PTEN(cytosol) + PIP₂(i) ↔ PTEN·PIP₂(i)

3. PI3K·Rec(i) + PIP₂(i) → PI3K·Rec(i) + PIP₃(i)

4. PTEN·PIP₂(i) + PIP₃(i) → PTEN·PIP₂(i) + PIP₂(i)

5. PIP₂(i) → PIP₂(j)

6. PIP₃(i) → PIP₃(j)

Index i represents a generic plasma membrane site and j one of its nearest neighbors. We have simulated the chemical kinetics on the inner face of the plasma membrane as a two-dimensional lattice gas coupled to an unstructured cytoplasmic reservoir. On each site i of the lattice we have assumed the presence of a number of molecules of the k-th species. Chemical reactions and diffusion have been simulated as random processes with intensities proportional to kinetic reaction and diffusion rates (see Tables 1 and 2).

The probability of performing a reaction on a given site is taken to be proportional to realistic kinetic reaction rates and local reactant concentrations (Tables 1 and 2).

Stochastic Simulation. The plasma membrane is represented as a sphere of radius R = 10 μm partitioned in N _s = 10,242 computational sites forming a honeycomb lattice with 12 pentagonal defects. Chemical concentrations are represented as integer variables giving the number of molecules of the k-th species present on site i. Reaction–diffusion kinetics is simulated according to Gillespie's method (8) generalized to the case of an anisotropic environment. Catalytic processes are described by Michaelis–Menten kinetics. The density of activated receptors is proportional to extracellular chemoattractant concentration. The probability of diffusion from a computational site to a neighboring one is assumed proportional to the difference in local concentrations, according to Fick's law. Time is advanced as a Poisson process. To provide complete reproducibility, the C code of the simulation algorithm is included as Data Set 1.

The density ρ of activated receptors controls PI3K recruiting to the plasma membrane, thus playing the role of a chemical potential (13). Although receptor activation is directly controlled in the experiments through the extracellular concentration of chemoattractant stimulus, it is difficult to control experimentally [PI3K]. We have therefore set [PI3K] to a realistic fixed value (Table 2) and varied ρ in the range 0–100 nM. Diffusion controls the behavior of the system in an obvious way, because for high values it tends to mix chemical species, but in conjunction with catalytic activity it also exerts an ordering action, transferring information between neighboring uniformly populated patches of the plasma membrane. Diffusion is therefore a second parameter, which strongly influences the system's dynamic and stationary state. We have considered D values in the range 0–5 μm²/s.

Order Parameter. Because PTEN localizes on PIP₂-rich regions of the plasma membrane, phase separation can be observed either at the enzyme or phosphoinositide level. In real cells, phosphoinositide clusters trigger actin polymerization and motility. For this reason we have studied symmetry breaking in phosphoinositide distribution. The physics of phase-separation dynamics has been thoroughly studied (14–18) and is known to give rise to a variety of effects such as self-organization, pattern formation, and pattern selection in many physical–chemical systems (refs. 6, 19, and 20; for a general reference on pattern formation in reaction–diffusion systems, see ref. 21). The degree of order of a chemical mixture undergoing phase separation can be quantified by means of an order parameter, i.e., a dimensionless observable assuming the value 0 in the symmetric, mixed state and a value of order 1 in the symmetry-broken state in which chemical species are separated (13). A convenient order parameter measuring the degree of phase separation of the phosphoinositide mixture is Binder's cumulant (9, 22) where φ = φ _i = [PIP₃] _i - [PIP₂] _i is a difference of local concentrations on site i, and 〈... 〉 denotes average over many different random realizations. Binder's cumulant measures the distance of the probability distribution of φ from a Gaussian distribution. When phosphoinositides are mixed, fluctuations around the uniform average value are Gaussian, and g tends to zero, whereas when phosphoinositides separate into distinct clusters (Fig. 1), the probability distribution of φ is characterized by two distinct peaks, and g becomes of order 1 (Fig. 1 d and e ). Maximal observed values of g for phase-separating systems are between 0.4 and 0.8, corresponding to the fact that the PIP₂ and PIP₃ distributions partially overlap.

Spontaneous phase symmetry breaking leads to the formation of PIP₂, PIP₃-rich clusters of different size. Cluster sizes can be characterized by harmonic analysis. For each realization, the fluctuations δφ = φ - 〈φ 〉 of the φ field can be expanded in spherical harmonics (23) as where u is a unit vector identifying a point on the spherical surface. Let us consider the two-point correlation functions where 〈... 〉_iso denotes average over ensembles (7) and over the action of the sphere rotation group, , and P_l are Legendre polynomials (23). A measure of the average cluster size is πR/2 〈l 〉, where . In particular, if most of the weight is concentrated on the l-th harmonic component, average phosphoinositide clusters extend over the characteristic length πR/2l.

When receptor activation ρ = ρ _i = [Rec] _i is not uniformly distributed we are interested in its correlation with localized PIP₂, PIP₃ clusters, which is measured by the components of the covariance matrix c _ρρ = 〈(ρ- 〈ρ 〉)² 〉, c _ρφ = 〈(ρ- 〈ρ 〉)(φ- 〈φ 〉) 〉, c _φφ = 〈(φ- 〈φ 〉)² 〉 and by the correlation coefficient (24).

Dynamic Phase Diagram. We have run ≈10 random realizations of the system for ≈400 (ρ_R, D) pairs with 1 μM < ρ < 100 μM, 0.1 μm²·s^-1 < D < 5 μm²·s^-1 for 2 h of simulated time [≈1 Tflop (trillion floating point operations) per realization]. For each (ρ, D) pair we have computed g as a function of time by performing surface and ensemble average at fixed time intervals. For each random realization we have started from a stationary homogeneous PTEN, PIP₂ distribution (10 s of simulated time evolution were enough to reach stationarity). To measure automatically phase-separation events and characteristic phase-separation times, we have computed the 10-min running average of g to cut off rapid fluctuations and selected the moment when it reached the threshold 0.4 and did not fluctuate below that value during the following 30 min of simulated time. Patterns observed at t = 2 h are then quite close to stationarity. Similarly, we have measured the harmonic components C_l . In this case we have identified the phase-separation time with the moment when the weight of the 10-min running average of a single harmonic component surpassed the 80% of the total weight and did not fluctuate below that value for the following 30 min of simulation time.

We let the system evolve to stationarity, in the absence of receptor activation, obtaining a homogeneous PTEN and PIP₂ distribution at the plasma membrane, identical to that observed in unstimulated cells. At time t = 0 receptor activation is switched on; either activated receptors are isotropically distributed or the isotropic distribution is perturbed with a linear term producing a 5% difference in activated receptor density between the North and the South poles.

In the isotropic case, we found that in a wide region of parameter space, the chemical network presents an instability with respect to phase separation (25, 26), i.e., the homogeneous phosphoinositide mixture realized soon after receptor activation is unstable and tends to decay into spatially separated PIP₂- and PIP₃-rich phases (Fig. 1 a–d ). Characteristic times for phase separation vary from the order of a minute to that of an hour, depending on receptor activation (Fig. 1e ). The dynamic behavior and stationary state of the system strongly depend on the values of two key parameters: the concentration ρ of activated receptors and the diffusivity D.

In the case of anisotropic stimulation, orientation of PIP₂ and PIP₃ patches clearly correlates with signal anisotropy (Fig. 1 g–i ), and, compared with isotropic stimulation, phase separation takes place in a larger region of parameter space and in times that can be shorter by one order of magnitude (Fig. 1j ).

Average phase-separation times are plotted in Fig. 2a , where red areas correspond to non-phase-separating systems. In the deep blue area, phase separation takes place in <5 min, whereas close to the boundary of the broken symmetry region phase separation can take times of the order of an hour (Fig. 1e ). Average cluster sizes at stationarity are plotted in Fig. 2b . In the red region, cluster sizes are of the order of the size of the system, corresponding to the formation of pairs of complementary PIP₂ and PIP₃ patches (Fig. 1). For diffusivities <0.1 μm²/s the diffusion-mediated interaction is unable to establish correlations on lengths of the order of the size of the system, and one observes the formation of clusters of separated phases of size much smaller than the size of the system (Fig. 2c ). For diffusivities >2 μm²/s the tendency to phase separation is contrasted by the disordering action of phosphoinositide diffusion (Fig. 2 a–c ).

Average phase-separation times for the anisotropic case are plotted in Fig. 2d . By comparing with the isotropic case (Fig. 2a ), it appears that there is a large region of parameter space where phase separation is not observed with isotropic stimulation, while a 5% anisotropic modulation of activated receptor density triggers a fast phase-separation process. Cluster sizes (Fig. 2e ) are in the average larger in the anisotropic case than in the isotropic case (Fig. 2a ). Fig. 2f shows that the orientation of PIP₂ and PIP₃ patches is strongly correlated with signal anisotropy. Therefore, anisotropy has two main effects: on one hand, it triggers the phase-separation process on shorter timescales and in a wider region of parameter space; on the other hand, its direction breaks the system's rotational symmetry and selects the stationary state of the system. For receptor activation ρ < 20 nM, the few PI3K enzymes bound to the plasma membrane are not sufficient to create extended PIP₃-rich regions; however, small intermittent PIP₃ clusters are still observed (Fig. 3a ).

Simulating a gradient in receptor activation similar to the one imposed in experimental assays, we are able (Fig. 4) to reproduce the input–output relationship observed on Dictyostelium (4). It is also worth noting that the characteristic timescales for phase separation emerging from our dynamical simulations, which have been performed by using realistic reaction and diffusion rates, are in agreement with the observed ones (4).

Physical Picture. The transition from a phase-separating to a phase-mixing regime results from a competition between the ordering effect of the interactions and the disordering effect of molecular diffusivity. The frontier between these two regimes varies continuously as a function of parameters. Importantly, we found that the overall phase-separation picture is robust with respect to parameter perturbations, because it persists even for concentrations and reaction rates differing from those of Table 2 by one order of magnitude. Moreover, both in isotropic and anisotropic conditions, signal amplification is completely reversible. Switching off receptor activation abolishes phase separation, delocalizes PI3K from the plasma membrane to the cytosol, and brings the system back to the quiescent state (see Movie 1, which is published as supporting information on the PNAS web site).

Thus, phosphoinositide diffusion is directly responsible for establishing correlations between neighboring sites, leading the system to a phase-separation instability. Although large diffusivity has a mixing effect, intermediate diffusivity cooperates with catalysis to order the system on large scales. Physically, this process can be understood as follows. Receptor activation shifts the chemical potential for PI3K, which is thus recruited to the plasma membrane. PI3K catalytic activity produces PIP₃ molecules from the initial PIP₂ sea. Initially, the two phosphoinositide species are well mixed. Fluctuations in PIP₂ and PIP₃ concentrations are, however, enhanced by preferential binding of PTEN to its own diffusing phosphoinositide product, PIP₂. Binding of a PTEN molecule to a cell membrane site induces a localized transformation of PIP₃ into PIP₂, resulting in higher probability of binding other PTEN molecules at neighboring sites. Such a positive feedback loop not only amplifies the inhibitory PTEN signal, but via phosphoinositide diffusion it also establishes spatiotemporal correlations strong enough to drive the system toward spontaneous phase separation. The time needed by the system to fall into the more stable, phase-separated phase, however, can be a long one if the symmetric, unbroken phase is metastable. In that case, a small anisotropic perturbation in the pattern of receptor activation can be enormously amplified by the system instability.

It is worth noting that, when the system phase separates, the final size of the clusters is limited only by the size of the system and the availability of cytosolic enzymes. In an infinite system, clusters would grow indefinitely. In a finite system instead, as long as PTEN molecules are recruited on the PIP₂ patch, they are no longer available to compete with PI3K on the residual PIP₃ patch, which is therefore stabilized. The net effect is that the cluster size saturates at a stationary value of the order of the size of the system.

Dimensional considerations suggest that the size of the patches should grow proportionally to , 1/k being the characteristic enzyme association–dissociation timescale, for low diffusivity values, and saturate to the size of the system for intermediate diffusivities. Simulations confirm the dimensional estimate and show that for higher diffusivities the species are mixed up and the size of clusters drops abruptly (Fig. 5).

The properties of the real reaction and diffusion processes described in Results and Discussion and Table 1 are better understood through the study of a one-dimensional model derived from them under simplifying assumptions (see Supporting Text and Fig. 6, which are published as supporting information on the PNAS web site), where a more complete analysis is possible. The simplified model presents a parameter region where multiple stable equilibria are possible, and stochasticity in the number of membrane-bound enzymes can trigger a transition from a less stable to a stabler state. The transition takes place through the formation of a small region of the stabler phase in the sea of the less stable one. Although small regions of the stabler phase are wiped off by diffusion, larger regions propagate with finite velocity in the less stable phase and would eventually take over the whole system if the number of available cytosolic enzymes was infinite. Because, however, the number of cytosolic enzymes is limited, the process slows down and eventually stops with the formation of a stationary front separating a PTEN-rich and a PI3K-rich region. The critical size of the nucleating region is determined by the relative strength of two intrinsically dynamic quantities: diffusivity and the velocity of front propagation.

Conclusions

Our results provide a simple physical clue to the enigmatic behavior observed in eukaryotic cells. As we have shown, there is a large region of parameter space where the cell can be insensitive to uniform stimulation over very large times but responsive to slight anisotropies in receptor activation in times of the order of minutes (Fig. 1 and Movie 1). Accordingly, by simulating shallow gradients of chemoattractant, we observed PIP₃ patches accumulating with high probability on the side of the plasma membrane with higher concentration of activated receptors, thus resulting into a large amplification of the chemotactic signal (Figs. 1 and 4). Moreover, we identified an intermediate region of parameters, where phase separation under isotropic stimulation is observed on average in a long but finite time. In this case, one would predict that on long timescales cells undergo spontaneous polarization in random directions and that the number of polarized cells grows with time. Intriguingly, this peculiar motile behavior is known as chemokinesis and is observed in cell-motility experiments when cells are exposed to chemoattractants in the absence of a gradient (27). An additional consequence of the tendency to phase separation is that for low values of receptor activation small intermittent clusters should form, because diffusion cannot establish a correlation between too-distant phosphoinositide patches and the number of enzymes producing the patches is too low. Indeed, intermittent phosphoinositide clusters are observed in our simulations (Fig. 3a ), and their formation has been recently described in Dictyostelium cells exposed to very low cAMP concentrations (28). It appears that D3-phosphoinositide patches serve as a spatial cue for pseudopod formation, which enhances the sensitivity and amplitude of chemotactic movement. Furthermore, cells migrating within tissues, such as neurons or immune system cells, may encounter multiple chemoattractant signals in complex spatial and temporal patterns that can potentially direct their path. Notably, cell polarization induced by multiple chemoattractant sources has been observed in vitro (4), and this experimental situation can be mimicked in silico by simulating the receptor activation pattern produced by multiple chemoattractant sources (Fig. 3b ).

In summary, the phase-separation scenario provides a simple and unified framework to different aspects of directed cell motility, such as large amplification of slight signal anisotropies, insensitivity to uniform stimulation, appearance of isolated and transient phosphoinositide patches, and stochastic cell polarization. It provides a link between known microscopic and macroscopic timescales. Finally, it unifies apparently conflicting aspects that previous modeling efforts (29, 30) could not satisfactorily reconcile (31), such as insensitivity to absolute stimulation values, large amplification of shallow chemotactic gradients, reversibility of phase separation, robustness with respect to parameter perturbations, stochastic character of cell response, use of realistic biochemical parameters, and space–time scales.

One of the characteristic features of our model is that it brings together into a unified picture several seemingly disconnected phenomena, such as response to anisotropic chemotactic stimulation, stochastic polarization, and small cluster formation for low activation levels. It would be useful to perform systematic quantitative observations on stochastic cell polarization under uniform stimulating conditions, both at high and low activation leves. Because nucleation can be arguably modeled as a Poisson process, in the presence of isotropic stimulation one should observe a number of stochastically polarized cells growing with time, with a rate comparable with the model predictions.