Activity-Induced Fluidization and Arrested Coalescence in Fusion of Cellular Aggregates

At long time scales, tissue spheroids may flow or appear solid depending on their capacity to reorganize their internal structure. Understanding the relationship between intrinsic mechanical properties at the single cell level, and the tissue spheroids dynamics at the long-time scale is key for artificial tissue constructs, which are assembled from multiple tissue spheroids that over time fuse to form coherent structures. The dynamics of this fusion process are frequently analyzed in the framework of liquid theory, wherein the time scale of coalescence of two droplets is governed by its radius, viscosity and surface tension. In this work, we extend this framework to glassy or jammed cell behavior which can be observed in spheroid fusion. Using simulations of an individual-cell based model, we demonstrate how the spheroid fusion process can be steered from liquid to arrested by varying active cell motility and repulsive energy as established by cortical tension. The divergence of visco-elastic relaxation times indicates glassy relaxation near the transition toward arrested coalescence. Finally, we investigate the role of cell growth in spheroid fusion dynamics. We show that the presence of cell division introduces plasticity in the material and thereby increases coalescence during fusion.


INTRODUCTION
A general understanding of the rheological properties of multicellular tissues is important to gain insight into the physics of morphogenetic processes during development. Furthermore, robust models of these materials allow for the design and characterization of generic unit operations, such as aggregation, dispersion, and fusion, which are used for the production of artificial tissues. Given its analogy to the merging of two liquid droplets, the fusion of tissue spheroids has received considerable interest as a model for soft tissue rheology. An analytical expression of the onset of coalescence of two equal viscous droplets under influence of their surface tension was first derived by Frenkel [1] and was improved and extended upon so that the dynamics of complete coalescence could be accurately modeled [2,3]. Furthermore, various extensions to this framework have been proposed to take into account specific properties of multicellular materials, for example differences in spheroid size [4], and the presence of biological processes such as proliferation, differentiation and apoptosis, which may conflict with the assumption of conservation of mass [4,5].
However, the liquid model cannot consistently reproduce in vitro tissue behavior. For instance, Kosheleva et al. showed that fusion dynamics did not correspond to liquid model predictions based on nano-indentation surface tension measurements [6]. Furthermore, arrested fusion has been observed between spheroids treated with Rho-kinase inhibitor affecting acto-myosin contractility [7,8]. Arrested fusion was shown to be controlled by cancer cell activity in tumors [9]. To explain these observations, Oriola et al. recently proposed to model spheroids during fusion as viscoelastic materials instead of simple liquid droplets [10]. Such an approach has already been applied to parallel plate compression [11], micro-pipette aspiration [12], and tissue-detachment/ fracture assays [13]. In general, a tissue may behave like an elastic solid at short time scales and like a liquid at longer time scales [14]. However, they can also appear solid-like at long time scales. This may occur when tissues change from a fluid-like to a solid-like state by undergoing a jamming transition [15]. In the jammed state, individual cells are caged by their neighbors preventing the cells from rearranging, resulting in glassy rheology, whereas in the unjammed state cells are able to move freely, resembling a fluid-like state.
Multiple computational approaches such as phase field models [16], Monte Carlo based models [17][18][19], and individual cellbased models [4,10,[18][19][20][21][22], allow for the characterization of tissue-scale rheological behavior as a function of cell-scale mechanical properties or to simulate fusion of spheroids in more complex geometries. Schötz et al. calibrated an individual cell-based model to better describe the random motion of cells in tissues close to a glass transition. Nevertheless, when investigating fusion with that model, it could still be described by a liquid model. On the other hand, simulations that resulted in a solidification during spheroid fusion have already been mentioned by Kosztin et al. [20]. This observation was attributed to high levels of cell adhesion, but was not further analyzed. Recently, using individual cellbased computational models, it was demonstrated that a divergence of the visco-elastic relaxation time can be observed within a transition region of arrested fusion, indicative of a jammed system [10]. The existence of arrested fusion could thus be traced back to the build up of elastic energy during fusion.
In this work, we derive a simplified analytical expression for the coalescence of visco-elastic tissue spheroids on the basis of elasticity caused by internal structure of the droplet [23,24]. Doings so we obtain similar results to [10]. We demonstrate the applicability of this expression by using an individual-cell based model approach to simulate the fusion process and, using this expression, we are able to compare the relaxation dynamics of fusion to the characteristic timescales involved in the preceding aggregation (or compaction) process in which the individual tissue spheroids have been formed. Finally, we extend this simulation framework to include a morphological model of the cell cycle. In analogy to active motility, cell division and growth may introduce excitation that may induce fluidization, as was observed in other systems [25][26][27][28]. Here, we investigate whether, in spheroid fusion, there is an additional active contribution of cell division beyond a mere correction for the increase in volume, and assess to what extent this contribution may unjam the cellular material [4,5] and recover tissue coalescence.

Individual Cell-Based Model
We follow an individual cell-based model approach as described in detail in Smeets et al. [29,30]. This model has already been used to study cell aggregation and compaction dynamics. In this paper we extend the model by taking into account cell proliferation. All details of the individual cell based model can be found in Supplementary Section S1. In brief, the model is based on a simulation framework introduced by Delile et al. [31], in which cells are simulated as self-propelled particles. Conservative active forces are exchanged between neighboring cells, similar to [10]. The connectivity network is based on a Delaunay triangulation of the cell center coordinates. For the edges of this network, a symmetric central potential is calculated, which is parameterized by adhesion w a and cortical tension w r . Protrusive active forces are responsible for active migration with velocity v t . These forces are calculated in the direction of the cell's polarization which is randomly diffusing with rotational diffusivity D r . Hence, activity from cell motility may be parameterized as D eff : v 2 t /(2D r ). Overdamped equations of motion are integrated to evolve the system over time. Finally, cells are able to grow and multiply based on a cell cycle model, which is explained in Supplementary Section S1.2. The model parameters are listed in tables in Supplementary Tables S1, S2.

Simulation Setup
The simulation pipeline is shown in Figure 1 and follows a sequence of generic unit operations to form a small tissue via fusion. In the first step we simulate the seeding and spontaneous aggregation of two aggregates for 24 h, each in their own microwell, similar to [30]. This guarantees that the shape of each spheroid is consistent with respect to its underlying mechanical properties. It should be noted that, depending on the mechanical properties that govern the aggregation process, this relaxed configuration may have a highly irregular shape. Any (few) remaining cells that did not get incorporated in the main aggregates, are removed from the simulation in analogy to a real fusion experiment. In the second step, the two aggregates are transferred to a larger micro-well and are brought into contact with each other, simulating another 50 h during which the fusion process naturally progresses. In the final step, we extract the spheroid contours and fit two circles to compute the contact angle θ for each simulation, as explained in (Supplementary Section S2). Each realization of a simulated fused spheroid is initialized independently. To calculate the average fusion dynamics, i.e., the average of sin[θ(t)] 2 across all repeats, we take into account that not all spheroids start fusing at the same time, because the initial contact between the spheroids can be weak. To account for this, we define the time at which fusion starts as the last time at which we observe two distinct objects in the extracted contours. For further analysis, the extracted shape measures are shifted in time using this offset.

Visco-Elastic Approximation of the Fusion Process
The fusion dynamics of two equal visco-elastic spheres can be approximated by the differential equation as derived in (Supplementary Section S3). Here, θ is the angle as shown in Figure 2, and is influenced by the surface tension Γ, the radius of the spheroid before fusion a 0 , the apparent viscosity of the tissue η and the shear modulus of the tissue G′. We combine these parameters into two characteristic time constants; τ Γ : a 0 η/Γ is the visco-capillary time, and τ ε : 4η/G′ is the visco-elastic time. The analytical solution of Eq. 1 is FIGURE 2 | Shape evolution during spheroid fusion adapted from [3]. The shape is described by two intersecting circles with equal radius a, and center A and B. O marks the center of the axis system, and depicts the initial contact point. Based on the intersections one defines the neck radius x and the angle θ. During fusion the centers A and B gradually move toward point O, the radius a increases and the doublet length L decreases. a 0 and a f is the radius of the spheroids before, and after complete fusion, respectively. L 0 and L f is the length of the spheroid doublet before, and after complete fusion, respectively.
Frontiers in Physics | www.frontiersin.org June 2021 | Volume 9 | Article 649821 in which the complex number C (Supplementary Section S3) can be obtained from the initial conditions θ(t 0) θ 0 In analogy to previous studies on spheroid fusion, the fusion dynamics are reported as sin 2 (θ) (x/a) 2 with a the radius of the spheroid and x the contact radius which both vary in time [4,10,18,20,21]. Other studies perform the analysis based on (x/a 0 ) 2 [5,6,17,19], but this is less consistent with the underlying theory (see Supplementary Section S4). Eq. 2 is used to calculate sin 2 (θ). When fitting, we use this equation with free variables: τ Γ , τ ε and θ 0 . In order to compare simulated fusion experiments, θ 0 is retained as a free variable since the discrete initial adhesion between the first contacting cell pair will permit a rapid relaxation toward a non-zero initial angle, θ 0 . Based on the fitted values, the predicted equilibrium angle θ eq can be obtained as 3 RESULTS

Arrested Coalescence Dynamics
The average dynamics of simulated tissue spheroid fusion are consistent with the derived visco-elastic material model, expressed in the time evolution of contact angle θ, Eq. 1, as shown in Figure 3B where we varied the cell activity D eff . At sufficiently low D eff , coalescence appears arrested and complete fusion is not attained. Given the correspondence to the viscoelastic model, the parameters τ Γ and τ ε can be interpreted as characteristic timescales of the multi-cellular visco-elastic material. From the fit of τ Γ , τ ε and the instantaneous initial angle θ 0 , we are able to calculate the equilibrium fusion angle, θ eq , using Eq. 4. When varying the cell activity D eff and the repulsive energy w r , two distinct regions can be recognized based on θ eq , Figure 3F. For low activity (D eff ) fusion is arrested, while at higher activities fusion is complete. Higher levels of repulsion w r require more cell activity to fluidize the material and hence to attain complete fusion. Similarly, the characteristic viscocapillary time τ Γ increases when cell activity decreases and when cell repulsion increases, Figure 3E. When reducing cell activity within the fluidized region, the visco-elastic time τ ε gradually increases and displays a divergence near the transition line toward arrested coalescence. This divergence of elastic relaxation time is indicative of an underlying glass transition, as was already pointed out in [10], although based on somewhat different analytical and computational models. The transition from arrested to complete fusion as characterized by θ eq , coincides with the glass transition as characterized by the divergence of τ ε , as shown in Figure 3D. Since this computational model is based on the same individual cell-based framework that was used to simulate the aggregate formation process, we are able to make a direct comparison between dynamics of aggregation and dynamics of tissue spheroid fusion. At sufficient cell activity, the aggregation process is consistent with the dewetting of a liquid film from a surface, as was demonstrated in [30]. Upon decrease in activity, this correspondence is lost. Instead, the aggregate density ρ a follows the dynamics of granular compaction, characterized by stretched exponential relaxation (KWW law): with initial density ρ 0 , final density ρ f , exponent β and characteristic timescale τ agg , see Figure 3A. In this equation the apparent packing density is expressed by ρ a 3V a /(4πR 3 a ) Here, V a is the apparent volume of the spheroid approximated as the sum of cell volumes using an average cell radius. R a is the radius of the spheroid, which is obtained based on the projected area on the bottom of the micro-well, assuming a circular geometry. During compaction, a similar glass transition as observed during the fusion process can be recognized from the divergence of τ agg [30]. For the same underlying cell properties D eff and w r , these two transitions closely align, as demonstrated in Figure 3B, although the transition in the fusion process consistently occurs at higher values of cell activity. We hypothesize that this discrepancy is due to the additional constraints in cell movement imposed by the geometry of the spheroid during the fusion process, compared to the loose network connectivity that characterizes the aggregation phase. Due to the sequential simulation of aggregation and fusion, which mimics the experimental procedure, additional compaction may occur during fusion, particularly at low activity when compaction proceeds slowly. We verified that the effect of this additional compaction on the predicted equilibrium angle θ eq is minor, as shown in Supplementary Figure S7, where we compare fusion of self-assembled aggregates to artificially compacted spheroid structures.

Increase of Coalescence Due to Cell Division
Next, we turn to the role of cell division in the arrest of coalescence. For this, we simulated the fusion dynamics in the presence of cell proliferation, see (Supplementary Section S1.2) Although cell proliferation is not explicitly accounted for in the derivation of our analytical model, Eq. 1 still fits the fusion dynamics of our simulated spheroid fusion in the presence of cell division well, Figure 4A. Figures 4B,C compares the characteristic visco-elastic time constant τ ε and the predicted equilibrium angle θ eq without (k div 2 × 10 − 8 h −1 ) and with (k div 0.06 h −1 ) cell division, when varying cell motility D eff at repulsive energy w r 2.09 min − 1 . This enables us to evaluate the effect of cell growth/division as an additional source of biological excitation compared to active cell motility. The simulated cell division rate corresponds to a cell cycle period of approximately 16.7 h, which is lower than many commonly used cell lines which are used for spheroid fusion. Yet, even at this relatively high division rate, we do not see a strong shift in critical activity beyond which the material appears fluid-like, as indicated by the coincidence of the peak in τ ε for varying D eff with and without cell division ( Figure 4B). However, we do observe that the presence of cell division greatly increases the overall viscoelastic relaxation time, indicating a decrease in the apparent elasticity of the multi-cellular material. Furthermore, cell division markedly increases the equilibrium angle θ eq in the arrested fusion region. Hence, in the simulated configuration, cell division appears to recover coalescence by increasing the plasticity of the tissue. However, it has no strong influence on the location of the fluidization transition.

DISCUSSION
In this work, we derived an expression for the arrested coalescence of tissue spheroid fusion, based on visco-elastic material properties. Simulations of a minimal individual cellbased model of the fusion process showed that this expression is able to describe the transition of liquid-like to arrested coalescence dynamics. Furthermore, a divergence in the viscoelastic relaxation time indicates the presence of jammed or glassy relaxation behavior near the transition toward arrested coalescence. These findings are highly similar to recent work from Oriola et al. [10], hence a brief comparison between these contemporary results is appropriate. First, the analytical expressions for the dynamics of θ (Eq. 1 in [10], compared to Eq. 1) are based on somewhat different assumptions. Our model is an extension of the model of Pokluda [3] by adding an elastic energy term in the equation which was suggested by [23,24]. This has the advantage that in the absence of elasticity, the model of Pokluda is retrieved. The downside of our approach is that this leads to an inconsistency in the strain and strain rate for viscous energy dissipation rate and the elastic energy rate (we note this difference as ε and ε). For the viscous dissipation term, we have, according to Pokluda [3] _ ϵ (1/a)(d/dt)[a(θ)cos(θ)]. On the other hand, the strain in arrested coalescence is suggested to be ε 1 − a(θ)(1 + cos(θ)/(2a 0 )) [23,24]. Therefore the strain rate is _ ε − (1/2a 0 )(d/dt)[a(θ)(1 + cos(θ))]. These expressions are equivalent in magnitude only for the onset of fusion i.e. θ ≈ 0 and a(θ) a 0 . In contrast to our hybrid model, Oriola et al. [10] consistently use the strain rate as suggested for arrested coalescence [23,24]. Furthermore, they introduce elasticity by considering the spheroid as an incompressible Kelvin-Voigt material, instead of adding a separate elastic energy term, hence obtaining a critical value for the Young's modulus for which fusion is inhibited. Moreover, whereas they include a "prestrain" to account for the initial fusion onset, we allow for an instantaneous initial angle θ 0 . Still, both models are parameterized by two essential characteristic timescales: the visco-capillary timescale τ Γ and the visco-elastic timescale τ ε . Secondly, one key difference between the individual cell-based model in this work and the one in [10] is the implementation of the active cell motility. Their model is based on "protrusions" defined on the level of cell-cell bonds, and effectuates persistence by means of a protrusion lifetime per bond. Our model is based on a polarization direction defined for each cell which diffuses over time. Nonetheless, the similarity of the main results underpins that the transition of liquid-like to arrested coalescence is a generic phenomenon that is not dependent on the precise assumptions of the visco-elastic description, nor on The estimated values for the visco-elastic time constant τ ε are based on fitting the average fusion dynamics in sin 2 (θ) using Eqs. 2, 3, the predicted equilibrium angle θ eq calculated by Eq. 4. The error bars for τ ε correspond to two times the standard error obtained from the fits. In these graphs, the error bars are only one time larger than the marker. These graphs show that cell division promotes spheroid fusion, although the system remains arrested at low levels of cell activity.
the implementation details of the individual cell-based model. Additionally, because arrested fusion is the main discrepancy for fusion of the jammed tissue, the simulation of shötz et al. [22] of cells close to the jamming transition can still adequately be described by a liquid model.
The glassy relaxation dynamics during fusion mirror the dynamics of the aggregation process during which the initial tissue spheroids are formed. A direct comparison between these two unit processes shows that there is a clear correspondence between, on the one hand, the transition from granular compaction to liquid dewetting during the aggregation phase, and on the other hand, the transition from arrested coalescence to liquid behavior during the fusion phase. However, this transition occurs for slightly smaller values of cell activity in the case of aggregate formation. Experimental confirmation of this correspondence can be found in studies involving Rho kinase inhibitor, which has been observed to cause arrested dynamics during aggregate formation of human periosteum-derived cells [30], as well as inhibit fusion of spheroids from human mesenchymal stem cells [8] and of embryonal chicken organoids [7]. The role of cell activity in fluidization has been experimentally shown [9]. They related cell activity as obtained from cell tracking to the outcome of fusion. The fusion of low cell activity tumor spheroids appears like a jammed tissue and results in arrested coalescence, while in high cell activity tumor spheroids, fusion is completed. They also noted that the dynamics of this arrested fusion process could not be captured with the classical model. The equation, as derived in this paper will therefore contribute to a better understanding of spheroid fusion in general.
In addition, we considered the effect of cell division on the dynamics of the fusion process. In the framework of arrested coalescence, we showed that the presence of cell division may recover coalescence of fusion. Other studies on the role of cell growth in biological active matter systems, for example in simulations of two-dimensional epithelial tissues [27], or in growing 3D cell aggregates [25], observed an increase of fluidization induced by cell division. More generally, an overview of the role of cell division in tissue rheology and mechanics is provided in [32]. However, in our simulations, the effect of cell division on the location of the fluidization transition was limited, at least for realistic cell division rates compared to the timescale of fusion. Instead, cell division appeared to increase plasticity of the glassy material and thereby improves coalescence in the arrested phase. However, it should be noted that the absence of fluidization as a result of cell division in simulations could be partly due to the minimal representation of cell shape, which introduces artificial energy barriers and thereby overly penalizes neighbor exchanges. This shortcoming could be addressed in more detailed tissue models, such as vertex models [33] or deformable cell models [34,35]. On the other hand, our representation of the cell cycle was limited to cell growth and cell division, while fluidization is often investigated in the presence of apoptosis when the number of cells is stationary or at a homeostatic pressure [28]. We expect that including apoptosis will further fluidize the tissue as this creates random vacancies around caged cells, offering a low energy route for local relaxation. Additional complexities on fusion dynamics could arise as cells enter a pressuredependent dormant state, increasing tissue heterogeneity. For example, the increased pressure in the core of tumor spheroids results in a jammed stage with dormant cells, while cells at the periphery grow and motility is super-diffusive [36].
In practice, technologies that involve the production of artificial tissues frequently incorporate subsequent steps of micro-aggregation and tissue assembly, where the latter often relies on the (partial) fusion of spheroids to create larger tissue constructs [37]. To complicate matters, all these steps are typically accompanied by biological processes such as cell division, production of extracellular matrix, cell differentiation or apoptosis. However, since the physical description of the underlying aggregation and fusion dynamics is highly generic, each of these steps may be parameterized in terms of its characteristic material properties, allowing for the comparison within and between different culture conditions and production formats. As such, continued efforts toward the characterization of structure, rheology and mechanics of these artificial tissues will become indispensable.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the Zenodo repository/repositories and accession number(s) can be found in the article/Supplementary Material.