A minimal model for structure, dynamics, and tension of monolayered cell colonies
AA minimal model for structure, dynamics, and tension of monolayered cell colonies
Debarati Sarkar, ∗ Gerhard Gompper, † and Jens Elgeti ‡ Theoretical Physics of Living Matter, Institute of Biological Information Processing and Institute for Advanced Simulation,Forschungszentrum J¨ulich, D-52425 J¨ulich, Germany
The motion of cells in tissues is an ubiquitous phenomenon. In particular, in monolayered cellcolonies in vitro, pronounced collective behavior with swirl-like motion has been observed deepwithin a cell colony, while at the same time, the colony remains cohesive, with not a single cellescaping at the edge. Thus, the colony displays liquid-like properties inside, in coexistence witha cell-free “vacuum” outside. How can adhesion be strong enough to keep cells together, while atthe same time not jam the system in a glassy state? What kind of minimal model can describesuch a behavior? Which other signatures of activity arise from the internal fluidity? We proposea novel active Brownian particle model with attraction, in which the interaction potential has abroad minimum to give particles enough wiggling space to be collectively in the fluid state. Wedemonstrate that for moderate propulsion, this model can generate the fluid-vacuum coexistencedescribed above. In addition, the combination of the fluid nature of the colony with cohesion leadsto preferred orientation of the cell polarity, pointing outward, at the edge, which in turn gives riseto a tensile stress in the colony – as observed experimentally for epithelial sheets. For strongerpropulsion, collective detachment of cell clusters is predicted. Further addition of an alignmentpreference of cell polarity and velocity direction results in enhanced coordinated, swirl-like motion,increased tensile stress and cell-cluster detachment.
Keywords: collective cell migration | liquid-vacuum coexistence | tensile stress | cellular velocity alignment | coordinated motion I. INTRODUCTION M any fundamental biological processes, like embryo-genesis, wound healing or cancer/tumor invasion requirecells to move collectively within tissues [1–3]. The physicsunderlying these processes ranges from understandingactin polymerization and tread-milling for force gener-ation [4, 5] and single cell migration [6, 7], to the collec-tive behavior of many migrating cells [8–10]. Here, wefocus on an observation from monolayers of migratingMadin-Darby canine kidney (MDCK) cells on surfaces,a prototypical model system for collective cell migra-tion. Experimental observations reveal large-scale collec-tive motion, like swirls, within the bulk of young mono-layers [11, 12], thus the display of fluid-like properties,before jamming occurs as the epithelial sheet matures[13–16]. Interestingly, as an initial colony expands, nocells detach from the boundary — even though the bulkof the tissue remains clearly fluid-like [9]. Cohesion isstrong enough that fingers of many cells can protrudeat the propagating tissue front without cell detachment.Even stronger-pulling “leader cells” do not detach [17–20]. Cells are thus in a ’liquid-vacuum’ coexistenceregime. Even more surprising, pioneering experimentshave revealed that these expanding colonies are undertensile stress [22, 23]. This raises the question how thisliquid-vacuum coexistence, in combination with strong ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] We use the term ’vacuum’ here somewhat loosely to represent aphase of extremely low cell density collective motion and tensile stress, can be captured andunderstood from a minimal physical model.An active Brownian particle (ABP) model [24–27] forcells with standard attractive Lennard-Jones (LJ) inter-actions has been proposed to study cell colonies [28–30]. However, only solid-vacuum (no fluidity of the con-densed phase), or liquid-gas (finite cell density in the di-lute phase) coexistence has been obtained. The coexis-tence of liquids with a very-low-density gas phase is ofcourse well known in many equilibrium systems. In thebiological context, for example, lipid-bilayer membranesare liquid in nature, but the critical-micelle concentra-tion is very low, so that lipids essentially never detachfrom the membrane. In the modelling and simulationof lipid membranes, a similar problem of attractive in-teractions and fluidity exists as for cell monolayers – toostrong attraction leads to solidification. In coarse-grainedsimulations of lipid bilayers, this problem was overcomeby an interaction potential with an extended range com-pared with the standard Lennard-Jones potential, whichprovides strong adhesion while still giving enough wiggleroom for the molecules that the membrane to remainsthe fluid phase [31]. In the spirit of minimal modeling,we propose an active Brownian particle (ABP) modelfor the cells, combined with a similar longer-range inter-action potential as employed successfully for the mem-brane lipids. We demonstrate that the LJ potential witha wider attractive basin indeed opens up a region in phasespace that displays liquid-vacuum coexistence. The sizeof the liquid-vacuum region expands as the basin widthof LJ potential increases. The fluidity of the condensedphase implies the emergence of several interesting be-haviors, like a tensile stress within the colony due to apreferred orientation of the boundary cells to the out- a r X i v : . [ phy s i c s . b i o - ph ] J un side, as awell as the detachment of cell cluster above asize threshold. When a coupling of cell polarity to theinstantaneous direction of motion – which is significantlyaffected by the interactions with the neighboring cells –is introduced, the formation of swirls and collective celldetachment is strongly enhanced. II. RESULTS AND DISCUSSIONA. Active Brownian particles with attraction
Liquid-vacuum coexistence requires strong inter-particle adhesion, so that cells can not detach from themain colony. Concurrently, the adhesion has to provideenough wiggling room that the cells remain locally mo-bile inside the condensed phase and provide fluidity tothe colony. A long-range coordinated motion of cells,like fingering or swirls, then already emerges to some ex-tent from the self-organized motion of cells which all varytheir propulsion direction independently and diffusively.However, pronounced correlations are found to requireadditional alignment interactions of cell orientation anddirection of motions. Here, the effect of neighbors push-ing (or pulling) a cell in a certain direction is assumedto induce a reaction in the cell to reorient and align itspropulsion direction with its instantaneous velocity di-rection.The ABP model, where each particle is a sphere (in3D) or a disc (in 2D) which undergoes rotational Brow-nian motion and additionally experiences a body-fixeddriving force of constant magnitude, was developed todescribe active motion on the microscale [15, 32]. Thismodel displays a rich phase behavior, most notablymotility-induced phase separation [24–26, 33], where per-sistence of motion and short-range repulsion induce clus-ter formation. Addition of a Lennard-Jones attractionleads to the formation of arrested clusters for smallpropulsion [28, 29, 34]. In order to obtain liquid-likeproperties at strong adhesion, we follow the spirit ofRef. [31] and propose an interaction potential with anextended basin of width ¯ σ , so that V m = (cid:15) (cid:2) ( σ/r ) − ( σ/r ) (cid:3) , < r ≤ / σ − (cid:15), / σ < r ≤ ˜ r (cid:15) (cid:2) ( σ/ ( r − ¯ σ )) − ( σ/ ( r − ¯ σ )) (cid:3) , ˜ r < r ≤ r cut (see Appendix A Fig. 6). Here, σ is the particle diameter,˜ r = (2 / σ + ¯ σ ), and (cid:15) is the interaction strength. Thismodified interaction provides a short-range repulsion orvolume exclusion for the particles with separation r < / σ , a force-free regime for 2 / σ < r ≤ ˜ r , and a long-range attraction for particle separation, ˜ r < r < r cut =2 . σ .For activity, each ABP is subjected to a constantactive force f along a body-fixed propulsion direction ˆn i = (cos θ i , sin θ i ). The orientation θ undergoes diffusivereorientation, and may additionally experience alignment forces. Time evolution follows a Langevin dynamics, m ¨ r i = − γ ˙ r i + F i ( r i ) + f ˆn i + √ Dη Ti , ˙ θ i = (cid:112) D r η Ri . (1)Here, F i = −∇ i V describes the interaction with othercells with the total potential V as a sum of all pair inter-actions, and f = v γ is the driving force which resultsin a self-propulsion velocity v for an isolated cell expe-riencing a drag force due to substrate friction with dragcoefficient γ , which is related to the thermal translationaldiffusion coefficient D = k B T /γ by the Einstein relation.Similarly, D r is the rotational diffusion coefficient. Thenoise forces η are assumed to be Gaussian white-noisevariables with (cid:104) η i ( t ) = 0 (cid:105) and (cid:104) η i ( t ) η j ( t (cid:48) ) (cid:105) = δ ij δ ( t − t (cid:48) ).However, note that this is an active system, and thusboth diffusion processes can in principle be independentactive processes with different amplitudes, and thus donot need to satisfy the Einstein relation or fluctuation-dissipation theorem. In order to emphasize the impor-tance of rotational over translational diffusion, we choose D r σ /D = 3.The cohesive nature of modeled cell colonies dependson the competition between adhesion and self-propulsion.A key parameter is the potential width ¯ σ , which controlsthe fluid-like consistency of the colony. Further detailsabout model parameters can be found in the supportinginformation. In the simulation results described below,all quantities are reported in dimensionless units based onthermal energy k B T , particle diameter σ , and rotationaldiffusion time τ r = 1 /D r . We characterize the system bythree dimensionless numbers, the P´eclet number, P e = v σ /D = 3 v τ r /σ which quantifies the activity of thesystem, the adhesion strength U = (cid:15)/ ( k B T ), and thepotential width ¯ σ/σ which determines the wiggle roomof the cells.In order to introduce local velocity-orientation align-ment, we assume alignment between propulsion directionand velocity for each cell individually [35–37]. In our sim-ple stochastic model, Eq. 1, the orientation dynamics inthis case is determined by˙ θ i = (cid:112) D r η Ri − k e D r ∂∂θ ( n i . v i ) (2)Here, k e is the strength of particle alignment. The align-ment force can be interpreted as arising from a pseudopotential V a = − ( k e / v · n ), acting only on the ori-entation n , but not on v . Unless noted otherwise, re-sults concern systems withouth velocity-alignment inter-actions (i.e. k e = 0). B. Liquid-Vacuum Coexistence
We begin our analysis by exploring the available phasespace spanned by activity
P e and adhesive strength U .The system is initialized with a circular cell colony with N = 7851 particles and a diameter 100 σ (packing frac-tion φ = 0 .
79) in a square simulation box of linear size150 σ . The resulting phase behavior as a function of P e and U is displayed in Fig. 1. Here, snapshots of parti-cle conformation after long simulation time ( t = 3300 τ r ),together with particle mobility, measured by the meansquared displacement d m = ( r i ( t + t (cid:48) ) − r i ( t )) (3)averaged over several reorientation times t (cid:48) = 12 τ r , areemployed to characterize the phases. Figure 1 shows thata line P e (cid:39) U separates a homogeneous gas phase at P e > U from a two-phase coexistence regime for
P e < U .For low activity (
P e (cid:28) U ), the condensed phase is solid,where particles do not show any significant movement,i.e. d m ≈
0. As activity increases and approaches
P e (cid:39) U , cells become mobile ( d m > P e (cid:38) (cid:29) U , attraction becomes negligible and con-ventional motility-induced phase separation is observed.A simple calculation, which equates the propulsion forcewith the maximum of the attraction force, reveals thatthe detachment of particle pairs occurs at P e = 2 . U ; forlarger P e the adhesive force is no longer strong enough tokeep particles together. Note that thermal fluctuationsare usually rather small in this study (because U (cid:29) P e (cid:46) U , cells are unable to detach from the colony,and the colony coexists with a cell vacuum outside. If P e (cid:28) U and U (cid:38)
8, the system is clearly kinetically ar-rested, but as activity increases, the “wiggle room” of thepotential allows particles to break the neighbor cage andmove, resulting in liquid-vacuum (L-V) coexistence. Thisstate of a single cohesive colony is not induced by the ini-tial conditions of a single circular patch, but also emergesfrom an initial random distribution of particles due toparticle aggregation and cluster coarsening. To quanti-tatively characterize and clearly distinguish mobile co-hesive colonies from the kinetically-arrested colonies, weemploy the ”mean squared particle separation” (MSPS).We choose random pairs of cells m and n inside the colonyat time t p , which are initially at contact with a center-center distance 1 . σ , and measure the squared separationof this pair over time. n average over N p such pairs atdifferent initial times t p yields M SP S ( t ) = 1 N p (cid:88) N p ( r m ( t p + t ) − r n ( t p + t )) . (4)A characteristic feature of the arrested dynamics in asolid phase is that particles do not exchange neighbors,so that the MSPS plateaus at M SP S < (1 . σ ) . Ina fluid phase, particles exchange neighbors at a constantrate, and MSPS increases linearly with time (see also Ap-pendix B, Fig. 8). Thus the MSPS is good indicator offluid-like behavior. Here, we choose M SP S > (1 . σ ) at time t = 12 τ r as a definition of fluid-like behav-ior. To quantify cohesiveness, we turn to a cluster anal-ysis, where particles are identified to be in the samecluster if their distance is less than the cutoff distance r cut . The condensed phase-vacuum coexistence is thensignaled by cluster number N cl = 1. Figure 2(A) dis-plays M SP S ( t = 12 τ r ) and N cl as a function of P e/U .For
P e (cid:46) . U , the system remains cohesive and solid.As activity increases, MSPS increases as well, but thecolony remains cohesive, clearly identifying the liquid-vacuum (L-V) coexistence region. Further increasing ac-tivity ( P e (cid:38) . U ) leads to the occasional detachmentof small clusters (above a threshold size) from the parentcolony (see discussion below). Interestingly, occasionalcluster detachment is not sufficient to disintegrate theparent colony, as detached cluster can rejoin the parentcolony, which thereby coexists with a gas of small clus-ters.Figure 2(B,C) display different cuts through the phasespace, to elucidate the region of stability of differentregimes. The results in Fig. 2(B) show that a minimumwidth ¯ σ/σ (cid:39) . σ/σ plays a crucial role to achieve a cell colony withfluid-like dynamics at strong adhesion. The importanceof P e/U as the relevant variable to distinguish two-phasecoexistence from a one-phase gas-like region, is empha-sized by Fig. 2(C), which demonstrates that the bound-aries between the different regimes occur at
P e/U (cid:39) . U (cid:38)
20. Note that all these bound-aries appear at
P e/U values, which are much smallerthan the unbinding threshold
P e/U (cid:39) . P e (cid:38) . U , small clusters areable to detach from the parent colony. This process canbe characterized by the cluster-size distribution P ( n pc ),see Fig. 2(D). The peak of the distribution for clustersin the size range from 10 to 100 indicates that particlesescape collectively. We do not observe the escape of anysingle cell from the colony in this regime. This can beunderstood from a simple argument, which considers asmall semi-circular patch of n pc particles at the bound-ary of the colony (see Appendix B, Fig. 9). The patchhas an interface with the colony of length proportionalto √ n pc . If all particle orientations point in roughly thesame direction (outwards), then the patch can unbindwhen n pc > n ∗ pc (cid:39) . U/P e ) , i.e. for sufficiently largepatch size, a size which decreases rapidly with increasing P e (see Appendix B for details, in particular Fig. 10).The probability for all particles in such a cluster to beroughly aligned depends on the P´eclet number, as parti-cles move toward the boundary with preferred outwardorientation [38]. However, the particle mobility in thefluid phase is very small due to the dense packing ofneighbors, so that the characteristic ballistic motion ofABPs for times less than τ r is completely suppressed (seeAppendix B, Fig. 7). Therefore, polar ordering is mainlyseen at the edge of the colony, see Fig. 3(A). Clusterformation therefore arises mainly from the increased mo-bility of pre-aligned particles at the boundary. FIG. 1: Phase diagram illustrated by snapshots at the end of simulation. Here we plot the mobility profile as a functionof attractive interaction U and activity P e . We start the simulation with initial circular patch. The colour code defines themagnitude of mobility. Blue means immobile and red means highly mobile colony. (¯ σ/σ = 0 .
3; overall packing fraction is0 . U = 40, P e = 30.
C. Stress Profile - Tensile Colonies
To gain a better understanding of the properties of theliquid-like cohesive colony, we analyze the polarizationfield of the active force and the stress profile. Figure 3(A)shows that the averaged local polarization is zero insidethe colony, but points outwards at the boundary. This isin contrast to what is typically found for motility-inducedclustering [24, 26, 39–41]. The reason is that the attrac-tive interactions keep outward-oriented particles at theboundary of the colony – which would otherwise moveaway – combined with the fluidity of the colony whichallows local particle sorting near the boundary, similarto the behavior of isolated self-propelled particles in con-finement [38, 42].The alignment of motility forces at the boundaryshould lead to an increase of tensile stress. For the liquidcolonies in coexistence with the vacuum phase, we findsignificant tensile stress in the center (see Fig. 3(C)). Asexpected from force balance, the stress is nearly constantinside the colony, but rapidly decreases in the bound-ary region where the tension is generated. The totalstress in the colony has three contributions: the inter-particle force, kinetic contribution, and swim stress. Atliquid-vacuum coexistence, the inter-particle contribu-tions plays the dominant role, whereas the swim stressis comparatively small (Fig. 3(B)). Figure 4 shows thedependence of the total central stress in the colony onthe activity
P e . In the solid regime, the central stressis negative due to the passive surface tension, resultingin a Laplace pressure proportional to U . Increasing ac-tivity leads to more liquid-like consistency, facilitating enhanced outward particle orientation at the edge, andhence tensile stress. Interestingly, the stress is found tobe a linearly increasing function of P e over the wholeinvestigated range, 0 < P e ≤
30, i.e. both in the solidand liquid regime of the colony. This indicates that theenhanced particle sorting occurs mainly near the edge,and an increased edge mobility exists already in the solidphase near the S-L phase boundary. A tensile stress incell colonies is observed similarly in experiments, wherethe average stress within a spreading cell sheet increasesas a function of distance from the leading edge [22]. Ina quasi-one-dimensional (rectangular channel) geometry,the total stress in the colony is obtained by integrationof the net active forces, and increases from zero outsideto a finite tensile stress in the center (see Appendix B,Fig. 11).
D. Coordinated Motion – Motion Alignment
In experimental observations, long-range velocity cor-relations are often visible in the bulk of spreading epithe-lial sheets [11, 12]. ABPs with adhesion display signifi-cant velocity correlations even without explicit alignmentinteractions (see Appendix Fig. 13, and Refs. [15, 43]).However, as ABPs display independent orientational dif-fusion, it is evident that realistic long-range correlationsrequire some type of velocity alignment. We employ alocal velocity-orientation alignment mechanism, in whichcell propulsion direction (=cell polarity) relaxes towardthe instantaneous cell velocity, resulting from the forcesinduced by its neighbors [35–37], as introduced in Eq. 2. (A) M SPS ( τ r ) Pe/U
MSPS S o li d - V ac uu m ( S - V ) L i qu i d - V ac uu m ( L - V ) L - V + D e t ac h e m e n t M a ny s m a ll c l u s t e r s < N c l > Pe/U
Detachment Many small clusters (B) - σ / σ Pe/U (C) S o li d - V ac uu m ( S - V ) L i qu i d - V ac uu m ( L - V ) L - V + d e t ac h e m e n t M a ny s m a ll c l u s t e r s U Pe/U (D) P ( n p c ) n pc ( LV + cluster detachment) -3 -2 -1 P ( n p c ) n pc FIG. 2: (A) ”MSPS” (left axis) and N cl (right axis) as a function of P e /U , for fixed U = 40 with increasing activity. changing P e value.
MSP S ( t ) is calculated at time separation t = 12 τ r . The orange area indicates S − V coexistence, the blue area L − V coexistence, the green area L − G coexistence with small detached cell clusters in the gas phase, and the red area ahomogeneous phase of many small clusters. Results are for ¯ σ/σ = 0 .
3. (B) Phase diagram as a function of the potential width¯ σ for fixed U = 40. (C) Phase diagram various U as a function of rescaled P´eclet number, P e/U , for fixed ¯ σ/σ = 0 .
3. (D)Cluster-size distribution of detached clusters, of size n pc , at ¯ σ/σ = 0 . U = 40, and P e = 32. Inset: Same data in log-logrepresentations, which also includes the parent colony along with the detached clusters.
0 10 20 30 40 50 60 70 80 90 (A) < p > r/ σ U=15, Pe=10U=40, Pe=30 -1 0 1 2 3 4 5 10 20 30 40 50 60 70 80 90 (B) S t r e ss r/ σ VirialSwimKEVirial+Swim+KE (C) S t r e ss r/ σ U=15, Pe=10U=40, Pe=30
FIG. 3: (A) Averaged polarization vector of the circular patch at ¯ σ = 0 . σ at L − V coexisting state for two different activityand (B) Different components of the stress at ¯ σ = 0 . σ and U = 40, P e = 30. ”Violet” represents the ”virial” stress profile,yellow color represents the stress profile due to activity, ”sky-blue” represents the kinetic contribution of the stress profile and”red” represents the sum of all these three contributions. (C) The total stress profile for different adhesive strength in L − V coexisting states at U = 15 ,
40. The stress starts generated at the boundary region of the colony and outside of it, the stressvanishes.
Without alignment, correlations arise from a small group of cells pointing in the same direction by chance, and thus -4-2 0 2 0 5 10 15 20 25 30 S t r e ss Pe U=40, - σ =0.3 σ FIG. 4: Total central stress calculated in the area startingfrom the center to a radius 30 σ as a function of activity, P e ,for U = 40 and ¯ σ = 0 . σ . The transition from S-V to L-Vcoexistence occurs at P e = 23 (compare Fig. 2(C)), whichessentially coincides with P´eclet number where central stresschanges sign. moving collectively more easily and furthermore draggingother cells along. The alignment interaction stabilizesand enhances this effect. In presence of velocity align-ment, the velocity field shows an enhanced coordinatedmotion with prominent swirls in the bulk of the colonyand fingering at the edge (see Fig. 5(A)).We quantify the spatial correlations by the velocity-velocity correlation function C vv ( r ) = (cid:42) (cid:80) r i δ v ( r i ) · δ v ( r i + r ) (cid:80) r i δ v ( r i ) · δ v ( r i ) (cid:43) , (5)as a function of distance r , where the brackets denotean average over all directions and time. Here, veloci-ties are measured relative to the average velocity ¯ v , ofthe whole colony, i.e. δ v ( r ) = v ( r ) − ¯ v , to avoid finite-size effects. The correlations decay exponentially with acharacteristic length scale, ξ vv (see Appendix B, Fig. 14).Figure 5(B) displays the correlation length ξ vv as a func-tion of alignment strength k e for various adhesive in-teractions. For fixed adhesion and activity, increasingalignment strength k e facilitates a transition from thesolid to the liquid state of the colony. Furthermore, thealignment coupling leads to stronger correlations, as in-dicated by the monotonic increase of ξ vv with k e v , andthus to swirls and fingers. Eventually, fingering is sostrong that clusters detach, and the colony is no longercohesive. However, correlation lengths up to ξ vv = 10 σ can be achieved, quite comparable to the 5 to 10 timescell size obtainable in experiments [11, 15]. Also, thetensile stress at the colony center increases (see inset ofFig. 5(B)) and becomes positive at sufficiently large k e v .The critical alignment strength k e v , where the colony isliquefied and the tensile stress becomes positive, increaseswith attraction strength U . σ (A) (B) ξ vv / σ k e v U=15, Pe=8U=25, Pe=15U=40, Pe=25 -1 0 1
0 0.5 1 1.5 2 C e n t r a l S t r e ss k e v U=40, Pe=20U=25, Pe=13
FIG. 5: Dynamics of cell colony with orientation-velocityalignment interactions. (A) Velocity fluctuation field of thefull cell colony for U = 40, P e = 25 and k e v = 1 . σ = 0 . σ , which displays prominent swirls in the bulkand finger-like structures at the edge. See also Appendix B,Fig. 13. (B) Characteristic correlation length ξ vv , extractedfrom velocity correlation, as a function of alignment strength k e , for different attractive interactions U and activities P e , asindicated, with ¯ σ = 0 .
3. The colony is in the fluid state for k e v = 0. Inset: Variation of total central stress as a functionof alignment strength k e , for adhesive strength U = 25 , σ = 0 . k e = 0, and transits to a liquid state at k e v (cid:39) .
65 and k e v (cid:39) . U = 25 and U = 40, respectively. III. CONCLUSIONS AND OUTLOOK
We have presented a minimal model for the fluidiza-tion of cell colonies, which consists of active Brownianparticles with adhesion. An attractive potential with in-creased basin width yields non-equilibrium structures,phase behavior and dynamics, which capture relevantfeatures of biological cell colonies. The main observa-tion is that for moderate adhesion and propulsion, thesystem exhibits liquid-vacuum coexistence, i.e. all par-ticles agglomerate into a single colony displaying liquid-like properties, while the outside remains devoid of anyparticles. This is reminiscent of in vitro experiments ofMDCK colonies, where cells show strong motion, whileremaining perfectly cohesive. Furthermore, the fluidityof the colony in our model results in outward ordering ofparticle orientations at the edge, thus leading to tensionin the colony. This is consistent with the results of trac-tion force microscopy, which show that MDCK coloniesare typically under tension [22]. Our model demonstratesthat no alignment interaction or growth mechanism needto be evoked to explain such tensile forces – the motilityof the cells combined with liquid properties of the colonysuffice. As motility force increases, particles start to de-tach from the parent colony, however not as single cellsbut collectively in small clusters of cells. Finally, we havedemonstrated how velocity-polarity alignment can fur-ther enhance fluidity, tension, and fingering of the colony,and collective cell detachment. Indeed, with velocity-polarity alignment, simulations look very reminiscent ofreal MDCK colonies, displaying strong fingering at theedge, high tension and long-ranged velocity correlations.Our model also provides a tentative explanation foranother biological phenomenon. When metastatic cellsdetach from a tumor, they typically detach collectively,as small groups of five cells or large aggregates [44–46],into the stroma and migrate to reach blood or lymphvessels. At the edge of the liquid-vapor region of ourmodel, particles show exactly this type of behavior; thecolony is no longer perfectly cohesive, but clusters of cellsbegin to detach.An interesting next question is how these results willbe affected by cell growth. Of course, if growth is slow,the dynamics will be independent of growth and the phe-nomenology will be unchanged. However, when timescales of growth and motion become comparable, novelphenomena may arise.
Acknowledgments
Financial support by the Deutsche Forschungsgemein-schaft (DFG) through the priority program SPP1726“Microswimmers – from single particle motion to col-lective behavior” is gratefully acknowledged. The au-thors also gratefully acknowledge the computing timegranted through JARA-HPC on the supercomputer JU-RECA [47] at Forschungszentrum J¨ulich.
Appendix A: Attractive Brownian Particles: Modeland Analysis of Simulation Data1. Simulation parameters
For numerical implementation of our model, we usethe LAMMPS molecular simulation package, with in-house modifications to describe the angle potential andthe propulsion forces, as described in the main text. Thesystem consists of N = 7851 particles (cells) in a 2Dsquare simulation box of size L x = L y = 150 σ withperiodic boundary conditions, unless noted otherwise. With velocity-orientation-alignment interaction, the sim-ulation is carried out in a box of size L x = L y = 250 σ .For the extended LJ interaction potential, see Fig. 6, weuse the cut-off distance r cut = 2 . σ . For numerical effi-ciency, we chose a finite mass m = 1 and drag coefficient γ = 100 such that the velocity relaxation time m/γ ismuch smaller than all physical time scales. The equa-tions of motion are integrated with a velocity Verlet al-gorithm, with time step ∆ t = 0 . × time steps, with rotationaldiffusion coefficient D r = 0 .
03 this corresponds to a to-tal simulation time longer than 3000 τ r , where τ r is therotational decorrelation time.
2. Polarization Vector
We define the spatial-temporal average polarization p for the quasi-circular cell colony by the projection of theorientation vector ˆn of the particle on the radial directionfrom the center of mass of the colony, i.e., (cid:104) p ( r (cid:48) , t ) (cid:105) = N (cid:88) i =1 ( ˆn i · ˆr (cid:48) i ) δ ( r (cid:48) − | r (cid:48) i | ) / N (cid:88) i =1 δ ( r (cid:48) − | r (cid:48) i | ) (A1)where r (cid:48) i = r i − r cm , and r cm is the center-of-mass posi-tion at a particular time t . Here, δ ( r ) is a smeared-out δ -function of width σ . (cid:104) p (cid:105) is further averaged over time.
3. Stress calculation
In ABP systems with short-range repulsion, it has beenshown that the pressure is a state function, dependingonly on activity, particle density, and interaction po-tential, but not on the interaction with confining walls[41, 48–50]. In comparison to passive systems, activityimplies a new contribution to pressure, which is calledthe swim pressure. The calculation of the local stress inan ABP system is a matter of an ongoing debate, whichmainly concerns the form of the active term. We followRef. [41], and define the stress in a volume ∆ V by∆ V Σ αα = N (cid:88) i =1 m (cid:104) ˙ r i Λ i (cid:105) − γγ R N (cid:88) i =1 (cid:104) v n i · ˙ r i Λ i (cid:105) N (cid:88) i =1 N (cid:88) j =1 (cid:104) λ ij r ij · F ij (cid:105) (A2)where ˙ r i , F i ( i = 1 , ..., N ) denote the velocity and forceof particle i , respectively. F ij represents the pair wiseinteraction between particle i and j and r ij = r i − r j .Here Σ αα are the diagonal stress-tensor components. γ R is the damping factor which is related to the rotationaldiffusion coefficient as γ R = 2 D r . Λ i determines the vol-ume ∆ V , where Λ i ( r ) is unity when particle i is within∆ V and zero otherwise. λ ij denotes the fraction of the ¯ σ . . . . . . . − − . . r/σ V m / (cid:15) V mod V LJ FIG. 6: Modified Lennard-Jones (LJ) potential (red coloredcurve) is modified by inserting a plateau of width ¯ σ at theminimum of this potential. The usual LJ potential is shownby blue crosses. line connecting particle i and j inside of the volume ∆ V .The first and the last term in (A2) are the classic kineticcontribution and the contribution of inter-particle inter-actions. The second term denotes the active-force con-tribution in the stress calculation. Notably, in the fluid-vacuum state we are focusing on in this work, the activestress component is negligible compared to the inter par-ticle interaction contribution. This is in line with resultsfor the pressure contributions in repulsive ABP systemsat coexistence between a high-density and a low-densityphase, where the swim pressure in the high-density phaseis negligible [41]. Appendix B: Supporting Considerations and Results1. Particle Mobility in Vacuum, Gas, and FluidPhases
A single, isolated ABP has a characteristic meansquare displacement, with ballistic motion (MSD ∼ v t )at short times t < τ r , and a diffusive motion (MSD ∼ v τ r t ) for t > τ r [51]. This is the behavior we ob-serve in the vacuum phase, see Fig. 7. A very similarbehavior is found in the gas phase, at packing fraction φ = 0 . τ r remains unaffected. How-ever, the behavior changes drastically in the fluid-likephase, where the collisions and attractive interactionscompletely suppress the ballistic regime, see Fig. 7.A very similar dynamic behavior is observed in themean squared particle separation (MSPS), see Fig. 8.The time dependence in the fluid phase is dominatedby linear diffusion behavior, while in the solid (jammed)phase it is strongly sublinear. t t ∆ r / σ t/ τ r MSD in fluid colonyMSD in gas phaseFree-ABP
FIG. 7: Mean square displacement (MSD) in three differentstates: (i) In the bulk of ”fluid-like” colony (purple circles)with U = 40 and P e = 30, (ii) in the ”gas” phase (red squares) U = 40 and P e = 41 and packing fraction φ = 0 . P e = 41. t t M SPS t/ τ r Pe/U=0.57Pe/U=0.50
FIG. 8: Mean squared particle separation (
MSP S ) as a func-tion of time for two different
P e numbers at U = 40 and¯ σ = 0 .
3. At
P e/U = 0 . MSP S remains constant at shorttimes, indicative of no neighbor exchange, and shows sub-diffusive behavior at longer time. However, at
P e/U = 0 . MSP S shows a long-time linear behav-ior, indicative of constant neighbor exchange and liquid-likebehavior.
2. Minimum Cluster Size for Detachment
For our model system, particles can only escape fromthe parent colony in the form of a small cluster at largeradhesive strength ( U ≥ n pc be the number of particles in the cluster.For packing in a roughly triangular lattice with latticeconstant a = σ + ¯ σ/
2, this implies a cluster radius L FIG. 9: A semi-circular cluster of particles formed at thecolony edge. All particles in this region are assumed to bealigned and to be oriented in the outward direction, normalto the interface, as indicated by the arrows. The length of theinterface between parent colony and detaching cluster is L . R cl /a = ( √ /π ) / n / pc . The length of the interface be-tween cluster and parent colony is L = 2 R , see Fig. 9.Along the interface, there are L/a bonds between par-ticles on both sides of the interface, which generate thesame maximum force as for the detachment of a singleparticle, where bond breakage occurs at
P e/U = 2 . n pc P e = 2 . L/a ) U = 4 . √ /π ) / U n / pc ,which implies n pc,min = √ π (cid:18) . UP e (cid:19) . (B1)Hence, with increasing activity P e at constant adhe-sive strength U , the minimal size of detached clustersis expected to decrease rapidly. Figure 10 shows simula-tion results for the number of particles present in thesmallest cluster n pc,min as a function of P e/U . Thisdemonstrates that beyond the liquid-vacuum region, n pc rapidly decreases with increasing activity P e and even-tually reaches a ”gas-like” phase, where single-particledetachment from the parent colony, i.e. n pc,min = 1, isobserved.Furthermore, we can use Eq. (B1) to estimate the clus-ter size when cluster break-off first becomes possible, at P e/U (cid:39) .
75, which is about n pc (cid:39)
20, in reasonableagreement with the lower cutoff of the cluster-size distri-bution in Fig. 2(D) of the main text.
3. Probability Distributions of Aligned ParticleClusters at Colony Edge
For randomly oriented particles, the probability to finda cluster in which the orientations of all particles havea positive projection into one chosen direction is 2 − N cl ,which is very small for clusters of size 10 or larger. How-ever, there is a sorting mechanism which strongly en- < N c l > Pe/U
65 at the bound-ary. Thus, properly oriented particles only have to dif-fuse laterally along the boundary to form clusters for de-tachment. This mechanism is supported by simulations,which allow the tracking of the history of cluster devel-opment.
4. Stress Calculation in a Quasi-One-DimensionalGeometry
As a further test to our stress estimates, we use forcebalance to obtain an independent measure of stress, sim-ilar to traction-force microscopy setups [22]. Under thephysical interpretation of our system that the particlesare cells which exert an active force γv on the substratein order to move against friction forces − γv , the trac-tion force of each particle is T = γv ˆn − γ v . In onedimension, force balance is closed, and we can calcu-late a change of stress via force balance. We simulate aquasi-one-dimensional geometry with a nearly flat inter-face. The system consists of a 2D channel of dimensions L x = 6 ∗ L y and L y = 20 σ , filled with N = 1200 particlesarranged initially to fill half the channel. This systemis subjected to periodic boundary conditions in both x and y directions. The stress within the cell colony is cal-culated via integration of force balance (assuming zerostress outside the colony).With the parameters ¯ σ = 0 . σ , adhesive strength U = 15 and activity P e = 13, the cell colony is in liquid-vacuum coexistence, see Fig. 11(top). The stress pro-0file in the direction ( x ) normal to the interface is shownin Fig. 11(bottom), while the tangential stress vanishes.Figure 11 also shows that the estimations of the stressprofile calculated from the traction forces and from thevirial expression agree quite well. FIG. 11: Top: Snapshot from the simulation in a quasi-one-dimensional rectangular channel, for a cell colony in liquid-vacuum coexistence. The parameters are ¯ σ = 0 . σ , adhesivestrength U = 15 and activity P e = 13. Bottom: Stress cal-culated from traction forces and virial contribution.
5. Origin of Tension in Attractive ABP Clusters –a Toy Model
Why are colonies of attractive ABPs under tension,while cluster of repulsive ABPs in the state of motility-induced phase separation are under pressure? A sim-ple toy model can elucidate the underlying mechanism.Consider two ABPs which are connected by a harmonicbond [52]. This bond represents the interaction betweenan ABP at the colony edge, and one (or more) ABPfurther inside the bulk. In case of a sufficiently highP´eclet number (large propulsion, slow rotational diffu-sion), the dumbbell quickly reaches a quasi-stationary,torque-free state, see Fig. 12(A,B). In this state, the forcecan be separated into a propelling component with direc-tion ( n + n ), normal to the instantaneous bond vector r ∼ ( n − n ), and a bond stretching component, seeFig. 12(B,C). The stretching force is f ext = f | n · ˆ r | = f | cos( θ ) | = f (cid:112) (1 − n · n ) / | sin(( ϕ − ϕ ) / | (B2)where ϕ i is the orientation angle of n i with respect tosome fixed axis. This stretching force has to be averagedover all orientations n and n , which yields (cid:104) f ext (cid:105) = 2 π f (B3) (A)(B)(C) FIG. 12: Arrangement and forces of ABP dumbbells. (A)Two connected ABPs with randomly oriented propulsionforces f and f . (B) A short time t later, with σ/v < t < τ r ,the ABP orientations remain essentially unchanged, but theparticles have rearranged to form a quasi-stationary, torque-free state, in which the bond is under tension. The dumbbellalso moves in direction n = n + n normal to the bond vector r . (C) In the quasi-stationary state, the ABP orientation vec-tors form angles θ and θ with the bond vector r ∼ ( n − n ),with θ = θ . Stress Σ αα is force per length, i.e. Σ αα (cid:39) (cid:104) f ext (cid:105) /a (cid:39)(cid:104) f ext (cid:105) /σ . We can thus use Eq. B3 to estimate the tensilestress as a function of P´eclet number. With f = γv , D = k B T /γ , D r = 3 D/σ , and P e = 3 v / ( σD r ), weobtain f = P e (in our dimensionless units). Thus, wepredict Σ = Σ + (2 /π ) P e from our toy-model calcula-tion. The linear dependence agrees well with the simu-lation results for the cell colony, see Fig. 4 of main text.However, the slope estimated from Fig. 4 is 0 .
16, abouta factor 4 smaller than the toy-model prediction. Twoobvious reasons for this overestimation of the slope inthe toy model are that (i) the bond vector takes all ori-entations with equal probability, but only orientationsroughly perpendicular to the interface contribute to thestress (factor 2), and (ii) the hard-core repulsion betweenABPs is neglected, which sometimes leads to a pressure(negative tension) (maybe another factor 2) – so that theoverall agreement is quite satisfactory.
6. Velocity Correlation Function
To quantify collective cell migration, we map out thevelocity field. Snapshots of the simulations in Fig. 13demonstrate that particles display strongly coordinatedmotion. To further quantify the correlations, we calcu-late the velocity correlation function C vv as described inEq. 5 of the main text. Figure 14 shows C vv ( r ). Onshort length scales, the velocity correlations decay expo-nentially with a characteristic length scale ξ vv . We esti-mate ξ vv by fitting the simulation data by exp( − r/ξ vv ).The dependence of the ξ vv on U and P e is discussed inthe main text, see Fig. 5(B).1 σ (A) σ (A) (B) σ (B) (C) FIG. 13: (Top) Velocity field (left) without and (right) withexplicit alignment interaction at U = 40, P e = 25 ¯ σ = 0 . k e v = 1 .
25. (Bottom) Velocity field in the finger-likestructure of the fluid-like in presence of alignment interaction( k e v = 1 . (A) C vv (r / σ ) r/ σ U=40, Pe=25, k e v =1.25, - σ =0.3 σ U=25, Pe=15, k e v =2.0, - σ =0.3 σ U=15, Pe=8, k e v =2.16, - σ =0.3 σ U=40, Pe=30, k e v =1.0, - σ =0.2 σ U=25, Pe=20, k e v =1.0, - σ =0.2 σ (B) C vv (r / σ ) r/ σ U=40, Pe=25, k e v =1.25, - σ =0.3 σ U=25, Pe=15, k e v =2.0, - σ =0.3 σ U=15, Pe=8, k e v =2.16, - σ =0.3 σ U=40, Pe=30, k e v =1.0, - σ =0.2 σ U=25, Pe=20, k e v =1.0, - σ =0.2 σ (C) C vv (r / σ ) r/ σ WO-align (k e v =0.0), U=15, Pe=8W-align (k e v =2.16), U=15, Pe=8 FIG. 14: Velocity correlation function C V V ( r/σ ) for a systemwith velocity-alignment interaction in the L-V region. (A)Spatial dependence for various attraction strengths U , P´ecletnumbers P e , and alignment interaction strengths k e v , andwidths ¯ σ of the attraction well, as indicated. (B) Same dataas in (A), with exponential decay demonstrated by log-linrepresentation. (C) Comparison of correlation functions with(w-align) and without velocity alignment (wo-align). Figure 14 shows C vv with and without alignment inter-action. Velocity-alignment interactions strengthen corre-lated motion and result in more swirl-like patterns, asindicated by a small negative minimum in C vv . [1] P. Martin and S. M. Parkhurst, Development , 3021(2004).[2] P. Friedl, Y. Hegerfeldt, and M. Tusch, Int. J. Dev. Biol. , 441 (2004).[3] V. Lecaudey and D. Gilmour, Curr. Opin. Cell. Biol. ,102 (2006).[4] K. Burridge and M. Chrzanowska-Wodnicka, Annu. Rev.Cell Dev. Biol. , 463 (1996).[5] G. Giannone, B. J. D. Thaler, O. Rossier, Y. Cai,O. Chaga, G. Jiang, W. Beaver, H. D¨obereiner, Y. Fre-und, G. Borisy, and M. P. Sheetz, Cell , 561 (2007). [6] D. A. Lauffenburger and A. F. Horwitz, Cell , 359(1996).[7] K. Keren, Z. Pincus, G. M. Allen, E. L. Barnhart,G. Marriott, A. Mogilner, and J. A. Theriot, Nature ,475 (2008).[8] Y. Matsubayashi, M. Ebisuya, S. Honjoh, and E. Nishida,Curr. Biol. , 731 (2004).[9] A. Puliafito, L. Hufnagel, P. Neveu, S. Streichan, A. Si-gal, D. K. Fygenson, and B. I. Shraiman, Proc. Natl.Acad. Sci. U.S.A. , 739 (2012).[10] M. Bindschadler and J. L. McGrath, J. Cell Sci. , 876 (2007).[11] T. E. Angelini, E. Hannezo, X. Trepat, M. Marquez, J. J.Fredberg, and D. A. Weitz, Phys. Rev. Lett. , 168104(2010).[12] T. E. Angelini, E. Hannezo, X. Trepat, M. Marquez, J. J.Fredberg, and D. A. Weitz, Proc. Natl. Acad. Sci. U.S.A. , 4714 (2011).[13] D. Bi, J. H. Lopez, J. M. Schwarz, and M. L. Manning,Nat. Phys. , 1074 (2015).[14] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning,Phys. Rev. X , 021011 (2016).[15] S. Garcia, E. Hannezo, J. Elgeti, J. F. Joanny, P. Sil-berzan, and N. S. Gov, Proc. Natl. Acad. Sci. U.S.A. , 15314 (2015).[16] M. Krajnc, S. Dasgupta, P. Ziherl, and J. Prost, Phys.Rev. E , 022409 (2018).[17] T. Omelchenko, J. M. Vasiliev, I. M. Gelfand, H. H.Feder, and E. M. Bonder, Proc. Natl. Acad. Sci. U.S.A. , 10788 (2003).[18] M. Poujade, E. G. Mongrain, A. Hertzog, J. Jouanneau,P. Chavrier, B. Ladoux, A. Buguin, and P. Silberzan,Proc. Natl. Acad. Sci. U.S.A. , 15988 (2007).[19] T. E. Angelini, E. Hannezo, X. Trepat, M. Marquez, J. J.Fredberg, and D. A. Weitz, Biophys. J. , 1790 (2010).[20] M. Vishwakarma, J. D. Russo, D. Probst, U. S. Schwarz,T. Das, and J. P. Spatz, Nat. Commun. , 3469 (2018).[21] We use the term ’vacuum’ here somewhat loosely to rep-resent a phase of extremely low cell density.[22] X. Trepat, M. R. Wasserman, T. E. Angelini, E. Millet,D. A. Weitz, J. P. Butler, and J. J. Fredberg, Nat. Phys. , 426 (2009).[23] M. Reffay, M. C. Parrini, O. Cochet-Escartin, B. Ladoux,A. Buguin, S. Coscoy, F. Amblard, J. Camonis, andP. Silberzan, Nat. Cell Biol. , 217 (2014).[24] Y. Fily and M. C. Marchetti, Phys. Rev. Lett. ,235702 (2012).[25] G. S. Redner, M. F. Hagan, and A. Baskaran, Phys. Rev.Lett. , 055701 (2013).[26] A. Wysocki, R. G. Winkler, and G. Gompper, EPL ,48004 (2014).[27] J. Elgeti, R. G. Winkler, and G. Gompper, Rep. Prog.Phys. , 056601 (2015).[28] G. S. Redner, M. F. Hagan, and A. Baskaran, Phys. Rev.E , 012305 (2013).[29] R. M. Navarro and S. M. Fielding, Soft Matter , 7525(2015).[30] V. Prymidis, S. Paliwal, M. Dijkstra, and L. Filion, J.Chem. Phys. , 124904 (2016).[31] I. R. Cooke, K. Kremer, and M. Deserno, Phys. Rev. E , 011506 (2005).[32] B. Smeets, R. Alert, J. Peˇsek, I. Pagonabarraga, H. Ra-mon, and R. Vincent, Proc. Natl. Acad. Sci. U.S.A. ,14621 (2016).[33] I. Buttinoni, J. Bialk´e, F. K¨ummel, H. L¨owen,C. Bechinger, and T. Speck, Phys. Rev. Lett. , 238301(2013).[34] B. M. Mognetti, A. Sari´c, S. Angioletti-Uberti, A. Cacci-uto, C. Valeriani, and D. Frenkel, Phys. Rev. Lett. ,245702 (2013).[35] B. Szab´o, G. J. Sz¨oll¨osi, B. G¨onci, Z. Jur´anyi,D. Selmeczi, and T. Vicsek, Phys. Rev. E , 061908(2006).[36] M. Basan, J. Elgeti, E. Hannezo, W. J. Rappel, andH. Levine, Proc. Natl. Acad. Sci. U.S.A. , 2452(2013).[37] K. N. T. Lam, M. Schindler, and O. Dauchot, New J.Phys. , 113056 (2015).[38] J. Elgeti and G. Gompper, EPL , 48003 (2013).[39] Y. Fily, S. Henkes, and M. C. Marchetti, Soft Matter ,2132 (2014).[40] P. Digregorio, D. Levis, A. Suma, L. F. Cugliandolo,G. Gonnella, and I. Pagonabarraga, Phys. Rev. Lett. , 098003 (2018).[41] S. Das, G. Gompper, and R. G. Winkler, Sci. Rep. ,6608 (2019).[42] Y. Fily, A. Baskaran, and M. F. Hagan, Soft Matter ,5609 (2014).[43] N. Sep´ulveda, L. Petitjean, O. Cochet, E. Grasland-Mongrain, P. Silberzan, and V. Hakim, PLOS Comput.Biol. , e1002944 (2013).[44] P. Friedl and D. Gilmour, Nat. Rev. Mol. Cell Biol. ,445 (2009).[45] J. J. Christiansen and A. K. Rajasekaran, Cancer Res. , 8319 (2006).[46] A. G. Clark and D. M. Vignjevic, Curr. Opin. Cell Biol. , 13 (2015).[47] J¨ulich Supercomputing Centre, J. Large-Scale Res. Facil. , A132 (2018).[48] S. C. Takatori, W. Yan, and J. F. Brady, Phys. Rev. Lett. , 028103 (2014).[49] A. P. Solon, Y. Fily, A. Baskaran, M. E. Cates, Y. Kafri,M. Kardar, and J. Tailleur, Nat. Phys. , 673 (2015).[50] T. Speck and R. L. Jack, Phys. Rev. E , 062605 (2016).[51] J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough,R. Vafabakhsh, and R. Golestanian, Phys. Rev. Lett. ,048102 (2007).[52] R. G. Winkler, Soft Matter12