Effect of contact inhibition locomotion on confined cellular organization

Experiments performed using micro-patterned one dimensional collision assays have allowed a precise quantitative analysis of the collective manifestation of contact inhibition locomotion (CIL) wherein, individual migrating cells reorient their direction of motion when they come in contact with other cells. Inspired by these experiments, we present a discrete, minimal 1D Active spin model that mimics the CIL interaction between cells in one dimensional channels. We analyze the emergent collective behaviour of migrating cells in such confined geometries, as well as the sensitivity of the emergent patterns to driving forces that couple to cell motion. In the absence of vacancies, akin to dense cell packing, the translation dynamics is arrested and the model reduces to an equilibrium spin model which can be solved exactly. In the presence of vacancies, the interplay of activity-driven translation, cell polarity switching, and CIL results in an exponential steady cluster size distribution. We define a dimensionless Péclet number Q—the ratio of the translation rate and directional switching rate of particles in the absence of CIL. While the average cluster size increases monotonically as a function of Q, it exhibits a non-monotonic dependence on CIL strength, when the Q is sufficiently high. In the high Q limit, an analytical form of average cluster size can be obtained approximately by effectively mapping the system to an equivalent equilibrium process involving clusters of different sizes wherein the cluster size distribution is obtained by minimizing an effective Helmholtz free energy for the system. The resultant prediction of exponential dependence on CIL strength of the average cluster size and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q^{1/2}$$\end{document}Q1/2 dependence of the average cluster size is borne out to reasonable accuracy as long as the CIL strength is not very large. The consequent prediction of a single scaling function of Q, particle density and CIL interaction strength, characterizing the distribution function of the cluster sizes and resultant data collapse is observed for a range of parameters.


Introduction
A distinctive feature of cell colonies is the emergence of collective organization, sensitive to the self propulsion characteristics of the individual constituent cells [1][2][3][4] .In general, the functionality of the tissues is determined by the ability of cells to self-organize into a diverse spectrum of spatio-temporal arrangement and patterns 1,2,5,6 .Collective cell migration is driven by the interplay of self-propulsion characteristics and the cell-cell interactions 1 .Drastic transformation of the underlying arrangement of the cells is associated with cell morphogenesis and cancer progression 1 .
For many cell types it has been observed that when individual cells come in contact with other cell, their direction of self propulsion tends to get re-polarized away from the neighbouring cells 7,8 .This process is also referred to as Contact Inhibition Locomotion (CIL) 1,2,7,8 .The CIL interaction arises due to the fact that when two cells collide, the cell front adheres to the colliding cell, which obstructs further cell movement 7,8 .This leads to repolarization of the cell's cytoskeleton.This in turn creates a new front away from the adhesion zone and facilitates the two cells to separate [7][8][9] .This inherent ability of the individual cell to reorient its direction of movement upon contact with other cells is an important determinant for multitude of cellular processes ranging from morphogenesis, tissue organization, wound healing and cancer progression among others [1][2][3] .The interplay of the distinct interactions at play at the level of individual cells including CIL provides the basis of the controlling mechanism which allows colonies of cells to organize themselves into diverse array of morphologies and functional characteristics 1 .
Phenomenological models for the collective behaviour of cell colonies and tissue reorganization have built on active matter approaches [10][11][12][13][14][15][16] .In this context, agent-based models have proven to be particularly useful for describing the collective behaviour of broad class of actively driven systems ranging from bacterial suspension to cell colonies and helped in providing crucial insights into the specific phenomenology of active matter 1,2,5,13,17,18 .On the other hand, active hydrodynamic models have served to highlight the unifying principles of active organization based on symmetry principles and conservation laws [19][20][21][22] .These approaches to understand actively driven systems have been complemented by active vertex models 23,24 and phase models 25 for some class of active systems.In particular, agent models which have explicitly incorporated the CIL interaction within their ambit, have been able to qualitatively explain tissue phenotypes 1,2 and theoretically predict transitions between different phases of collective organization in the cell colonies 1 .While many of these studies have had reasonable success in providing a qualitative understanding of several types of cell transitions, exploration of full phase parametric space for different cell types has remained an important open problem, owing to increased complexity arising from more detailed specification of cell-cell interaction rules.In this context, it may be noted that, while diverse agent based continuum models have been studied extensively, discrete driven lattice gas modeling approaches 16,[26][27][28] have remained relatively unexplored in the context of cell colonies and tissues.It is worthwhile to point out the potential of such minimalist modeling approach.In particular, variants of discrete driven lattice gas models have been used to describe transport across bio-membranes 29 , bidirectional cellular cargo transport 30,31 , collective transport of motor proteins on biofilaments 32 , and growth process of fungal mycelium [33][34][35] .
We also adopt a similar modeling approach to gain insight on the role of CIL in determining the organization of cells in quasi 1D settings, motivated by experiments performed with cells on 1D collision assays that have been designed using micropatterning techniques 7,8,21 .The restricted nature of cellular movement in such assays allows for more precise identification of collision event of cells apart from enhancing the efficiency of CIL response since the collisions between the cells is head on 7,8 .In particular, this experimental technique has been used to study and quantify the effect of CIL in cultured Xenopus Neural crest cells which are confined to move on micro-patterned fibronectin lines with very narrow width.The narrow width of these assays forces the cells to move along the narrow fibronectin line and to undergo a repolarization by 180 • due to CIL interaction 21 .
We put forward and study an Active spin model, belonging to the class of driven lattice gas models, which mimics the movement of individual cell and binary cell-cell interaction mediated by CIL in 1D assay.We use this model to investigate the dynamics of collective organization of cells and the clustering characteristics of cells that are subject to CIL interactions in such geometries.

Model and Methods
We consider a discrete 1D lattice consisting of L sites and represent individual cells as particles.The individual lattice sites can either be empty or be occupied by one particle.Each particle possesses discrete states of the individual polarization vector ⃗ P, associated with their direction of movement on the lattice.The polarization state of the particle at site i maybe described in terms of a variable σ i which can take two different values, ±1, depending on whether the polarization vector points towards right or left direction on the lattice, respectively.A general configuration of the system at a given time t maybe represented as, Here, (→) corresponds to a particle moving to the right, while(←) corresponds to a particle moving to the left, and 0 corresponds to an empty lattice site.
The primary characteristic of CIL is the propensity of the cells to align the direction of movement away from other neighbour cells.To mimic the effect of CIL in 1D channel, we consider a short range interaction between the particles described by a nearest neighbour interaction potential,
For this choice of H, the configurations of a pair of neighbouring particles on the lattice would have the following energy: where we have set k B T = 1.Thus particle configurations for which polarization vectors of neighbouring particles face each other are disfavoured, while configurations for which polarization vectors are oppositely aligned are favoured.We define a cluster as a continuous array of at least two particles that are bounded by vacancies at both ends.We consider that the switching dynamics of the particles between the different polarization states occurs only for clusters comprising of two or more particles, and is dictated by the interaction potential given by Eq. 1, such that the switching rates of particle's polarity, k s ∝ exp(−∆E), where ∆E is the energy difference between the two configurations.Thus the switching dynamics of the polarity state of a particle in the bulk of a cluster ( a particle bounded by other particles at its both end) maybe represented by, Here, ∆J = J 1 − J 2 and b is the switching rate of particle in a cluster between polarization state which have no energy difference.
Similarly, the switching dynamics of the polarity of the particle at the cluster boundary maybe represented by, As far as the translation dynamics is concerned, the particles at site i hops to the adjacent site in the direction of its polarization vector, with a rate a, provided that site is vacant, i.e., if σ i = +1 then the particle hops to site i + 1 with rate a, if it is vacant.On the other hand if σ i = −1, then it hops to site i − 1 with rate a, if it is vacant.These rules of particle movements maybe summarized as,

Connection and mapping to other models
Its worthwhile to point out that in the absence of CIL interaction within clusters, which corresponds to setting J 1 and J 2 to zero, the model is very similar to Persistent Exclusion Process (PEP) 13,31 with the crucial difference that while for our case a single particle ( which is not part of a cluster) does not undergo switching of their polarities, PEP allows for switching of a single particle in the lattice.

3/15
The Hamiltonian introduced in Eq.1 can also be decomposed into an Ising like alignment term and an interaction term which mimics the effective repulsion interaction arising due to reorientation of the particle polarities.Specifically, where, ni,i+1 is an unit vector along r i − r i+1 , corresponding to the position vector of the i th and (i + 1) th site in the lattice.This is indeed similar to the interaction potential invoked in hydrodynamic modeling of CIL phenomenology 21 .The second term of this interaction potential is also similar to interaction potential arising out of coupling between polar order and nematic splay in nematic liquid crystal 38 .For the symmetric case, J 1 = J 2 , i.e., ∆J = 0, the contribution of Ising-like alignment term vanishes and the interaction term associated with repulsion due to CIL essentially does not contribute to the energy in the bulk of the system and it features only as a boundary term .As a consequence, for a fully packed lattice comprised of particles with such interaction, the specific heat of the system would be zero at all temperatures and the particle-particle correlation length, ξ = 0. Of course, in the presence of vacancies, which allows for active transport of particles on the lattice, this interaction term plays a vital role in determining the nature of organization on the lattice.
For the asymmetric case, J = J 1 = −J 2 , the model reduces to a 1D Antiferromagnetic Potts model for J > 0 and to a 1D ferromagnetic Potts model for J < 0. However, the choice of J < 0 for the model Hamiltonian does not correspond to physically realistic scenarios for CIL.CIL-motivated interactions require the signs of J 1 and J 2 to be positive in order to ensure that configurations where the polarization vectors of neighbouring cells point towards each other are disfavoured, while the opposite holds when the polarization vectors of neighbouring cells point away from each other.Such choice necessarily breaks the Z 2 symmetry of the underlying Hamiltonian, and consequently distinguishes itself from the underlying symmetry of the Potts Hamiltonian.

Simulation Details
We perform Monte Carlo (MC) simulations of the system of N particles on L lattice sites starting with random configuration of particles uniformly distributed on a 1D lattice with fixed number density (ρ = N/L).For the initial configuration of particles on the lattice, we choose the polarization states of individual particles randomly with equal probability.We adopt a random sequential update procedure by choosing a site randomly with equal probability.For a site which is occupied by a particle, we perform the MC move for a particular process ( translation or directional switching) with the prescribed relative rates for the different processes.In order to ensure that the system settles to steady state, we wait for an initial transient time of 1000 L w steps , where w stands for the lowest rate ( among translation and switching).Thereafter we collect the statistics for the cluster sizes and other cluster characteristics, time averaging typically over at least 5000 samples.These samples are collected with a time spacing of 10 L w to ensure that the samples are uncorrelated.

The fully packed state: A reduced equilibrium model
In the experiments performed with MDCK cells on fibronectin coated strips, when all the adherent cells fill the strip corresponds to a confluent state 21 .In the confluent state, the dynamics at the boundaries of the adherent cells can lead to motility 21 .However, from the perspective of our minimal model, where we treat the particles as entities with hardcore repulsion, in the absence of vacancies, the state of the system is such that the translation dynamics of the particles is arrested and the dynamics of the particles is restricted to switching process between different states of polarization of the particles.In this regime, the system behaves as an effective equilibrium system, whose properties are determined by the Boltzmann weight associated to the Hamiltonian of Eq. 1 an effective temperature prescribed by the MC scheme.The corresponding form of the partition function, Z c , maybe expressed as, where the corresponding transfer matrix for Z c reads, Using standard techniques of Transfer Matrix method, in the thermodynamic limit of N → ∞, we obtain The corresponding expression for the average energy is, Fig. 2(a) displays the average energy per particle ε vs ∆J from analytical expression of Eq. 4 and its comparison with the simulation results.
While the average polarization is zero, the correlation function for the polarization as a function of particle distance, r, assumes the form, The corresponding expression for the correlation length is, In the limit of exp(∆J/2) >> 1, the correlation length ξ → 1 2 exp(∆J/2).Thus as would be expected for equilibrium 1D systems, there is no long range correlation of the polarization.However, for any finite lattice size system of size L, as long as ξ > L, the confluent state would tend to exhibit a polarized state.
In the presence of an external field h which couples with the individual particle polarization, the corresponding Hamiltonian assumes the form, The corresponding expression for the Gibb's Canonical Partition Function Z g is, Using this, the expression for the average polarization of the particle, σ is, Fig. 2(b) shows the excellent agreement of the variation of the average polarization with external field h obtained from Eq. 8 and with the simulation.

Clustering in the presence of vacancies
We next focus our attention on the nature of clustering behaviour of the particles that arises out of the interplay of translation dynamics of particles and the switching dynamics of particles.
The CIL interaction is controlled by the parameters J 1 and J 2 .While J 1 is the energy cost associated with the polarities of neighbouring particles pointing towards each other, J 2 is the reduction of energy associated with the polarities of the neighbouring particles pointing away from each other.In the absence of CIL between particles in a cluster, the switching rate from (→) to (←) is b and it is identical to the switching rate between (←) to (→).We define a dimensionless quantity Q, which is the ratio of the translation rate and which behaves as a measure of cell activity, a and the switching rate, b, and set b = 1 without loss of generality.In general, clustering will be controlled by the strength of CIL and Q, although the overall vacancy density is a conserved quantity, and it will also affect the cluster size distribution.In Fig. 3, we display the spatio-temporal evolution of the clusters for different activity strength quantified in terms of Q, for a fixed value of CIL interaction strength J 1 and holding J 2 = 0.At relatively high values of Q ( See Fig. 3b and Fig. 3c), the system segregates into alternating domains of dense clusters and a low density gas region.The mean sizes of these dense clusters increases on increasing Q.For the dense 6/15 clusters, the dynamics of their sizes is governed by the processes at the cluster boundaries.Typically, for large Q regime, the composition of the dense cluster is such that an array of right pointing (→) particles occupy the left end of dense cluster domain while the right end comprises of an array of left pointing (←) particles.Thus the internal structure within the bulk of such dense cluster comprises of defects-pairs of (→) and (←) particles.It may also be noted that such dense cluster for which the polarity of the particles at both ends are pointing inwards are immobile.At any instant of time, a cluster size can increase if a (→) particle from the adjoining gas region joins the left end of the dense cluster or if a (←) particle from the adjoining gas region joins the right end of the dense cluster.On the other hand, a cluster size can decrease if a single particle at the cluster ends switches it's polarity and subsequently leaves the cluster.
We next investigate the effect of CIL and strength of activity on the cluster size distribution.First we choose J 2 = 0, and modulate the CIL strength by variation in J 1 .Fig. 4 displays the cluster size distribution for different CIL interaction strength (J 1 ) and strength of activity Q.The limit -Q << 1, corresponds to a situation of low cell activity when particle translation rate is much slower than its polarity switching rate.For this case, in the absence of CIL interaction (J 1 = J 2 = 0), the probability of m-particle cluster is simply proportional to ρ m (1 − ρ) for N → ∞, where, ρ is the number density of particles on the lattice 36 .Consequently, the normalized probability of m-particle cluster is , P(m) = ρ m−1 (1 − ρ) and which maybe expressed in the exponential form, where ξ = | ln ρ | −1 , and the mean cluster size reads ⟨m⟩ = (1 − ρ) −1 .Eq. 9 is identical to the cluster size distribution of a Totally asymmetric exclusion process (TASEP) and Symmetric exclusion process (SEP) 36 .For low activity, increasing the CIL strength does not change the exponential nature of the cluster size distribution (See Fig. 4a).Although there is significant deviation from the exponential decay for small size clusters at large activity levels, Q >> 1, the asymptotic decay remains exponential as indicated by the log-plot of the distribution in Fig. 4(b), irrespective of the strength of CIL.As shown in Fig. 4(a) and Fig. 4(b), increasing the activity (Q) leads to an increase in the average cluster size, and to a broadening of the cluster size probability distribution both in the presence and absence of CIL.This generic trend comes from the accumulation of active particles in regions where they move slowly 11,15 .Due to excluded volume interaction, motility vanishes when a particle is obstructed by a neighbouring particle leading to formation of a cluster.The restoration of the mobility in a cluster is related to the cluster disintegration, which in turn depends on the switching rate of the particle's polarity at the cluster edges.Consequently, lowering the switching rate (hence increasing Q) increases the average cluster size.
To unravel the role of CIL on the average cluster size, we scrutinize separately the impact of J 1 and J 2 .Within a cluster, while J 1 favours polarity alignment, J 2 favours polarity anti-alignment of neighbouring particles, and accordingly, the interaction term of the Hamiltonian associated with CIL breaks the symmetry associated with independent rotation of the polarity vector.
When J 1 = J 2 , i.e., ∆J = 0, the CIL term in the Hamiltonian is identical to the CIL interaction considered on a phenomenological hydrodynamic theory for 1D cell assemblies 21 .In this regime, increasing CIL disfavours large cluster formation, leading to a monotonic decrease of the average cluster size with increasing CIL irrespective of the activity strength, Q.In this limit, the overall behaviour of the cluster size as a function of Q and J 1 maybe quantified in terms of the contour plot displayed in Fig. 5(a).
If J 1 ̸ = J 2 , which corresponds to an asymmetric CIL interaction, the effect of J 1 on the cluster size is more involved.Increasing J 1 at fixed J 2 has two effects on the clusters.While J 2 enhances the tendency of particles at the cluster edge to flip outwards, J 1 favours a polar alignment among neighbouring particles.This coupling favours overall that particles point inwards.Hence, J 2 promotes decrease in cluster size while J 1 has the opposite effect.Therefore, we can expect larger clusters to be stabilized when J1 > J2.In general, the competition between the opposite effects associated to J 1 and J 2 lead to a non-monotonic behaviour of the average cluster size, ⟨m⟩, as a function of J 1 beyond a threshold value of activity Q. Fig. 5(b), which displays the contour map of ⟨m⟩, for a fixed value of J 2 = 4, illustrates how the re-entrant like behaviour features for a wide range of Q and J 1 .The corresponding non-monotonic behaviour of ⟨m⟩ as a function of J 1 for different sets of Q, is displayed in Fig. 5(c).

Approximate expression for average cluster size
When the particle translation rate is much faster than the polarity switching rate, (Q >> 1), the system segregates into alternate regions of dense clusters (c) and a low density gas phase (g).In this regime, the interaction between the dense clusters is weak, and the stationary state is achieved as a balance between the incoming particle flux into the dense clusters from the surrounding gas region and the outgoing particle flux from the boundaries of the dense cluster region due to the particle switching their polarity at the cluster boundaries 13 .This balance can be estimated from an equivalent equilibrium process where the size distribution of dense clusters is determined by minimizing the effective Helmholtz free energy.It is worthwhile to point out that in the absence of CIL, the corresponding cluster size distribution is obtained by minimizing the system configurational entropy, as has been done for persistent exclusion process (PEP) 13 .

Minimization of effective Helmholtz Free energy (F)
Without loss of generality, for simplicity we incorporate the effect of CIL through J 1 , and set J 2 = 0.However, the procedure outlined to obtain the form of cluster size distribution can easily be generalized to J 2 ̸ = 0. From Eq. 4 we can express the mean energy of a cluster of length l as The configurational entropy S of the dense cluster phase corresponds to the number of ways by which Ω clusters can be arranged such that the clusters of same length are indistinguishable and are subject to the constraint that the total number of sites occupied by the clusters, N c , is constant.Then by fixing Ω, the entropy can be expressed as, where G c (l) is the number of clusters of length l in the cluster (c) phase and λ and γ are Lagrange multipliers, as already proposed for a PEP process 13 .Accordingly, the free energy, F, of the dense phase can be derived by accounting for the previous configurational entropy and the cluster internal energy due to CIL, leading to The corresponding expression for the variation of the free energy, δ F, reads and its minimization, δ F = 0 for independent variations of δ G c (l), provides cluster size distribution which has an exponential shape.
A similar argument for the gas(g) phase yields G g (l) -the number of clusters of length l in the gas(g) phase, with the corresponding form being, The corresponding parameters , A g , A c , l g and l c , that characterize univocally the coexisting cluster size distributions of the gas and dense phase can be determined imposing (a) The total number of clusters in the gas phase must equal total number of clusters in the dense cluster phase.This implies, (b) The total number of sites occupied by clusters in the gas phase together with the total number sites occupied by clusters in the dense cluster phase must equal total number of lattice sites, (c) If φ c and φ g are the number density of particles in dense cluster region and gas region, with ⟨l c ⟩ and ⟨l g ⟩ being the average length of the dense cluster and gas region, then the overall particle density ρ should obey, (d) If the hopping rate is much larger than switching rate, the gas region has typically a very low density of particles , hence φ g << 1.When a switching event at the boundary of a dense cluster occurs, a particle is emitted into the gas region.This leads to production of a dimer within the gas region which are the dominant clusters in the gas phase.Accordingly, we invoke the condition that in the limit of Q >> 1, the steady state is determined by the condition of matching the dimer production and disintegration rates 13 .The dimer production in the gas occurs when a particle at the edge of the adjoining dense cluster (with its polarity pointing inwards towards the bulk of cluster), has flipped its polarity and thus breaks away into the gas region, and within the time τ that it takes to reach the adjacent dense cluster, the particle at the edge of the adjacent dense cluster region has flipped its polarity.The average time that it takes for a particle in edge of the dense cluster to reach the edge of the adjacent dense cluster reads 9/15 Figure 6.Variation of the average cluster size, ⟨m c ⟩, in the cluster phase as a function of Q: Comparison of the theoretical expression, Eq. ( 26) ,with MC simulation.(a) J 1 = 0 (No CIL).(b) J 1 = 1 and (c) J 1 = 2.In all cases, J 2 = 0 and ρ = 0.8.Inset figures are the corresponding log-log plots.Solid lines correspond to expression in Eq. ( 26) while the circles correspond to MC simulations for each case.The dashed lines in the inset log-log plot corresponds to a slope of 1/2.MC simulation done with L = 2000 and averaging was done over 3000 samples.⟨τ⟩ = ⟨l g ⟩/a.Assuming that the neighbours of the edge particle of the dense cluster are pointing inwards, the rate of flipping of the particle at the edge is approximately equal to b.Hence, the overall dimer production rate, W p 2 , may be approximated by, where the contribution from trimer disintegration is neglected.
A typical configuration of a dimer would be (→)(←).The disintegration of a dimer would occur when either of the particles of the 2 particle dimer cluster switches their direction of polarity.While in the absence of CIL the switching rate is b, in the presence of CIL the switching rate is be J 1 /2 .Therefore In the steady state we equate Eq. 17 and Eq.18 to obtain the condition Since the particle flux from the gas into the cluster region, aφ g /2, must equal the particle flux from the dense cluster region, b (due to particle switching at the cluster boundary), at steady state, we arrive at

Mean cluster size
We approximate the sums in Eq.( 14), Eq.( 15) and Eq.( 19) by integrals to obtain an approximate expression for the mean cluster size.We set the integration limit for the dense cluster phase between 2 and ∞ since, for having a dense cluster, there needs to be at least 2 particles.For the gas phase, there has to be at least 1 site, which sets the lower integration limit.The upper integration 10/15 limit is taken to ∞ for N → ∞.With these conversions, we have, Substituting these expressions in Eq. (14-19), we obtain the algebraic relations involving A c , A g , l g and l c , A c l c e −2/l c = A g l g e −1/l g (21) (2 + l c )A c l c e −2/l c + (1 + l g )A g l g e −1/l g = N (22) Qe J 1 /2 A c e −2/l c = (1 + l g )A g l g e −1/l g (24) The previous approximation leads to, Using these relations, we arrive at the mean cluster size prediction Fig. 6 compares the analytic prediction of ⟨m c ⟩ with the simulation results for a range of CIL strengths.While in the absence of CIL (Fig. 6a), or weak CIL strength (Fig. 6b), the approximate analytical form matches the MC simulation results reasonably well, for moderate CIL strength (J 1 = 2), the deviation from the simulation result is more pronounced although even for this case ⟨m c ⟩ ∼ Q 1/2 when Q >> 1 (See Fig, 6c).However for very strong CIL interaction strength, the approximate analytical curve fails to reproduce average cluster size characteristics.

Cluster Polarisation
CIL has a strong impact not only on host particle aggregate in clusters, but also on how they align relative to each other.The net polarization, of a cluster of size m, is quantified by P m = ∑ m i=1 σ i .However, since clusters generically do not develop a net polarity, we analyze their root-mean square (RMS) fluctuations, In the absence of CIL when Q << 1, particle's polarization are uncorrelated from the neighbouring ones.In this regime, as displayed in Fig. 7a(i), S F = m −1/2 , corresponding to the RMS of a binomial distribution.For large activity, Q >> 1, in the absence of CIL, for relatively small cluster sizes, there is a negative deviation of S F compared to m −1/2 , although for sufficiently large clusters, it merges with the form corresponding to binomial distribution as displayed in Fig. 7a(ii).In the presence of CIL for which J 1 is high, and J 2 = 0, for range of cluster sizes, S F does not scale as m −1/2 and instead remains virtually unchanged irrespective of whether activity is high or low as displayed in Fig. 7a(iii) and 7a(iv)) respectively.However the cases of high and low activity are distinguished by the fact that the formation of large size clusters is much more prevalent when activity is high as compared to the case when activity is low.When J 2 ̸ = 0 and J 1 = 0, CIL induces "anti-ferromagnetic" like order within the cluster.This ordering is captured through the "anti-ferromagnetic"-like order parameter defined in terms of fluctuation of the difference of polarization of two sublattices, A and B within a cluster of size m, 11/15 where the summation over i for A sublattice is over odd sites and for B sublattice it is over the even sites within a cluster.In the absence of CIL interaction, when Q << 1, the polarization of the individual particles are uncorrelated with the polarization of the other particles constituting the cluster.Hence in this limit, both the ferromagnetic and anti-ferromagnetic order parameter, e.g., S F and S AF are identical to RMS fluctuation of a binomial distribution with S F = m −1/2 as displayed in Fig. 7b(i) and 7b(ii).On the other hand when Q >> 1, when J 2 is high (and J 1 = 0), the anti-ferromagnetic order parameter, S AF has a positive deviation from m 1/2 for all ranges of cluster size, while the ferromagnetic OP (S F ) has a negative deviation as shown in Fig. 7b(iii) and 7b(iv).

Effect of external field on cluster characteristics
We incorporate the effect of an external field h o whose effect is to aid translation of right end directed (→) particles and oppose motion of left end directed (←) particles.Effectively the hopping rate of the (→) particle becomes a + h o , while for (←) particles the hopping rate is a − h o .As long as h o < a, the natural direction of movement of both (→) and (←) particles is retained.When h o = a, the − ended particles do not move.For h o > a, both (→) and (←) particles move in the same direction, as prescribed by the external field.For relatively low activity, when h o < a, increasing h o has the effect of marginally increasing the average cluster size.In general, the cluster size is not significantly altered ( See Fig. 8a).
As shown in Fig. 8(b), for large activity the average cluster size varies non monotonically with h o .Initially, the average cluster size decreases, but after a threshold, the mean cluster size starts increasing with the external field magnitude.However for sufficiently high strength of external field h o , the average cluster size again starts increasing and the corresponding distribution function of the cluster size starts broadening with h o .As h o approaches the hopping rate a, average cluster size sharply increases.Finally when h o = a, for which the (←) particles stop moving, the cluster size distribution becomes broadest, with average cluster size increasing drastically (around 2.5 times the average cluster size for h o = 0).As shown in Fig. 8(c), when h o is increased further, h o > a, not only do (→) particles and (←) particles move in the same direction, the corresponding average cluster size monotonically decreases with increase in h o .

Discussion
In this paper we have presented a minimal discrete driven lattice gas model which mimics the phenomenology of CIL interactions between cells and the movement of the cells in confined in 1D channel.In the absence of vacancies ( akin to dense packing of cells in 1D array), the translation dynamics is arrested.In this limit, the model reduces to an equilibrium spin model which does not possess Z 2 symmetry like Ising model.We solve this model exactly both in the presence and absence of an external field which couples to the individual polarization of the cell.In the presence of vacancies, the interplay of translation dynamics and CIL interaction between particles results in a steady size cluster size distribution which is exponentially distributed for the large size cluster.The typical cluster size and the distribution function is controlled by CIL strength and activity.The effect of increasing activity at fixed CIL strength invariably leads to an increase in the average cluster size.However the effect of varying CIL interaction ( through J 1 at constant J 2 ), can result in a non-monotonic dependence of the average cluster size as a function of CIL strength J 1 .In the high activity regime, Q >> 1, an analytic form of average cluster size can be obtained approximately by effectively mapping the system to an equivalent equilibrium process involving of clusters of sizes wherein the cluster size distribution is obtained by minimizing an effective Helmholtz free energy.The resultant prediction of exponential dependence on J 1 of the average cluster size and Q 1/2 dependence of the average cluster size is borne out to reasonable accuracy as long as J 1 is not very large.At very high values of J 1 , the prediction breaks down, indicating the failure of this approximation scheme.The polarization characteristics within an cluster is indicated by emergence of a local ferromagnetic like order parameter within a cluster due to the effect of CIL, when J 1 > J 2 .In the presence of an external field h o which couples to translation dynamics of individual system, system, such that it aids translation of right end directed (→) particles and opposes motion of left end directed (←) particles, we find that the the average cluster size exhibits a non-monotonic dependence on h o .While many studies of active particles systems have predicted and reported existence of phase transitions due to Motility induced Phase separation(MIPS) 15,16 , we have not observed any such phase transitions for our model for the range of parameters that we have explored.
While in our present work we have focused solely on the interplay of cell movement and CIL interaction between cells, and as such treated the cells as hard objects, its worthwhile to point out that apart from re-polarization events, adhesion of the colliding cells and the relatively rare event of cells walking past each other has also been observed 7,8 .Cells are contractile entities which exert forces on each other and on the underlying substrate.It has been observed that cells on contact with each other make use of trans-membrane adhesion molecules such as cadherins to adhere to each other 9 .Presence of these adhesion molecules serve to generate an effective attractive force between the cells.Thus a more holistic understanding of the collective organization of cells would require inclusion of attractive interaction between the cells, additional alignment mechanisms that might be at play apart from the effect of CIL that we have discussed.In our future work we seek to incorporate these aspects in the framework of our modeling approach.Finally we also seek to generalize our discrete lattice gas modeling approach to study the effect of CIL in context cellular organization on 2D substrate.

Figure 1 .
Figure 1.Schematic representation of the relevant dynamical process: (a) Translation process.(b) Polarity switching of a particle at the boundary of a cluster.A single particle bounded by vacancies does not switch its polarity.

Figure 2 .
Figure 2. (a) Variation of the energy per particle, ε = E/N as a function of ∆J = J 1 − J 2 .ε is expressed in units of k B T , with T being an effective temperature of the active system.(b) Variation of average polarization per particle, σ as a function of external field h for ∆J = 2. Solid lines correspond to the analytical expression while the points corresponds to results obtained using MC simulations for N = 1000.

7 / 15 Figure 5 .
Figure 5. (a) Contour plot of average cluster size ⟨m⟩ as function of J 1 and Q when J 1 = J 2 i.e. ∆J = 0. (b) Contour plot of average cluster size ⟨m⟩ as function of J 1 and Q for constant J 2 (J 2 = 4): Re-entrant like behaviour is observed for Q > Q c with Q c = 12.5.(c) Plot of ⟨m⟩ with J 1 corresponding to (b).For all cases, MC simulations where done with L = 1000 and averaging was done over 2500 samples at ρ = 0.8.

Figure 7 .
Figure 7. (a) RMS Fluctuation of polarization, (S F ) with cluster size (m): (i) No CIL and low Q; Q = 0.1, (ii) No CIL and high Q; Q = 30, (iii) J 1 = 7 and low Q; Q = 0.1, (iv) J 1 = 7, high Q; Q = 30.Here.J 2 = 0.The binomial distribution corresponds to solid black line.(b) RMS Fluctuation of local antiferromagnetic order parameter,S AF , and ferromagnetic order parameter, S F , in a cluster of size m vs Cluster size (m) : (i) S F for No CIL and low Q; Q = 0.1, (ii) S AF for No CIL and low Q; Q = 0.1, (iii) S F for J 2 = 7 and high Q; Q = 30, (iv) S AF for J 2 = 7, high Q; Q = 30.Here, J 1 = 0. MC simulations were done with L = 1000 and averaging over 2000 samples.