Segregation in binary mixture with differential contraction among active rings
Emanuel F. Teixeira Carine P. Beatrici Heitor C. M. Fernandes Leonardo G. Brunnet





emanuel.teixeira@ufrgs.br
Instituto de Física, Universidade Federal do Rio Grande do Sul, CP 15051, CEP 91501-970 Porto Alegre - RS, Brazil


carine@if.ufrgs.br
Instituto de Física, Universidade Federal do Rio Grande do Sul, CP 15051, CEP 91501-970 Porto Alegre - RS, Brazil


heitor.fernandes@ufrgs.br
Instituto de Física, Universidade Federal do Rio Grande do Sul, CP 15051, CEP 91501-970 Porto Alegre - RS, Brazil

leon@if.ufrgs.br
Instituto de Física, Universidade Federal do Rio Grande do Sul, CP 15051, CEP 91501-970 Porto Alegre - RS, Brazil




§ ABSTRACT 
Cell cortex contraction is essential for shaping cells, enabling movement, ensuring proper division, maintaining tissue integrity, guiding development, and responding to mechanical signals — all critical for the life and health of multicellular organisms. Differential contractions in cell membranes, particularly when cells of different types interact, play a crucial role in the emergence of segregation. In this study, we introduce a model where rings composed of active particles interact through differential membrane contraction within a specified cutoff distance. We demonstrate that segregation arises solely from differential contraction, with the activity of the rings functioning similarly to an effective temperature.
Additionally, we observed that segregation involves cluster fusion-diffusion process. However, the decay exponent of the segregation parameter we found is close to λ∼ -1/3, which differs from the λ∼ -1/4 predicted by previous theoretical approaches and simulations.




Segregation in binary mixture with differential contraction among active rings
    Leonardo G. Brunnet
    September 9, 2024
==============================================================================




In cellular systems made up of different species, spontaneous sorting is a common emergent behavior.
During embryonic development, cells undergo differentiation, which leads to the spontaneous segregation of cell types in tissue formation. Based on strong experimental evidence <cit.> and motivated by physical systems like binary mixtures, where spontaneous separation of two liquids is observed,  Steinberg proposes <cit.> that the general mechanism for cell segregation lies in the difference in adhesion between cells of  different  types.
This proposal became known in the literature as the Differential Adhesion Hypothesis (DAH) or Steinberg Hypothesis. 




Harris <cit.> questions the foundation of the DAH hypothesis, highlighting that mere maximization of intercellular adhesion does not necessarily lead to the observed effects of cell sorting. He suggested that within a heterogeneous cell aggregate, variation in surface contraction, driven by the active regulation of the acto-myosin cortex, could be the driving force behind tissue engulfment and cell sorting in vivo.
The surface contraction of a particular cell appearing more pronounced when it comes into contact with a cell of a different histological type. 
Besides that, cells from different tissues exert different surface contraction when they contact the medium. This sorting mechanism was named Differential Surface Contraction Hypothesis (DSCH). 

In line with Harris's ideas, Brodland <cit.> successfully described segregation using finite element simulations that incorporated adhesion and surface contractions, both contributing to the total interfacial tension, a mechanism known as differential interfacial tension hypothesis (DITH).

In fact, these hypotheses have been tested using various other extended numerical models, including Cellular Potts <cit.> and Vertex models<cit.>. 
All theses models define an effective interfacial tension between different units but do not explicitly separate adjacent cell membranes or detail the forces involved, such as adhesion and surface contraction. 
However, from an experimental perspective, the role of membrane fluctuations was clearly highlighted by Mombach <cit.>, and the distinction between adhesion and cortical tension was emphasized in the work of Krieg et al.<cit.> and Manning et al. <cit.>.


In this work, we present a model of active rings that interact through differential membrane contraction when within a specific cutoff distance. By utilizing two interacting membranes, this model provides a more detailed representation of biological processes, offering insights into how different layers interact and differentiating the roles of fluctuations, adhesion, and cortex contraction. To our knowledge, this is the first time Harris's mechanism has been simulated in isolation.




Boromand and collaborators <cit.> introduced a ring system composed of passive particles connected by springs. 
Building on this model, we incorporated active properties into the particles forming the ring in previous articles <cit.>. We model a 2D system with N rings, each representing a cell formed by n active particles (see Fig. <ref>).
The system is set up as a binary mixture of active rings confined in a circular arena with repulsive walls and radius R_0. 
The set of coupled overdamped equations governing the dynamics of each particle are: ṙ_i,j = v_0 𝐧_j -μ 𝐅_i,j and ṅ_j = √(2 D_R) ξ_j×𝐧_j,
where 𝐫_i,j denotes the position of the  i-th particle within the j-th ring at time t, μ represents its mobility and v_0 is the magnitude of active velocity with its orientation given by 𝐧_j. 
The second term, 𝐅_i,j = -∇_i,j E represents the total force acting on the i-th particle within the j-th ring.
The direction of the active force, described by unit vector 𝐧_j, experiences angular Gaussian white noise ξ_j = ξ_jê_z with correlation ⟨ξ_j(t_1)ξ_k(t_2) ⟩ = δ_jkδ (t_1-t_2). The noise term D_R represents the rotational diffusion coefficient, which defines a characteristic timescale given by τ_R = 1/D_R.

The energy function has contributions from (i) a perimeter energy (springs connecting ring neighboring particles), (ii) area conservation, (iii) contact-dependent contraction term, (iv) core repulsion among non-neighboring particles from any ring and inter-cellular adhesion, 

    E    =   ∑_j=0^N{ϵ_P/2∑_i=0^n( |l⃗_i,j|/l_0 - 1 )^2 + ϵ_A/2( A_j/A_0 - 1 )^2. 
       +   . ∑_i=0^nΛ_αβ |l⃗_i,j| } + ϵ_c/2∑_r_ik≤σ( r_ik/σ - 1 )^2
       +   ϵ_adh/2∑_σ < r_ik≤  l_adh( r_ik/σ - 1 )^2,

where l⃗_i,j = r⃗_i,j - r⃗_i-1,j is the vector connecting consecutive particles in ring j, ϵ_P is the elastic energy of the spring controlling perimeter fluctuations, and l_0 is the equilibrium distance in the ring. 
The elastic energy related to area control is ϵ_A, A_j is the inner area of the j-th ring (not considering the area of particles), and A_0 is the equilibrium area. The last two terms in Eq. <ref> represent the core repulsion between non-neighboring particles  and the adhesion between particles of different rings, respectively. Parameter ϵ_c denotes the characteristic energy of the core repulsion interaction, while r_ik is the distance between particles i and k. The equilibrium cut-off distance, σ, effectively defines the particle diameter. The adhesion energy and interaction distance are characterized by ϵ_adh and l_adh, respectively.




A key distinction between the ring model and other extensive models (Finite Elements <cit.> and Vertex-Voronoi model  <cit.>) lies in the rings contact interface, which comprises two contact membranes. 
This bring us to the third term in Eq. <ref> which incorporates differential surface contraction. In this term, parameter Λ_11 represents the line tension between particles belonging to type 1 rings (red rings in Fig. <ref>), Λ_22 represents the line tension between type 2 rings and Λ_12=Λ_21 represents the line tension between rings of different types. All these tensions acting when they are within a cutoff distance l_Λ.
A central parameter in this work is Λ = Λ_12/Λ̅, with  Λ̅≡ (Λ_11+Λ_22)/2.
Fig. <ref> shows a schematic representation of the active ring system, illustrating  the involved tensions. 
 
To handle the interaction with the walls, which can be seen as the medium in our context, we specify that the particles of a ring increase their contraction tension once they reach a distance σ/2 from the wall. 
This prevents the rings from preferring to accumulate on the wall. 
The contraction between rings of type 1 (red) and the “medium” (wall) is defined as Λ_1M/Λ̅ = 12, while for rings of type 2 (green), we use Λ_2M/Λ̅ = 0. 
This choice  satisfies the criterion Λ_12< Λ_1M-Λ_2M (valid for all values of Λ_12 in this work) for the engulfment of one type (red) by the other (green) <cit.>.


We measure spatial coordinates and time in units of l_0 and τ_R, respectively. 
As a result, the model equations may be written in terms of the Péclet number <cit.>, defined as Pe ≡v_0 τ_R/l_0 , which controls the level of activity. Details of the numerical integration method and the remaining parameter values are provided in the Supplemental Material.

We begin our simulations with an initial random distribution of rings (see Fig. <ref>a). For Λ > 1, we observe a spontaneous segregation process among the rings, with larger Λ values resulting in a smaller interface between ring types, indicating improved segregation. This is illustrated at the bottom of Fig. <ref>a, where Λ = 10. For Λ = 1.5 (middle of Fig. <ref>a), the system still segregates, but a larger asymptotic interface remains. In contrast, for Λ < 1 (top of Fig. <ref>a), an organized mixed state emerges, forming a checkerboard-like pattern where each ring tends to be in contact with a ring of a different type.


 



To quantify the level of segregation in the system, we use the parameter γ, as introduced by Belmonte and collaborators <cit.>. This parameter is defined as the mean fraction of neighboring rings of type 2 (green) surrounding rings of type 1 (red): γ = ⟨n_2/n_1+n_2⟩_1, where ⟨..⟩_1 denotes the average over all rings of type 1, with n_1 and n_2 representing the number of first neighbors of type 1 and type 2, respectively. In the literature, it is well established that the segregation parameter is directly correlated with the interface I between the types, expressed as γ∼ I <cit.>.

We examined the evolution of the segregation parameter γ for various Λ values while keeping the Pe value fixed. As shown in Fig. <ref>b), the asymptotic value of γ reaches low levels (approximately γ∼ 0.1 with 10^3 rings) for the segregated state and increases to γ∼ 0.9 as Λ decreases. The line separating the segregated and checkerboard patterns occurs at Λ∼ 1, where the initial value of the mixed state (γ∼ 0.7) is maintained throughout the evolution.





Furthermore, we measured mean values of steady-state γ as a function of Λ for several values of Pe (see Fig. <ref>a).

 Here we emphasize the role of parameter Pe, which broadens the transition as its value increases, and the importance of initial configurations to avoid stay trapped in a meta-stable state. 
While at low Pe values (Pe ≤ 0.2 in Fig. <ref>a) and Λ <1, the system evolves to the checkerboard pattern, for Λ>1 it remains trapped close to the initial configuration, so we employed the segregated state as the initial condition to test the effect of Pe. 
In brief, Pe functions like temperature, enabling the system to overcome potential barriers to approach the minimum energy state. 

Simultaneously, it is responsible for deviating from the optimal value due to the introduced fluctuations.


In Fig. <ref>b we present a (Λ× Pe) diagram, colors indicate the mean value for γ. 

At low Pe values, the segregation region occurs just above Λ = 1, but  at high activity we observe a mixed state well above this limit. 
This indicates that the segregation criterion <cit.> Λ_12 > (Λ_11+Λ_22)/2 is only valid at low Pe. 
Similar arguments apply to the checkerboard state in the region defined by Λ<1.




Additionally, we conducted an analysis similar to previous works <cit.> and observed how the steady state of largest cluster of type 1 rings (red) changes with Pe and Λ.
We define the fraction of these rings that are inside the largest cluster as c_in = N_int/N_1. N_1 and N_int are the total number of type 1 rings and the number of rings in the largest asymptotic cluster, correspondingly. 
Using c_out = 1 - c_in, we can construct a binodal curve that delineates regions of segregation and mixing as function of Λ^-1 and Pe (see Figures <ref>c-d). 
Below the intersection we have the segregated state (c_in∼ 1 and c_out∼ 0, one large cluster) and above, the mixed state (several clusters of similar size). 
Fig. <ref>d shows the effect of the activity parameter Pe on the binodal curves. 
Depending on the value of Λ, there is a  value of Pe beyond which the system transitions from segregated to mixed state.




 
 
 
 
 
 
 
 

 
 
 
 










Our focus now shifts to understanding how the system evolves from a randomly mixed configuration to a segregated one. 
This is particularly important because the long time behavior before saturation of the segregation parameter and the mean size growth of clusters can provide information about the underlying mechanisms occurring during tissue formation. 
The mean cluster's size of rings of type 1, M(t), is obtained using the cluster-counting algorithm used by Beatrici et al. <cit.>.

Here, we analyze the case where contact between  different rings type is unfavorable (Λ = 10). 
Initially, the system is in a mixed configuration with a corresponding value of γ≡γ_0=0.7, consequent of the (30:70) proportion of ring types. 
In Fig. <ref>a, we show the evolution of the segregation parameter normalized by its initial value, γ̅(t) = γ(t)/γ_0. 


The activity Pe determines the duration of the transient period, with higher activity implying shorter transients.
Subsequently we observe an algebric decay, γ∼ t^λ with an exponent close to -1/3 for at least two decades, independent of the activity value Pe. 
Additionally, in this regime, the mean domain size M grows with an exponent close to 2/3 (see Fig. <ref>c). 
Thus, we observe an inverse relationship, γ∼ M^-2, consistent with previous findings <cit.>.




In Fig. <ref>b, we observe that significantly unequal proportions slow down the domain growth process, as it relies on the coalescence of minority clusters diffusing within the majority type. However, when the proportions are equal (50:50), the dynamics change — the system quickly percolates, and the process of interface reduction becomes dominated by rounding.
In this case, we expect exponent ∼ -1/2 <cit.>. For other proportions, we observe asymptotic time exponents close to -1/3 and 2/3 for γ and M, respectively, while maintaining the inverse relationship, γ∼ M^-2 (see Fig. <ref>b-d).
 
 

















In conclusion, we present a model where a cell is represented by a ring of active particles. Each ring is explicitly equipped with its own membrane, allowing it to accommodate negative curvatures and replicate all cell sorting mechanisms, whether in isolation or combined. This represents a significant advancement in the class of extended cell models, offering a path to more realistic depiction of cell behavior as observed in experiments.

By employing this comprehensive cellular model and incorporating a term for differential contraction among rings of different types, we conducted a numerical evaluation of Harris' proposed hypothesis on cellular differential surface contraction <cit.>. 
We found that by keeping all cells with identical attraction forces and solely employing adequate differential interfacial contraction, the system segregates.
We observed a monotonic relationship between the differential contraction parameter Λ and the steady-state average value of the segregation parameter ⟨γ⟩_t, indicating that segregation occurs for Λ > 1 and a checkerboard pattern emerges for Λ < 1. 

The rate at which the system transitions from an initially mixed state to its asymptotic state, as well as the ultimate value it reaches, also depends on the activity level Pe.
In the case of segregation (Λ>1), increased activity has a tendency to disrupt the configuration of minimum interfacial energy.
 Binodal curves, showing a transition between the mixed and the segregated states, corroborate the inverse correlation between the contraction interaction Λ and the activity Pe. 
 These findings suggest that the activity Pe plays a role similar to temperature in thermodynamic equilibrium systems <cit.>.
  
  Furthermore, we found that the segregation parameter γ(t) and the mean cluster size M(t) exhibit a  power-law regime with an exponent  close to -1/3 and 2/3, respectively. 
  The activity Pe changes the typical timescale at which the asymptotic regime begins but does not modify the growth exponent of the domains. 
Similarly, the proportions do not alter the characteristic exponent, except in the (50:50) case where the system starts close to percolation and the evolution appears to be dominated by rounding. We observed an inverse relationship between the segregation parameter and the typical size, γ∼ M^-2.


Finally, the observed segregation exponent of λ=-1/3 is at odds with the literature. The cluster-cluster aggregation mechanism seen in our simulations would typically correspond to an exponent of λ=-1/4, as predicted by surface diffusion models <cit.> or mean cluster models <cit.>, particularly in the absence of alignment interactions. An exponent of λ=-1/3 would instead be expected from the Cahn-Hilliard equation if the mechanism were evaporation-condensation, which we do not observe. The reasons for these discrepancies will be explored in future studies.


We express our gratitude to the Brazilian agencies CAPES, CNPq, and FAPERGS for their financial support. H.C.M.F. and L.G.B. acknowledge the support from the National Council for Scientific and Technological Development – CNPq (procs. 402487/2023-0 and 443517/2023-1). E.F.T. acknowledges ICTP-SAIFR/IFT-UNESP. The simulations were conducted using the https://pnipe.mcti.gov.br/laboratory/19775VD Lab cluster infrastructure at IF-UFRGS.



§ SUPPLEMENTARY MATERIAL



§ A - BOUNDARY CONDITIONS - REPULSIVE CIRCULAR WALLS

We configure the system as a binary mixture of N active rings confined within a circular arena of radius R_0. The force exerted by the wall on a specific particle within a ring is given by


    F⃗_w =  -F_w (|r⃗_i,j-r⃗_cm|/R_0 - 1  )  ( r⃗_i,j-r⃗_cm  )/|r⃗_i,j-r⃗_cm|

where ( r⃗i,j-r⃗cm )/|r⃗i,j-r⃗cm| is the unit vector connecting the center of particle i within ring j to the center of mass, and F_w is the characteristic force exerted by the wall on the particle.


 §.§ B - Forces



We derive the analytical expressions for the forces acting on each particle i within a ring j. These particles are subjected to forces resulting from both shape and interaction energies, which are determined by the area, perimeter, differential line tension, particle-particle overlap, and adhesion terms as defined in the main text. The force on particle i in ring j is obtained by taking the vector derivative of the total energy E with respect to its coordinates.


    r⃗_i,j   =    x_i,jê_x + y_i,jê_y,  
    
    F⃗_i,j   =    -∂ E/∂r⃗_i,j≡ -∂ E/∂ x_i,jê_x-∂ E/∂ y_i,jê_y.



  §.§.§ Perimeter force


The perimeter energy, as defined in the main text for a configuration of N rings labeled by j = 1, …, N, with n particles labeled i = 1, …, n, is given by

    E_P = ϵ_P/2∑_j=0^N∑_i=0^n( |l⃗_i,j|/l_0 - 1 )^2.

The force on particle i due to deviations in the segment length |l⃗_i,j| = |r⃗_i,j - r⃗_i-1,j| from its preferred value l_0 is given by

    F⃗^i,j_P   =    -∂ E_P/∂r⃗_i,j = -∂ E_P/∂ x_i,jê_x - ∂ E_P/∂ y_i,jê_y, 
    
    F⃗^i,j_P   =    -ϵ_P/l_0 [  (|l⃗_i,j|/l_0-1   )∂ |l⃗_i,j|/∂r⃗_i,j  -  (|l⃗_i+1,j|/l_0-1   )∂ |l⃗_i+1,j|/∂r⃗_i+1,j ], 
    
    F⃗^i,j_P   =   ϵ_P/l_0 [  (|l⃗_i+1,j|/l_0-1   ) l⃗_i+1,j/|l⃗_i+1,j| -  (|l⃗_i,j|/l_0-1   )l⃗_i,j/|l⃗_i,j| ],

where we use the following relations:

    ∂ |l⃗_i,j|/∂r⃗_i,j   =   ∂ |l⃗_i,j|/∂ x_i,jê_x+  ∂ |l⃗_i,j|/∂ y_i,jê_y = -l⃗_i,j/|l⃗_i,j|, 
    
    ∂ |l⃗_i+1,j|/∂r⃗_i+1,j   =   ∂ |l⃗_i+1,j|/∂ x_i+1,jê_x+  ∂ |l⃗_i+1,j|/∂ y_i+1,jê_y = l⃗_i+1,j/|l⃗_i+1,j|.



  §.§.§ Area force

The force related to area energy on particle i in ring j is given by

    F⃗^i,j_A= -∂ E_A/∂r⃗_i,j = -ϵ_A/A_0 ( A_j/A_0 -1  )∂ A_j/∂r⃗_i,j,

being E_A the area energy defined by

    E_A = ϵ_A/2∑_j=0^N( A_j/A_0 - 1 )^2.

The ring area A_j is calculated using the vector product property, given by


    A_j   =   1/2∑_j^n |(r⃗_i,j-r⃗_cm)×l⃗_i,j| 
    
       =   1/2∑_j^n (x_i,j- x_cm)(y_i,j-y_i-1,j) 
    
       -   (y_i,j-y_cm)(x_i,j-x_i-1,j),

where the factor 1/2 is included to avoid double-counting the area. Therefore, we obtain:


    ∂ A_j/∂ x_i,j   =    ( y_i+1,j -y_i-1,j )/2, 
    
    ∂ A_j/∂ y_i,j   =    ( x_i-1,j -x_i+1,j )/2,

thus, 

    F⃗^i,j_A   =    -∂ E_A/∂r⃗_i,j = -  ϵ_A/2A_0 ( A_j/A_0 -1  ) 
    
       ×    [  ( x_i+1,j - x_i-1,j ) ê_x -  ( y_i+1,j - y_i-1,j )ê_y ].



  §.§.§ Interaction forces


The core repulsion for non-neighboring particles within the same ring and the adhesion between particles of different rings are accounted for by the interaction energy E_int, which is described by a truncated harmonic potential,


    E_int   =   {[   ϵ_c/2( r_ik/σ - 1 )^2,                                            r_ik≤σ; ϵ_adh/2( r_ik/σ - 1 )^2,                                   l_adh≥ r_ik > σ;                       0,                                     r_ik > l_adh. ].

Therefore, the corresponding force

    F⃗^i,j_int   =    -∂ E_int/∂r⃗_i,j, 
    
    F⃗^i,j_int   =   r̂_ik{[   -ϵ_c/σ( r_ik/σ - 1 ),                                          r_ik≤σ; -ϵ_adh/σ( r_ik/σ - 1 ),                                 l_adh≥ r_ik > σ;                      0,                                    r_ik > l_adh ].

where r̂_ik = r⃗_ik/r_ik is a unit vector connecting particle i and k.


  §.§.§ Differential line tension (contraction)

The contraction energy represents a line tension that acts between two particles from different rings when they are within a distance cutoff l_Λ. This energy term is given by


    E_Λ =  ∑_j=0^N∑_i=0^nΛ_αβ |l⃗_i,j|.

The force acting on particle i inside ring j is

    F⃗^i,j_Λ   =    -∂ E_Λ/∂r⃗_i,j, 
    
    F⃗^i,j_Λ   =   Λ_αβ [  l⃗_i+1,j/|l⃗_i+1,j| - l⃗_i,j/|l⃗_i,j| ].

Therefore, the force F⃗^i,j_Λ on particle i within ring j depends on the type α of its ring and the type β of the ring of the particle it interacts with. If the particle is within the cutoff distance and interacts with multiple particles from other rings, the tension value will be the mean of all the individual tensions. Specifically, this mean is given by

    Λ_αβ = 1/n_c∑_k^Λ_α k,

where the sum k runs over all contacts (both α-α and α-β interactions) and n_c represents the total number of contacts.




 §.§ C - Control parameters


The value of the total equilibrium area A_r is the sum of the equilibrium area imposed by the area elastic energy term   plus half the area of each particle composing the ring, 
 
    A_r = A_0 + nπσ^2/8.

We define a packing fraction ϕ for the system, relative to a circular region 
of radius R_0,  through the relation 

    ϕ = NA_r/π R_0^2.

We fix ϕ =  0.895  allowing for neighbor exchange and neighbor interaction simultaneously and we use a fixed dimensionless perimeter-area equilibrium ratio p_0 = P_0/√(A_0) = 4, where P_0 = nl_0 is the equilibrium perimeter. 
The shape parameter p_0 defines the degree of stiffness of the ring. 
The choice of this value is based on previous works with Vertex <cit.>, Voronoi <cit.>, and ring models <cit.>, where an excess of cell perimeter is found for p_0 > 3.81 as well as the emergence of a liquid-like behavior.  To prevent overlap among rings, we reached a compromise by setting comparable values for ϵ_c, ϵ_P and ϵ_A/n, while assigning a much lower value to ϵ_adh. 
So, we simulate the ring system keeping the following parameters fixed:  ϵ_c/ϵ_P = 1, ϵ_A/ϵ_P = 35, (F_wσ)/ϵ_P = 1, ϵ_adh/ϵ_P = 5.10^-4, n = 10, l_0= 1, σ = l_0, l_Λ = l_adh = 1.5 l_0 and μ = 1. 
This choice of parameters ensures that ϵ_P and ϵ_A are sufficiently large to maintain the bond length close to l_0 and the equilibrium area close to A_0.  The adhesion forces between different rings are kept equal for all rings. 
Furthermore, to emphasize the effects of differential contractions, we adopted an adhesion value that is considerably lower compared to the other energy terms involved. We integrate the equations of motion, using the Euler-Maruyama algorithm with a time step Δ t = 0.01.
 


 §.§ D - Neighbors criterion for measurements

In this work, we use the criterion for defining whether two rings are neighbors based on the distance between their centers of mass, d_jk = |R⃗_j - R⃗_k| ≤P_0/2 = 5σ. This relation ensures that even two completely flattened rings in contact will be considered neighbors.

apsrev4-2



