Melting of two-dimensional biological tissues containing active Ornstein–Uhlenbeck cells

The solid–liquid transition of biological tissues is numerically investigated in the presence of Ornstein–Uhlenbeck noise. We demonstrate that the melting scenario of the system is controlled by three parameters: temperature, the persistence time that controls the nonequilibrium properties of the system, and the target shape index that characterizes the competition between cell–cell adhesion and cortical tension. An increase in the persistence time always causes the system to transition from disordered (liquid state) to ordered (solid state). For stiff cells (small target shape index), on increasing temperature, the system undergoes the first order melting for short persistence time, while it undergoes a continuous solid–hexatic transition followed by a discontinuous hexatic–liquid transition for long persistence time. For soft cells (large target shape index), the melting always occurs via a continuous solid–hexatic transition followed by a discontinuous hexatic–liquid transition and the parameter range where the hexatic phase occurs increases with the persistence time. These behaviors are confirmed by the evolution of the density of topological defects. The phase diagrams of the system are also presented based on three parameters (temperature, the shape index, and the persistence time). Our study may contribute to the understanding of melting in two dimensional systems with many-body interactions and deformable particles.


Introduction
For decades, numerous studies have been done on two-dimensional melting, but there are still some puzzles. In two-dimensional systems with short-range interactions, the melting occurs through a first order solid-liquid transition or through a two-step process with subsequent solid-hexatic and hexatic-liquid transitions [1][2][3][4][5][6]. The two-step melting scenario was presented by Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY), and it suggests that the transition is continuous and defect-mediated [6][7][8][9][10]. The two-step melting via continuous solid-hexatic and hexatic-liquid transitions [11][12][13], or the mixed one, where a continuous solid-hexatic transition followed by a discontinuous hexatic-liquid one [14]. From perfect crystal to the intermediate hexatic phase, the bound dislocation pairs dissociate into free dislocations and the long-range translational order disappears, but quasi-long-range orientational order remains. From the intermediate hexatic phase to liquid, the dislocations unbind into isolated disclinations and the quasi-long-range orientational order is also destroyed. However, the first order melting scenario also exists when the core energy of the free dislocation is less than a critical value [15][16][17][18][19]. Researchers show great interests in the melting scenario of more systems. For example, in the XY model, the melting may also be induced by local disorder and phase coexistence [20]. In the two-dimensional Lennard-Jones systems, people controvert where the boundary of the hexatic phase is [21][22][23]. In hard disks system [24,25], a continuous solid-hexatic transition is followed by a discontinuous hexatic-liquid one.
Recently, researchers have observed experimentally an intermediate mesenchymal-epithelial phase of confluent tissues in various biological processes, including embryonic development, cancer metastasis and wound healing [26,27]. The phase resembles the intermediate hexatic phase in two-dimensional melting [28][29][30][31][32]. By using the cellular Potts model, Durand and Heu demonstrated that the thermally driven melting of soft cellular systems is a continuous transition following most of the predictions of KTHNY theory [31]. Benefit from the successful description of confluent tissue by the Voronoi model, Li and Ciamarra [32] have investigated the role of cell deformability in the two-dimensional melting and found that the soft cells tend to undergo the continuous transition while the stiff ones tend to undergo the first order transition.
Although scientists have done extensive research on two-dimensional melting, the melting scenario for active matter system is still an open problem. Non-Gaussian nature or the persistence of self-propulsion is the fundamental of active matter [33][34][35][36][37][38][39][40][41]. In the previous work [30], Pasupalak and co-workers studied the effect of non-Gaussian properties (active Brownian particles) of active cellular system on the melting scenario and found that the melting transition to occur via two continuous transitions. How the persistence of self-propulsion (another nonequilibrium factor which usually induced by Gaussian colored noises) affects the melting scenario remains unclear. To clarify this problem, we study the solid-liquid transition of two-dimensional biological tissues in the presence of Ornstein-Uhlenbeck noise. The melting behaviors are determined by temperature, the persistence time, and the target shape index. An increase in the persistence time always causes the system to transition from disordered to ordered. For stiff cells, the system melts through a first-order solid-liquid transition with no hexatic phase for short persistence time, while the system undergoes a continuous solid-hexatic transition followed by a discontinuous hexatic-liquid transition for long persistence time. For soft cells, the system undergoes a two-step melting and the hexatic phase range increases with the persistence time.

Model and methods
We consider a confluent tissue consists of N cells propelled by Ornstein-Uhlenbeck noise in a rectangular box with side length ratio L y /L x = 2/ √ 3, under periodic boundary conditions. In contrast to particle-based models, the cells in the Voronoi-based cellular model are confluent and can change their shape to completely fill space, thus their interaction is shaped based. The Voronoi model describes a confluent tissue as a network of polygons. Each cell is characterized by its position r i ≡ (x i , y i ), and cell shape is determined by the resulting Voronoi tessellation. The many-body cellular interaction is captured by the total tissue mechanical energy [30][31][32][42][43][44][45]: where the first term indicates cell incompressibility and resistance of the monolayer to height fluctuations and it is governed by a quadratic in cell area A i . A 0 and K a denotes the preferred cell area and area stiffness moduli, respectively. The term quadratic in cell perimeter P i results from the competition between cortical tension and cell-cell adhesion. P 0 and K p are the preferred cell perimeter and perimeter stiffness moduli, respectively. The effective dimensionless target shape index p 0 = P 0 / √ A 0 is an important parameter that controls the elastic behavior of the cells.
The effective energy of equation (1) leads to a mechanical interaction force on cell i by F i = −∇ i E. In contrast to particle-based models, F i is a multi-body interaction force that cannot be expressed as a sum of pairwise force between cell i and its neighboring cells. Each cell has its self-propelled velocity v i , which means that the system is out of equilibrium. The dynamic of cell i in the overdamped limit follows the Langevin equation: where μ is the mobility.
In order to pinpoint the fundamentals of active matter, scientists tend to study simplified models, in which only the persistence or the non-Gaussian nature is retained. In this paper, we consider the former case (referred to as active Ornstein-Uhlenbeck particles), which has been used to model the collective dynamics of cells [47] as well as for the motion of passive tracers in an active bath of bacteria [48]. In this case, the self-propelled velocity v i is described by the Ornstein-Uhlenbeck processes [39,46,49]: where ξ i (t) are Gaussian white noise with zero mean and unit variance. D = μk B T controls the noise amplitude and the persistence time τ controls the nonequilibrium properties of the system. k B is the Boltzman constant and T is temperature. For vanishing τ , v i reduces to the Gaussian white noises and the system corresponds to the one in thermal equilibrium (the case similar to passive Brownian dynamics). For finite τ , the velocities form a set of zero-mean Gaussian colored noises with correlation v ix (t) v jy (0) = D τ e −|t|/τ δ ij δ xy [46,49]. We point out the relevance of the ratio D/τ , i.e. the variance of v i , whose square root gives the typical average value of the self-propelled velocity, so D/τ is the scale of the self-propulsion speed. Since the temporal correlations of the Ornstein-Uhlenbeck processes are not matched by similar correlations for the drag, the system does not satisfy the standard generalization of the Stokes-Einstein relation to systems with memory, so the system is out of thermal equilibrium. The active Ornstein-Uhlenbeck particle is the simplest Gaussian model which displays at long times the same average and two-time activity-activity correlation function as active Brownian particles [50]. In order to focus on the effect of the persistence time on melting, we fix the noise amplitude D, instead of the self-propulsion speed D/τ , and vary the persistence time τ .
To get more insights into the structural change during the melting process, we focus on the two order parameters which characterize the orientational and translational symmetries of the system, respectively. The global orientational order parameter is defined as where ψ 6 (r i ) = 1 n n m=1 e i6θ im is the local orientational order [1,[7][8][9][10]. n is the number of Voronoi neighbors of cell i. θ im stands for the angle of the bond between centers of cell i and m relative to a fixed reference axis. The global translational order parameter is given by where ψ t (r i ) = e i G·r i is the local translational order [1,[7][8][9][10]. G is one of the first Bragg peaks, identified by the static structure factor [23]. Note that in all calculations involving positional order and correlation, we use the G from the actual positions of the Bragg peaks, instead of from the expected peak-positions of a perfect triangular lattice. The melting transition is accompanied by the destruction of the translational order and the orientational order. It is a vital evidence for hexatic phase that the system has short-range translational order but quasi-long-range orientational order. To further characterize the transition, the spatial decay of the correlation function of the local orientational order parameter is focused within the KTHNY theory [6,7]. When the system in solid phase, the decay exponent η → 0 and C 6 (r) approaches a constant. A hexatic phase is detected by 0 < η 1 /4 and C 6 (r) decays algebraically [51]. In the liquid phase, η > 1 and C 6 (r) decays exponentially.

Results and discussion
Equations (2) and (3) are numerically integrated using a stochastic Runge-Kutta algorithm. The integration time step was chosen to be 0.01 and the total integration time was more than 10 5 (this time is sufficient to ensure that the system can reach a steady state). The model can be nondimensionalized by introducing the characteristic length √ A 0 and time 1/(μK a A 0 ). Unless otherwise stated, we set k B = 1 and K p = 1. For all simulation runs, we start with the same ordered hexagonal tiling (honeycomb) and wait for equilibration before recording data. We explore the melting behaviors by varying the reduced temperature T = k B T K a A 2 0 , the persistence time τ , and the target shape index p 0 . In order to check for the finite seize effects on the melting, in figure 1(a), we have plotted the global order parameter |Ψ 6 | as a function of the reduced temperature T * for different N. It is found that the same global order parameter |Ψ 6 | is obtained irrespective of the system size N and the finite-size effects can be negligible. Therefore, we only consider the case of N = 3960 in the follows. Figure 1(b) depicts the temperature dependence of |Ψ t | and |Ψ 6 | for tiff cells (p 0 = 3.0). Both |Ψ t | and |Ψ 6 | decay with the increase of T and finally tend to zero, which means that the translational and orientational order of the system disappear. When the persistence time tends to zero (e. g., τ = 10 −2 ), the Ornstein-Uhlenbeck noises reduce to the Gaussian white noises. In this case, both |Ψ t | and |Ψ 6 | drop abruptly to zero at the same temperature, which indicates that the system undergoes the first order transition. For the finite persistence time (e. g., τ = 1), the temperature of phase transition increased obviously and |Ψ t | falls slightly before |Ψ 6 |. It suggests that there may be an intermediate hexatic phase with quasi-long-range orientational order and short-range translational order, which agrees with the KTHNY scenario. Therefore, for stiff cells, the presence of the persistence would cause the system to change from a first order melting to the continuous melting.
To further characterize the type of phase transition, we calculate the decay exponent η of C 6 (r) in figure 2(a). When τ = 10 −2 , η rises abruptly and C 6 (r) goes suddenly from being no decay to decaying exponentially, suggesting the first order melting of the system. However, when τ = 1, η increases gently and there exists a parameter regime where C 6 (r) goes through an algebraical decay with 0 < η 1 /4, which indicates the existence of an intermediate hexatic phase. The decay of C 6 (r) is illustrated in figure 2(b) for τ = 10 −2 and figure 2(c) for τ = 1.0. It proves that the melting of the biological tissue in the absence of activity is the first order phase transition, which agrees with the results in the previous work [32]. However, the two-step melting with an intermediate hexatic phase occurs in the system in the presence of Ornstein-Uhlenbeck noise. We can explain the above phase transition behavior via topological defects in the Voronoi partition. According to the KTHNY scenario [1,[7][8][9][10], the two consecutive transitions are mediated by topological defects. The solid-hexatic transition is induced by dissociation of bound dislocation pairs (5-7-5-7) into free dislocations (5-7), while the hexatic-liquid transition is driven by the unbinding of dislocations into disclinations (5 or 7). Figure 2(d) displays the evolution of the number of defects (including dislocations and disclinations) with T . When τ = 1, the population of dislocation (5-7) soars slightly before than disclination (5 or 7). If there are dislocations but not disclinations in the system, it means the presence of hexatic phase. However, both dislocation and disclination rise simultaneously at τ = 10 −2 .
We further confirm whether the solid-hexatic and the hexatic-liquid transition are continuous at large τ by evaluating the probability distribution of local number density P ρ Ω i . The local number density is where N Ω i is the number of cells contained in the tiny domain Ω i of cell i, and A Ω i is the area of the domain. As shown in figure 3, it is a unimodal distribution for the case near solid-hexatic transition (T = 0.7) while it is a bimodal function for the case near hexatic-liquid transition (T = 0.85). Therefore, the system undergoes a continuous solid-hexatic transition followed by a discontinuous hexatic-liquid transition for long persistence time.
We have presented the typical configurations for different phases in figure 4. There are no defects in the perfect hexagonal lattice (shown in figure 4(a)). In the deformed solid phase (shown in figure 4(b)), defects are mainly bounded dislocation pairs (5-7-5-7 quartets), where the orientational order is long ranged and the translational order is quasi-long ranged. In the hexatic phase (shown in figure 4(c)), there are isolated dislocation (5-7 pairs) and no isolated disinclination, where the long-range translational order disappears, but quasi-long-range orientational order remains. In the liquid phase (shown in figure 4(d)), defects are both isolated disclinations and isolated dislocations, where both the orientational order and the translational order are short ranged. Now we discuss the role of the persistence time on the melting of the system. We consider the case of T = 0.1, where the system shows the liquid phase in thermal equilibrium system. The variation of |Ψ t | and  |Ψ 6 | with τ is shown in figure 5 at T = 0.1 for different p 0 . On increasing τ from 10 −2 , both |Ψ t | and |Ψ 6 | increases from zero and finally tends to 1. In a word, an increase in the noise persistence always causes the system to transition from disordered to ordered. This is can be explained as follows. When the persistence time τ → 0, the Ornstein-Uhlenbeck noises reduce to the Gaussian white noise and the system corresponds to the one in thermal equilibrium system, so the system is in the liquid state. When the persistence time τ → ∞, the typical self-propulsion speed D/τ tends to zero, the system is in the solid state (the initial configuration is chosen to be a perfect hexagonal lattice). Therefore, we can obtain the liquid-solid  transition by changing the noise persistence. In addition, the increase of |Ψ 6 | prior to |Ψ t | indicates that there exists a hexatic phase in the transition. Figure 6(a) illustrates the decay exponent η as a function of τ . As τ increases, both curves of η decline monotonically where contains the regions of algebraical decay of C 6 (r) with 0 < η 1 /4. The existence of the hexatic phase can be confirmed. Figures 6(b) and (c) display the details of decay tendency of C 6 (r) for different state at different p 0 , including the case η ≈ 1 4 . Another prove of the hexatic phase is the dependence of number of defects on persistence time τ described in figure 6(d). The dislocations disappear later than the disclinations, which indicates the appearance of the hexatic phase. What's more, the liquid-solid transition τ for p 0 = 3.3 is lower than that for p 0 = 3.5. Now we will give the physical interpretation of the above results. When τ increasing, the self-propulsion velocity of the system decreases and the puny energy is not enough to produce topological defects at a sufficiently large τ . Besides, a smaller p 0 means the system tends to solid-like state [44] so that the critical τ of liquid-solid transition moves forward.  Investigating the melting by increasing T and τ for stiff cells (p 0 = 3.0), we obtain the phase diagram in figure 7. The boundary between the solid and the hexatic phase (the red line) and the one between the hexatic phase and liquid (the blue line) are obtained by analyzing the decay of C 6 (r) and further confirmed by the number of defects. At small τ values the system melts through a first-order solid-liquid transition with no hexatic phase. However, a two-step melting with a continuous solid-hexatic and a consecutive first-order hexatic-liquid transition appears for large persistence time. The range of T values where the hexatic phase occurs expands as the persistence time increases, moves toward higher T .   Figure 8 displays the dependence of |Ψ t | and |Ψ 6 | on the shape index p 0 for different T at τ = 1. When p 0 increasing from 3.0, both |Ψ t | and |Ψ 6 | drop to zero, but the decay of |Ψ t | slightly before |Ψ 6 |. When T = 0.2, the critical shape index is smaller than the value at T = 0.1. According to the variation of |Ψ t | and |Ψ 6 |, we can infer the transition is a two-step melting.
The confirmation of the intermediate hexatic phase in based on the regions of algebraical decay of C 6 (r) (shown in figure 9(a)). It is obvious that the decay exponent η ≈ 1 4 when p 0 3.39 (T = 0.2) and When p 0 increasing, cellular viscoelasticity dominates the competition with tension of cells, so the mobility of cells is increased and the topological defects will be formed more readily. Moreover, a higher T means that the system has a higher average kinetic energy, and the particles easily get out of the cage of neighbors and leave a topological defect. Therefore, the critical temperature of transition moves forward under a larger T .
To study in more detail the dependence of the melting scenario on the shape index p 0 and the persistence time τ , we plot the phase diagram as a function of τ and p 0 at T = 0.1 in figure 10. Both boundaries are determined by analyzing the decay of C 6 (r) and further confirmed by the number of defects. The system undergoes the first order melting when p 0 and τ are small, while a two-step melting appears when p 0 and τ are large. The range of p 0 (τ ) values where the hexatic phase occurs expands as p 0 (τ ) increases, and moves toward larger p 0 (τ ) values.
Note that in our study we fixed the noise intensity D, instead of the scale of the self-propulsion speed D/τ , so our system is frozen in its initial state for large τ and the state of the system depends on its initial state when τ → ∞. For example, the system tends to jam in a disordered state if it starts from a liquid like configuration. A similar study performed in a system of hard disks [52] shows that the system transients from the solid to the liquid phase in the limit τ → ∞, which is different from that in our system. This can be explained as follows. In this work [52], the mean drift velocity v is fixed, when τ → ∞, the persistence distance l = vτ is large enough to form topological defects and melting appears, so the system transients from the solid to the liquid phase. However, in our system, when τ → ∞, the scale of the self-propulsion speed D/τ tends to zero, so the system is frozen in its initial state (crystalline).

Concluding remarks
In summary, we have presented a detailed study of the solid-liquid transition of the confluent tissue consisting of active Ornstein-Uhlenbeck cells in rectangular box with periodic conditions. Each cell is represented by the polygons obtained from Voronoi tessellations and driven by the Ornstein-Uhlenbeck processes. We found that the solid-liquid transition of the system is determined by the three parameters: temperature, the persistence time and the target shape index. The system is frozen in its initial state (crystalline) when the persistence time is very long. For stiff cells, the system undergoes the first order melting for short persistence time, while it undergoes a two-step continuous melting for long persistence time. For soft cells, the system always undergoes a two-step continuous melting and the parameter range where the hexatic phase occurs increases with the persistence time. These behaviors are confirmed by the evolution of the density of topological defects. To study in more detail the dependence of the melting scenario the system parameters (p 0 , τ , T ), we have presented the phase diagrams in the τ -T and p 0 -τ plane.
Active Ornstein-Uhlenbeck particles are usually used to model the collective dynamics of cells [47] as well as for the motion of passive tracers in an active bath of bacteria [48]. The persistence time corresponds to the correlation time of the Ornstein-Uhlenbeck process of cells, which has been presented foremost by Fürth [53] and was shown to be applicable to fibroblast crawling by Gail and Boone [54]. The stochastic motion of epithelial cells, such as the Madin-Darby canine kidney cells [47,55], can be described by the Ornstein-Uhlenbeck process. The correlation time of the Ornstein-Uhlenbeck process can be controlled in experiments. For example, the persistence time can be controlled by the blebbistatin in the cytoplasm of living cell experiments [56]. Our study may contributes to the understanding of solid-liquid transitions in two dimensional systems with many-body interactions and deformable particles and these transitions are biologically relevant to understanding of embryonic development, wound healing, and cancer.