Coarsening of topological defects in 2D polar active matter

We study the following model of 2D polar active matter.

Eq.1 describes the evolution of the 2D polar field $\mathbf{p}\equiv(p_{x},p_{y})$, coupled to a concentration field $c$. Eq.2 is a continuity equation ensuring the conservation of the total number of active particles transported through active ($\mathbf{j_{a}}$) and diffusive ($\mathbf{j_{d}}$) currents. The first two terms on the right hand side of Eq. 1 arises from the free energy functional of the 2D XY model [22]. The third, active, term couples $\mathbf{p}(\mathbf{r},t)$ to $c(\mathbf{r},t)$ through the activity parameter (called contractility), $\zeta$. In Eq.2, the active term couples the two fields through the self-propulsion speed $v_{0}$, which we take to be 1. This model is useful in studying many biological systems, such as the actin cytoskeleton, where the concentration and the orientation of the active constituents are coupled to each other [22]. We simplified Eq. 1 further by taking advection coefficient $\lambda=0$ and noise $\mathbf{f_{p}}=\mathbf{0}$, which does not change the generality of the results. Simulation parameters and methodology are described in the SI.

To study the steady state behavior of the model, we numerically solved the equations for $10^{7}$ iterations ($dt=10^{-3}$) for different values of $\zeta$, keeping all other parameters fixed. We initialized the equations with $\mathbf{p}(\mathbf{r})=(1,0)$ for all $\mathbf{r}$ and used a disordered initial condition for $c(\mathbf{r})$. The initial configuration rapidly evolved to a state with many defects, which then coalesced with each other to converge to a final state, where the number of defects, $n_{0}$, did not change appreciably over several million iterations (Fig. 1a-d). $n_{0}$ is nonzero for all $\zeta$ values. However, the approach to $n_{0}$ is different above and below $\zeta=5$ (Fig. 1d-e). Below this threshold, the initial ordered state becomes linearly unstable through laminar structures at wave numbers $\mathbf{q_{d}}$, whose magnitude $q_{d}$ scales as $\zeta^{0.5}$. For $\zeta>5$, the ordered state transforms to a lattice of asters with each row alternating between inward (in-) and outward (out-) pointing asters. Because of the coupling of the $c$ field with the $\mathbf{p}$ field, $c$ is depleted from the out-asters and enhanced at the in-asters (SI movies). Hence, local peaks in $c$ field correspond to the location of the in-asters. After the initial transient, the nonlinear effects become important, which leads to the nucleation of asters followed by their coalescence for $\zeta\leq 5$, leading to a nonmonotonic evolution of asters towards $n_{0}(\zeta)$. For $\zeta>5$, the number of asters decays to $n_{0}(\zeta)$ through coalescence of smaller asters to larger asters. The zeta-dependence of $n_{0}$ is nontrivial. For $\zeta\leq 5$, $n_{0}\sim\zeta^{1}$, which agrees with the prediction from the linearized theory, because the nucleation of the defects happen along the laminar structures, implying the scaling $n_{0}\sim q_{d}^{-2}$. In contrast, for $\zeta\geq 20$, $n_{0}\sim\zeta^{1/3}$, implying that the minimum separation between the defects scales as $\zeta^{1/6}$ (SI). The effect of nonlinearity becomes important for $5\leq\zeta\leq 20$, where $n_{0}$ shows smooth transition between the two regimes. The difference in the steady states in the different regimes is readily observed in the defect arrangements, which organize along 1D channels for $\zeta\leq 5$ and in liquid-like 2D structures for $\zeta\geq 20$, showing intermediate arrangements in between (SI movies).

The dynamics of defects is fundamentally different from its equilibrium counterpart. In equilibrium, at zero temperature, all defects are annealed out to lower the free energy [4]. The annealing kinetics is slow and shows logarithmic relaxation [8]. For example, below $T_{KT}$, the defect density, $n(t)$, scales as $n(t)\ln n(t)\sim t^{-1}$, which goes to zero as $t\rightarrow\infty$ [8]. The relaxation kinetics of $n(t)$ in our model is clearly different from thermal equilibrium. To make the difference in kinetics evident, we compare how $n(t)$ varies with $\Upsilon(t)=-(n^{-2}dn/dt)^{-1}$ [8]. When the equilibrium kinetics hold, $n(t)\sim e^{-A\Upsilon}$, where $A$ is a system-dependent constant. In contrast, in our model, we find that $n(t)-n_{0}\sim\Upsilon^{-1}$, which confirms the fundamental difference between the active and the passive cases.

To get a more microscopic understanding of the coarsening kinetics, we studied the coalescence of two in-asters. In between every two asters, there is a saddle point, which carries -1 topological charge (Fig 3). The saddle point connects two out-asters along a separatrix that separates the basin of attraction of the two in-asters (Fig 3a). The coalescence is facilitated through a large-scale reorganization of the $\mathbf{p}$ field. In the immediate neighborhood of the two in-asters of interest, this reorganization results in the disappearance of the separatrix and the overlap of their basin of attractions (step 1, Fig 3b). After the overlap, one of the in-asters rapidly annihilate the saddle (step 2, Fig 3c-d), and the surrounding material relaxes to a rotationaly symmetric structure around the remaining defect (step 3,Fig 3e). What is surprising is that this process does not continue indefinitely. Once the defect density reaches $n_{0}$, the defects do not coalesce any more. Our initial hypothesis was that these are kinetically arrested states and would require longer simulations to continue the coalescence. However, even after running the simulation ten times longer than required to reach the steady state (Fig. 1), we could not observe a single coalescence event. This observation is surprising because, in equilibrium, defects interact with each other through 2D coulomb interactions, which implies that at $T=0$ the coalescence should continue until no defects are present. This observation, therefore, is consistent with the hypothesis that activity screens the effective interaction between the defects. The screening in the presence of activity is inherently different than that observed for $T>T_{KT}$ in equilibrium, where proliferation of many defects leads to the screening of the interaction and maintenance of the defect plasma. In contrast, here, the screening increases with $\zeta$, implying that the nonreciprocity of the active terms may be the origin of the screening, which has also been shown in the absence of number conservation ¹¹1Unpublished data from Pankaj Popli and Sriram Ramaswamy. Understanding the active screening process requires the equation of motions of defects in 2D polar active materials, which are currently unavailable. Therefore, we resort to an indirect confirmation of this hypothesis by measuring the aster size distributions.

We have defined the size of the basin of attraction (BoA) of an aster as the size of the aster. However, true BoA is hard to calculate for our system and a simpler approximation must be found to generalize across the different system sizes and the parameters that we explore here. To this end, we have found that a disk centered on the aster core, whose radius is the distance to the nearest saddle, is a good approximation of the BoA and is our measure of the size, $A$, of the aster (Fig. 4a). The aster size distribution, $p(A)$ shows lognormal tails with power law head (Fig. 4b). We can quantitatively reproduce $p(A)$ from a model of aggregation which only assumes that the aster interactions are screened beyond a cutoff distance $R_{c}$ and that the asters are immobile (SI). An aster increases its size by coalescing with the nearest aster until it exhausts all asters within $R_{c}$ or gets absorbed into another aster. For a wide range of $R_{c}$, $p(A)$ obtained from the simulation and the model are identical except at the tail (Fig. 4b), which bolsters our hypothesis that the equilibrium long-ranged interactions between the asters is screened in the presence of activity.

In experiments, such as those with actomyosin complexes, identifying BoA is difficult, if not impossible. In contrast, the intensity of materials around the defects can be easily measured through optical microscopes, which is what we define here as the intensity of asters. Among all the defects considered here, the in-asters are the sinks, whereas the out-asters are the source. Hence, the $c$-field concentrates around the in-asters. Therefore, by integrating the $c$ field in the disk to measure $A$, we can calculate the intensity of the in-asters. Because $c\mathbf{p}$ vanishes exponentially outside the circle, the quantity of materials present outside the disk is negligible (Fig. 4c). Hence, using our method, we can easily and accurately calculate the aster intensity for various parameters. The mean aster intensity $\langle I\rangle\sim n_{0}^{-1}$. Surprisingly, the aster intensity distribution, $p(I/\langle I\rangle)$, itself shows nontrivial transition with $\zeta$ (Fig. 4d). For $I_{sc}=I/\langle I\rangle<1$, $p(I_{sc})\sim I_{sc}^{k(\zeta)}$. The power law exponent, $k(\zeta)$, interestingly, scales as $k\sim(\zeta-\zeta_{0})^{1/2}$ for $\zeta>\zeta_{0}\approx 12$, and as $(\zeta-\zeta_{0})^{0}$ for $\zeta<\zeta_{0}$²²2The scaling of $k$ with $\zeta$ is suggestive and requires more careful investigation to ascertain the scaling exponent.. The observed scaling indicates that the power law exponent is the order parameter of some underlying transition, which we cannot ascertain here. The value of $\zeta_{0}$ is also unclear due to noise in the data. A comprehensive investigation of this transition is beyond the scope of this manuscript.

Experiments on the actin cytoskeleton provides direct experimental confirmation of our results. For example, in vitro experiments on actomyosin complexes show coalescence of actin asters that are visually similar to our model: asters coalesce rapidly, the smaller asters get absorbed by larger asters, and the asters never coalesce into one single aster [25, 26]. We expect similar dynamics to be present in live cells also. In fact, measurement of residence times of peripheral membrane proteins on the plasma membrane provides indirect evidence of the heterogeneous distribution of actin asters [24]. To make this observation more concrete, we compared the steady state distribution of aster intensities from our model with that obtained from live cell experiments on mouse embryonic stem cells (mESCs) [27]. To make exact comparison, we converted the aster densities in our model to experimental units and compared it with the experimental aster densities (Fig. 5a). For untreated cells and for cells treated with the drug cytochalasin D (CytoD) the average density of the asters was approximately $0.2\mu m^{-2}$. In our model, the aster density at $\zeta=60$ is comparable to the above experiments, and we found that the aster intensity distribution for this $\zeta$ matches almost exactly with the experimental data (Fig. 5b). Treating the mESCs with the drug SMIFH2 reduces the aster density to $0.15\mu m^{-2}$, which corresponds to $\zeta$ between 20 and 30 in our model. We found that $p(I)$ for $\zeta=20$ matches excellently with the experimental data (Fig. 5c), which predicts that SMIFH2 reduces the density of asters through the reduction of contractility. This observation is surprising because SMIFH2 is an inhibitor of Formin, a catalyst for actin polymerization, whereas contractility arises due to the effect of the motor protein myosin on the ACS. Interestingly, recent experiments suggest that SMIFH2 also inhibits myosin [30], which may explain our observation. Finally, we compared our model with CK666-treated cells, for which $\zeta=5$ seemed the best match. However, the predicted $p(I)$ distribution did not match well with the experimental data for $\zeta=5$, but matched well with the distribution for $\zeta=30$ (Fig. 5d, see SI for the comparison). This mismatch suggested that CK666 lowers aster density not by lowering the contractility of the ACS. Indeed, CK666 is an inhibitor of the actin nucleating proteins Arp2/3, which has been shown to regulate the ACS architecture [27]. Therefore, our model needs to be modified appropriately to understand the effect of Arp2/3-dependent aster nucleation.

From the experimental data [26], it is evident that the asters never coalesce into a single aster even after many hours, which are orders of magnitude longer than the maximum timescales of our simulations. Therefore, it remains unclear whether the heterogeneous multi-aster state is a steady state or an arrested state. It is an important theoretical question, which will help us characterize the coarsening of topological defects in active polar system and contrast it with the equilibrium systems. Specifically, it was recently shown that a spontaneous proliferation of defects into multi-aster state destroys order in 2D active polar fluids without number conservation [31]. In that context, if the multi-aster state is an arrested state, we expect the order of the flocking state to be restored at long times. Whereas if the multi-aster state is a true steady state, then the ordered state is truly metastable. For our model, understanding this question will be important to understand the coarsening kinetics. In particular, we need to confirm whether the power law relaxation to the final state is a true relaxation to the steady state. The details are important and remains to be worked out in the future.

Pragmatically, the distinction between an arrested metastable state and a steady state should be judged based on the relevant timescales of the problem. For a biological system, the important timescales are not the timescale to reach the final state, but the timescales of signaling, which are often much shorter than the time to relax to the final state. Similar considerations may apply even when we consider a material made out of active polar filaments, where the timescales of perturbations maybe much shorter than the relaxation timescales. As we have shown, there is excellent correspondence between experiments on actomyosin cortex and our model. Therefore, various predictions from our model, such as the phase transition in the intensity distribution, can be easily tested in existing experimental systems. The theory itself will be useful to understand general 2D active polar systems. Especially, the dynamics of multiple topological defects will prove to be useful to understand the mechanical and the rheological properties of these materials. These questions will be addressed elsewhere.