We perform Monte Carlo simulation to explore the energy landscape and the spatial organization of the force dipoles emerging from the interactions between force dipoles (Eq. 5). We consider a square lattice of size L × L where each lattice site (r _i ) contains a point force dipole characterized by a polarization magnitude P(r _i ) and an orientation angle ϕ _i . The polarization magnitude is given by the product of the dipole length and the dipole force, P (r _i ) = dF (r _i ). In our simulations, the dipole length d is taken to be the same for all dipoles. We also choose d (≪ l), the lattice spacing, to mimic the point dipole in Eq. 5. The dipole force, F (r _i ) is chosen from a distribution, as detailed below.

Several biological cell lines are known to have their aspect ratios distributed according to a k − Γ distribution. We replicate this in our simulations by drawing polarization magnitudes (P) from a distribution characterized by parameters k and θ, and described by probability density function f P = 1 Γ k θ k P k − 1 e − P / θ . (7) The orientation angles (ϕ) are randomly chosen between 0 and π. The interaction between force dipoles is given by Equation 5, where a dipole’s nematic tensor q α β ( r i ) = n ̂ α n ̂ β can be obtained by writing the dipole orientation vector as n ̂ = ( cos ⁡ ϕ , sin ⁡ ϕ ) . The total energy is obtained by summing the long range interactions between all pairs of dipoles with periodic boundary condition on all sides.

To analyze the energy landscape, we access a multitude of metastable states by thermalizing the system at a high temperature and then running the simulation at a temperature that is low compared to the energy of the system: k B T E ≈ 1 0 − 6 . For convenience, we will often use a dimensionless form of the energy given by ε = E l 3 N ⟨ P ⟩ 2 d 2 , where E is the dimensionful energy (Eq. 5), l is the lattice spacing, N is total number of dipoles, ⟨P⟩ is the mean value of polarization magnitude, and d is the dipole length.

During each Monte Carlo step, the two properties of each dipole, polarization and orientation, are updated following two different rules. For the orientation, we use model A (non-conserved) dynamics (Hohenberg and Halperin (1977)), where a new orientation angle is proposed randomly between 0 and π. The polarization magnitudes are updated following model B (conserved) dynamics (Hohenberg and Halperin (1977)), where the polarization magnitudes of two randomly chosen dipoles are exchanged. This dynamics guarantees that the overall distribution of the polarization does not change during the simulation, however, the spatial organization of the magnitudes evolves along with the orientation of the dipoles. We draw the polarization magnitudes of our force dipoles from a k-distribution for each of our initial conditions. Individual Monte Carlo runs then lead to a sampling of the energy landscape corresponding to a particular k-distribution. Statistical properties can then be inferred by averaging over initial conditions. Acceptance rate of a proposed update is determined by the Boltzmann factor exp ( − Δ E K B T ) , ΔE being the difference of energy of the proposed state and the present state, and T is the temperature.

The simulation is run for a long enough time to ensure that the system reaches a metastable, local minimum, of the energy function: the energy fluctuates around a well-defined average value. The time evolution of the system is monitored by observing three different quantities - i) the total energy of the system, E, which undergoes a sharp decrease from the initial high-energy state, and then fluctuates around a much lower value, ii) the average nematic scalar order parameter, S = ⟨ cos 2ϕ⟩, which evolves from ≈ 0 in the initial state to a value close to unity in the low energy state, iii) the average weighted order parameter, Ω = ⟨P _i cos 2ϕ⟩/⟨P⟩, which shows a similar trend, growing from ≈ 0 in the initial state to a large positive value ( ≥ 1 ) in the final state. The time evolution of these three quantities, for a sample run, are plotted in Figure 2.

Figure 2 shows initial and final configurations of a sample run. The initial state (A) shows dipoles that are placed on the sites of the square lattice, with random orientations and polarization magnitudes drawn from a k − Γ distribution. The final state (B), which has a much lower energy, has a majority of the dipoles aligned along the horizontal direction. We also observe that the dipoles are sorted into well-defined domains by their polarization magnitudes.

In Figure 3, we take a closer look at the spatial organization of the order parameters S and Ω in the final, low-energy states. The three rows represent results from simulations with three different values of θ (Eq. 7). Each row represents adifferent representation of the same end state, in order to highlight different measures of the ordering. The leftmost figure in each row shows the two defining properties of every dipole - orientation (arrows) and polarization magnitude (heatmap). We see clear formation of domains of large polarization in all three sets of data. The middle figure in each row shows local values of scalar order parameter S which again organizes the system into domains of nematic order, oriented in different directions. In the rightmost figure of each row, the heatmap shows local values of the weighted nematic order parameter Ω. From the figures it is clear that Ω provides a better characterization of the “order” than the pure nematic order parameter, S, and clearly identifies regions where the nematically ordered domains overlap with domains of large polarization magnitudes. For clarity we have also plotted the polarization magnitudes using contour lines.

These figures clearly demonstrate, one of our primary findings, that domains of largest polarization magnitude overlap with domains of highest nematic order. Although not surprising in retrospect, this aspect of the mechanistic model is missing from pure liquid crystalline models of tissues where the magnitude of P is the same for all cells, and points to the need for characterizing cell-polarity in tissues by both the magnitude of the polarization and its orientation. Our results imply that because of mechanistic interactions between cells, defects in orientational order would likely be localized to regions where the cells are not significantly distorted from an isotropic shape. The consequences of this for tissue development and morphogenesis needs to be explored further. Therefore, further exploration of the mechanisms that regulate cell orientation and the consequences of defects in cell orientation for tissue development and disease is crucial. This could involve using advanced imaging and computational techniques to study the dynamics of cell orientation in tissues, as well as investigating the role of mechanical and biochemical factors in regulating cell orientation. Ultimately, a better understanding of these processes could lead to new therapeutic strategies for treating a wide range of diseases and disorders.

In a separate set of simulations, of Model 2, we keep the dipoles fixed at the same positions as in the initial disordered state. Dipoles do not exchange polarization magnitudes and are allowed to update only their orientation angles. Please see Supplementary Material for results.

The dynamics of cells within solid-like biological tissues is, of course, much more complex than can be captured by Model 1 or 2, discussed in this work. To connect to observations in a much more realistic model of tissues, we present results from simulations of the vertex model, subjected to external shear. Most biological tissues undergoing shape change due to development, growth, morphogenesis, wound healing, etc., Are under external stresses. We have not extended our mechanistic model to explore the effect of external stresses. Below, we compare the order-parameter correlations that develop in a vertex model under shear to results from Model 1.

The Voronoi-based implementation Bi et al. (2016) of the vertex model Farhadifar et al. (2007); Li et al. (2019, 2018); Yan and Bi (2019); Mitchel et al. (2020); Das et al. (2021) is employed to model a 2D cell layer, with the cell centers r _i serving as the degrees of freedom and their Voronoi tessellation dictating the cellular structureBi et al. (2016). The mechanics of the cell layer are described by the energy function Staple et al. (2010). E = ∑ i = 1 N K A A i − A 0 2 + K P P i − P 0 2 (8) , where N is the number of cells, K _A and K _P are the area and perimeter elastic moduli, respectively, and A _i and P _i are the area and perimeter of the ith cell. A ₀ and P ₀ are their corresponding equilibrium values. The origin of the first term in the expression, which is quadratic in the cell areas A _i , can be attributed to the incompressibility of the cell volume. As a result, a 2D area elasticity constant K _A and a preferred area A ₀ are produced, as discussed in Farhadifar et al. (2007); Staple et al. (2010). The second term in the expression, which is quadratic in the cell perimeters P _i , is a result of the contractile nature of the cell cortex and is described by an elastic constant K _P Farhadifar et al. (2007). The target cell perimeter P ₀, which represents the interfacial tension between adjacent cells arising from the competition between cortical tension and adhesion Farhadifar et al. (2007); Staple et al. (2010). Here, we focus on the case where all cells share the same single cell parameters: K _A , K _P , P ₀, A ₀. We also choose A 0 = A ̄ , the mean cell area, which defines the unit of length. For all results presented in this work, we used N = 400 cells.

To study tissue mechanical response in the vertex model, we subject the model tissue to quasistatic simple shear using Lees-Edwards boundary conditions Allen and Tildesley (1989); Huang et al. (2022). We start from a strain-free state, the strain γ is increased in increments of Δγ = 2 × 10⁻³, while cells are subject to an affine deformation given by Δ r i = Δ γ y i x ̂ . At each strain step, the tissue energy is minimized to find a mechanically stable configuration.

We plot the data obtained from vertex model simulations in Figure 4. Figure 4A shows initial state (γ = 0) of the tissue with the cellular aspect ratios indicated by the heatmap and cell orientations indicated by arrows. Figure 4B shows a configuration of the tissue in the quasistatic plastic flow regime at larger strain values. Here the cell shapes are oriented at ∼ 45° to the x-axis due to shear. We observe alignment of cells and formation of domains of large polarization magnitudes parallel to the direction of the external shear. Figures 4C, D shows the corresponding heat maps for the order parameter Ω, and contour lines indicating magnitude of polarization. Here again we observe that the domains of large polarization coincides with domains of large nematic order, similar to the Monte Carlo simulations. This system reaches steady state following a very different mechanism from the Monte Carlo simulations. In the steady state, the cells undergo elongation and plastic failure Huang et al. (2022) multiple times as indicated in the time series plot of the mean polarization magnitude. Still we find end states that are qualitatively similar to those of annealed Monte Carlo runs. The emergence of order in the system is indicated by the time series plots of S and Ω, both of which starts at a value close to 0 in the initial state and plateaus to a large positive value as time progresses.

We are specifically interested in understanding the differences/similarities between (a) correlations, and (b) energy landscapes that result from our mechanistic model and the vertex model. These comparisons are discussed in the context of Figure 5, and Figure 6, respectively. Our perspective is that both model cell shapes in tissues but ours is a purely mechanistic model, and comparison between the two can distinguish between mechanistic and geometrical influences on the spatial organization of cell shapes. Both Monte Carlo and vertex model simulations show the common features of i) formation of ordered domains, ii) formation of domains of large polarization, and iii) overlap of these two types of domains.

In order to quantify these observations we calculate spatial correlations of (a) the polarization magnitudes (C _P ) and (b) the order parameter (C _Ω) for the two models (Figure 5). In the initial states, the correlation of both P and Ω die out very fast, indicating that there is no correlation in these variables to start with. Since we observe formation of anisotropic domains with a long and a short axis in the final states, we calculate correlations in two independent directions - i) parallel to the largest domain and ii) perpedicular to it. Because of the externally imposed shear, the “parallel” direction is always at 45° to the x-axis in case of the vertex model simulation. In the case of Monte Carlo simulations, for each final state, we align the largest domain along the x-axis, for convenience, and determine the average correlations over all final states. In the final states of both the Monte Carlo simulations and the steady states of the vertex model simulations, C _P and C _Ω decay simultaneously indicating a strong correlation between the two quantities. For the Monte Carlo runs we see a significant increase in correlation lengths of both quantities to almost 1/5-th of the system size in perpendicular direction and system size in parallel direction. This sharp contrast between longitudinal and transverse correlations is also evident in Figure 3C depicting nematic order in Armengol-Collado et al. (2022). There it was remarked that these “chain-like” correlations are reminiscent of force-chains in granular media. We want to emphasize that the origin of this type of correlation in the polarization that we observed is the same as that observed in stresses in granular media (Nampoothiri et al., 2020). In Figures 5A, B the heatmap of Ω is analogous to a heatmap of grain-level stresses shown in Figure 3 of Nampoothiri et al. (2020), since Ω is a measure of P _xx . The origin of the longer-ranged correlations of the dipole tensor P _αβ in our model tissues and stress components in jammed granular solids is the imposition of the local constraint of force balance on each cell Bischofs et al. (2004) and grain (Nampoothiri et al., 2020; Nampoothiri et al., 2022). As discussed at length in the context of granular solids, this constraint leads to a Gauss’s law type constraint in the continuum elasticity theory (Nampoothiri et al., 2020; Nampoothiri et al., 2022), leading to a pinch-point singularity in stress-stress correlations that implies much longer-ranged correlations in the longitudinal direction. The Greens function in Eq. 5 encapsulates the force-balance constraint, and is directly responsible for the observed difference between longitudinal and transverse correlations. The visual appearance of “force-chain” like structures is a reflection of these correlations.

As in the Monte Carlo simulations, the correlations observed in the steady-state of the vertex model (Figure 5F) are stronger in the parallel direction than those in the perpendicular direction. But the correlation lengths are systematically smaller in the vertex model. We believe that this is due to a difference in cell dynamics in the two models. Cell dynamics in the vertex model simulation is subject to geometric constraints which prevent “long-range exchanges” of cell aspect ratios (polarization magnitues), allowed in Model 1 of the mechanistic model. These Monte Carlo moves in Model 1 facilitate formation of large ordered domains of polarization magnitudes, spanning the system size. As observed in Figure 5B, the steady states of the vertex model are instead characterized by multiple domains of large polarization magnitudes. This reduces the correlation length in the vertex model simulations compared to those in the Monte Carlo simulations. If instead of Model 1, we analyze the results of Model 2, where we spatially quench all the force dipoles in the Monte Carlo simulation, we observe several small domains, and much shorter correlations lengths (See Supplementary Figure S1 in Supplementary Material). We expect the dynamics in a physical tissue to lie somewhere in between the two extreme scenarios represented by Model 1 and Model 2, because there are additional constraints on how cells can migrate or change their geometrical shape. It is more likely that there is a characteristic length scale over which cells can exchange polarization, which will define a domain size.

Lastly, we analyze the energy landscape explored at low temperatures by looking at the energies of the final states in the simulations. In the Monte Carlo simulations, independent runs achieve different values of the final energy. We observe, however, that if the simulation is run with a fixed set of P but with different spatial realizations in the initial state, they attain final states with energy values that are virtually indistinguishable from each other. This observation indicates that the energy landscape at low temperatures is not sensitive to the initial configurations but is controlled by the particular realization of P, which are drawn from a k − Γ distribution. To quantify this feature, we plot the energy as a function of the average field h _αβ (Eq. 6) in the final states. Each point in Figure 6 represents one configuration, with its y coordinate indicating total energy and x coordinate indicating value of average field h α γ = 1 N ∑ i ∑ i < j G α β , γ β ( r i − r j ) P ( r j ) . We see a strong correlation between energy and the average local field h _xx + h _yy in both Monte Carlo and vertex model simulations. As seen from Eq. 6, the average local field is controlled by the total polarization, ∑ _i P _i . We have checked that replacing average field h _αβ by the total polarization produces plots that capture the same correlation as seen in Figure 6. The vertex model runs also show a strong correlation between final energy and the transverse field h _xy because of the shearing in that direction, but this trend is absent in the Monte Carlo simulations, where the values of h _xy are also much smaller because of the isotropy of the system.

The other aspect of the model that we explore in this study is the energy landscape explored in Model 1. Figure 7A shows the time dependence of energy for ten independent Monte Carlo runs. Each of our simulations is initiated with a set of polarization magnitudes drawn from a k − Γ distribution which is then kept fixed throughout that run. The initial configuration is completely random and hence is a high energy state. As the simulation progresses, the system approaches a low energy configuration. It is evident from the plot that each system approaches a different value of the final energy. The energy of the final state depends on the average value of polarization magnitudes of the dipoles. We performed three different sets of simulations fixing the value of k to be 2.5 but with θ = 1.0, 2.0, 3.0. Since the systems are of finite size, each set of polarization magnitudes drawn from, nominally, the same k − Γ distribution is characterized by k, θ values that are slightly different, leading to different values of the final energy.

If we normalize the final energies from all three sets of runs by the square of the average polarization of the respective run, then they fall in a pretty narrow range of values as shown by the distribution in Figure 7B.

In Figure 8A, we plot the normalized final energies of each system against their respective k, θ values. In all three sets of data, the range of the normalized energy is approximately the same. We divide the data points into three sections by their energy values and plot correlations in two mutually perpendicular directions r ̂ ‖ and r ̂ ⊥ . r ̂ ‖ is defined as the direction in which the largest polarization domain is aligned and correlation C _‖ is calculated by taking pairs of points lying on thin strips along this direction. Similarly for C _⊥ we take pairs of points on thin strips along the perpendicular direction. We observe that both P and Ω have longer-ranged correlations in the parallel direction, as compared to the perpendicular one. In addition, as we go to higher values of absolute energies (lower in the energy landscape) the systems become more strongly correlated in the parallel direction whereas energy has no significant effect on the correlations in the perpendicular direction. Unsurprisingly, the conclusion is that lower energy states are more ordered.