Arrested Coalescence: A Tool to Explore Tissue Rheology

Abstract

Tissue spheroids are self-organised 3D cellular aggregates that serve as a versatile platform in tissue engineering. While numerous high-throughput methods exist to characterise the cellular function of tissue spheroids, equivalent techniques for the mechanical characterisation are still lacking. In this review, we focus on tissue fusion— a simple, fast, and inexpensive method to characterise the rheology of tissue spheroids. We begin by discussing the implications of tissue rheology in development and disease, followed by a detailed explanation of how the phenomenon of arrested coalescence can be used to explore the rheology of tissue spheroids. Finally, we present different theoretical models that, when combined with experimental data, allow us to extract rheological information.

1. Introduction

Tissue spheroids are small, self-organised 3D cellular aggregates that mimic the architecture and function of native tissues. Unlike traditional two-dimensional cell cultures, spheroids provide a more physiologically relevant environment. This 3D arrangement is crucial for studying complex biological processes, as it closely resembles the natural organization found in living organisms. Their small size (∼100 $\mathsf{\mu }$m) compared to bulk 3D tissues (∼10 mm), makes them attractive for scalable tissue fabrication in a wide range of fields, such as bioprinting [1], drug discovery [2,3], or regenerative medicine [4]. In addition, by further culturing these 3D tissues under appropriate conditions, they can be used in organoid generation [5,6] as well as disease modelling [7] (Figure 1). Multiple high-throughput omics techniques exist to characterise the molecular state [8] as well as the mechanical state of cells [9,10]. For the mechanical characterisation of 3D microtissues [11,12,13] a variety of methods exist, such as tissue surface tensiometry [14,15], atomic force microscopy [16], micropipette aspiration [17,18,19], oil droplets [20,21], microgel elastic probes [22,23], or magnetic beads [24]. Nevertheless, the majority of the previous techniques are low-throughput, and the high-throughput mechanical characterisation of microtissues remains a challenge. In this review, we provide a detailed description of how tissue coalescence can be used as a high-throughput, non-contact method for the rheological characterisation of tissue spheroids. The review is structured as follows: we first motivate the study of tissue rheology by discussing its implications on development and disease, we next provide a detailed explanation of the tissue coalescence method, and finally discuss how computer simulations of tissue fusion can help to bridge the cellular and tissue scales.

2. The Role of Tissue Rheology in Development and Disease

In recent years, several studies have brought about an understanding of the material properties of living tissues by using an analogy to passive, inert materials [25]. Such an association can be made by regarding cellular collectives as disordered materials that can transition from fluid to solid states [26,27,28]. In fluid-like states, cells are able to freely rearrange upon the application of external forces by remodelling their architecture. Conversely, in solid-like tissues, cell movement is restricted, enabling the material to bear elastic stresses. Both behaviours are crucial for tissue form and function [29]. Experimentally, the rheological properties of tissues are measured by monitoring the strain or deformation of the tissue in response to force application and removal, mainly through creep and recovery assays [11,30]. By choosing an appropriate constitutive model describing the deformation response over time, one can estimate tissue-scale parameters such as the effective elastic modulus, viscosity, and yield stress [11,30,31]. One of the first pieces of evidence of fluid-to-solid phase transitions in living tissues was found in epithelial monolayers in vitro. In this case, a proliferation-driven increase in cell density beyond a critical threshold was found to arrest cell movement, leading to glass-like dynamics [32]. Subsequent studies have demonstrated that cell-cell interactions [33,34,35,36,37,38,39,40,41], supracellular stresses [42], as well as active fluctuations [43,44,45] can also drive phase transitions in tissues. The previous work has motivated the representation of tissue states in a jamming phase diagram [25,26,31]. Notably, fluid-to-solid phase transitions in tissues are not only interesting from a purely biophysics perspective, but they have been shown to play a key role in development and disease.

2.1. Tissue Rheology in Development

Animal shapes are highly diverse with organs of unique morphologies and functions. Such complex multicellular structures originate during embryonic development through morphogenetic signals that alter the cellular dynamics, ultimately patterning and molding the different tissues in the embryo. Recently, several studies have reported gradients in the mechanical properties of tissues [16,18,24,31]. In particular, these gradients have been proposed as one of the main ingredients to mold tissues [26,30] in conjunction with localized cell proliferation and active stresses [46,47]. A well-studied model organism where tissue rheology has been shown to play a key role in early development is zebrafish [18,31,48]. During epiboly, the blastoderm undergoes fluid-like spreading driven by cell-cell contact disassembly during mitotic cell rounding [48]. Further studies have demonstrated that the blastoderm undergoes a rigidity phase transition due to a small reduction in adhesion-dependent cell connectivity below a critical value [18]. In a later developmental stage during tailbud formation, an N-cadherin-dependent gradient in yield stress has been shown to drive body axis elongation by fluidising the posterior part of the tailbud [31]. The previous evidences suggest conserved physical mechanisms of morphogenesis across different organisms, with tissue rheology playing a crucial role [47].

2.2. Tissue Rheology in Disease

Tissue spheroids resemble small micro-tumours and hence they are a valuable tool to study cancer. Historically, cancer detection methods rely on morphological and molecular screening of the tissue sample in a procedure known as histopathology [49,50]. Despite significant advances in this direction, new early detection techniques are needed, and hence recently, attention has also been drawn to the mechanical properties of cancerous tissues [51]. So far, it has been observed that the tumour stroma undergoes abnormal changes, including increased stiffness of the extracellular matrix and the accumulation of stress gradients within its mass [52,53]. By doing so, their invasive and metastatic potential can increase [54,55], especially in colon cancer-derived cell lines, as they show high sensitivity to substrate stiffness and can initiate metastasis-related phenotypes [56,57]. Rheology-based diagnostics have been confirmed feasible by successfully discriminating healthy from cancerous tissues using atomic force microscopy or shear rheometry [51]. Moreover, in contrast to healthy tissues, cancerous tissues have been observed to become significantly more dissipative at increased shear strain, that is, more viscous [51]. Recent studies in brain tumours also suggest the potential of this feature as a new marker in tissue pathology [58,59,60,61]. Unfortunately, rheology is still difficult to assess due to the compromising accessibility to these tissues in the human body and has not been proven to predict the stage of the disease yet, which would be mandatory in clinical oncology [51]. Finally, the previous studies are mainly based on macroscopic rheometry on millimetre-sized tissue samples [51,58]. A promising alternative approach is to generate micro-tumours using patient-derived cells in vitro (see Figure 1) [62], which could improve the throughput of rheological measurements, for example, by means of microfluidic methods [19].

3. Arrested Coalescence to Explore Tissue Rheology In Vitro

Tissue coalescence is a well-known process that has been used in vitro for a variety of purposes in the fields of developmental biology, tissue engineering, and biophysics. Some examples are tissue maturation [63,64], 3D tumour invasion [65,66,67], bioprinting [68,69,70], and tissue rheology [67,71,72,73]. Here we focus on the last application, that is, to obtain rheological information from the coalescence of tissue spheroids. Indeed, tissue fusion is arguably one of the simplest methods to obtain mechanical information from tissues (see

Table 1. Comparison of tissue fusion (TF) to common techniques for the mechanical characterisation of 3D microtissues. Abbreviations: Tissue Surface Tensiometry (TST), Micropipette Aspiration (MA), Droplet-based Sensors (DS), Microgel Elastic Probes ($\mu$GEP), and Magnetic Beads (MB). Despite the fact that MA is usually low-throughput, in a recent study, the technique has been adapted to high-throughput [19].

Method	Measure		Throughput		References
	Absolute	Relative	High	Low
TST	×			×	[14,15,76]
MA	×		×	×	[17,18,19,77,78]
DS/$\mu$GEP	×			×	[20,21,22,23]
TF		×	×		[71,72,74,79]
AFM	×			×	[16,51]
MB	×			×	[24]

). Unlike other tissue rheology methods such as droplet actuation [21] or micropipette aspiration [17], the rheological timescale is not imposed, but rather it is controlled by the system itself through cell-cell interactions. Additional advantages of the method are the fact that there is no need for a calibrated probe, it is a non-contact method, and it can be high-throughput [12,26,74,75]. The main limitation of the method is the fact that only relative measures can be obtained, i.e., ratios between mechanical parameters [26,75].

Apart from the coalescence timescale, which provides information about the fluidity of the tissue, the steady state shape of the fused assembly also provides valuable mechanical information related to the solidity of the tissue. Indeed, incomplete fusion of biological tissues has been observed in human mesenchymal stem cells [80], human mamma and cervix carcinoma spheroids [66], as well as in mouse embryonic stem cell aggregates [79]. In the soft matter field, such a process is commonly referred to as arrested coalescence or partial coalescence [81,82,83]. The stable anisotropic shapes it can produce have been exploited extensively to produce emulsions in a wide range of industries like food, cosmetics, petroleum, and pharmaceutical formulations [84,85,86,87]. Interestingly, this phenomenon has also been observed in various active matter systems such as ant [88] or bacterial [89] aggregate colonies. Here we propose to exploit this phenomenon to mechanically characterise tissue spheroids in vitro in a high-throughput manner. Next we describe how this can be achieved by combining high-throughput imaging together with mathematical modelling.

High-content imaging systems enable the acquisition of large image datasets in an automated manner (see Figure 2). They are commonly used for drug screening purposes [90] or phenotypic characterisation [91]. Recently, some attempts have been made to study tissue coalescence in these high-content systems [65,79]. Next, we specify the requirements and typical workflow for the formation and subsequent analysis of the fusion of 3D cellular aggregates [66,72,79]. Common cell culture techniques to generate spheroids for tissue fusion are pellet culture [72], static suspension cultures (also known as the liquid overlay technique) [92], the hanging drop method [93], or cell aggregation using U-bottom multiwell plates [66,79,94]. Fusion dynamics is critically affected by the size of the aggregates [69,79], thus this is a critical parameter to control in the experiments. The size distribution is typically small for the cell pellet, hanging drop, and cell aggregation methods, while it can be very broad for the static suspension culture [95]. Quantitative analysis of tissue fusion dynamics typically requires image acquisition every ∼10 min for a timespan of several hours or days [66,72,79]. Low-magnification brightfield images are sufficient for image analysis. To quantify the dynamics of tissue fusion, one typically measures the fusion contact angle [66,69,96], neck radius [72,97], or the end-to-end length of the assembly [79]. Typical cell aggregate diameters range from ∼100 to 500 $\mathsf{\mu }$m. Below, a diameter of ∼100 $\mathsf{\mu }$m, it was shown that the continuum limit breaks down [79]. Statistical averaging of fusion events is of crucial importance to reduce the variability in spheroid formation and fusion; thus, multiwell systems are ideal for this purpose. Finally, given the amount of fusion events to analyze (∼100 events for multiwell plates), it is essential to optimize the segmentation pipeline, for example, using software such as MOrgAna [98].

4. Theoretical Modelling of Tissue Fusion

In order to extract rheological information from the dynamics of tissue fusion, a theoretical model describing the mechanics of the aggregate is needed (see Figure 2). Continuum models of tissue fusion are a suitable choice to obtain rheological information with a reduced number of effective tissue-scale parameters [69,79]. Alternatively, agent-based models allow a connection between cell behaviour and tissue-scale properties [79,99], usually at the cost of a larger number of microscopic parameters. In the following sections, we discuss both modelling approaches in detail, highlighting their strengths and weaknesses.

4.1. Continuum Modeling of Tissue Fusion

From a fluid mechanics perspective, cellular aggregates can be regarded as complex materials that generally behave as viscoelastic drops [66,69,71,72,79,97,100,101,102]. Intuitively, one can understand how tissues fuse by considering the interplay of tissue surface tension ($\gamma$), tissue viscosity ($\eta$), and tissue elasticity ($\mu$). In this very simplified picture, tissue fusion is driven by surface tension, dampened by viscosity, and resisted by elasticity [79,81,82,83,97,102]. The typical fusion speed at the onset of tissue fusion is given by the viscocapillary velocity $v_c=\gamma /\eta$ [69,71], while the degree of fusion is characterised by the so-called elastocapillary length ${\ell }_e=\gamma /\mu$ [81,82,103]. Despite the elastocapillary length ${\ell }_e$ being below the atomic scale in ordinary hard solid materials, this lengthscale becomes non-negligible in soft materials such as elastomers, gels, or biological tissues [103,104]. In particular, in the case of tissues with typical ranges $\gamma \sim$ 1–10 mN/m [12,14,105,106] and $\mu \sim$ 10–10³ Pa [12,69,79,107], we obtain ${\ell }_e$ values ranging from the $\mathsf{\mu }$m to the mm scale, comparable to the typical size of tissues. In other words, at the multicellular scale, surface tension forces are comparable to elastic forces and thus play a key role in deforming tissues.

The study of how two contiguous cell aggregates fuse has received notable attention from the biophysics as well as tissue engineering communities in the last decades [69,71,72,108], mainly due to the simplicity of its geometry. Experimental studies show that the shape of the system can be well approximated to that of two contiguous spherical caps [69,71]. This fact allows us to simplify the three-dimensional problem to studying the dynamics of the contact angle $\theta$ between the two aggregates (see Figure 2). To derive the dynamics of $\theta$, studies have followed the approach initially proposed by Frenkel & Eshelby [109,110] and modified by Bellehumeur [111] to study the sintering of highly viscous molten drops. Those studies simplify the flow field inside the droplet assembly by considering a homogeneous biaxial extensional flow during fusion [69]. Next, we provide a detailed derivation of the dynamics of the contact angle $\theta$ by following a minimisation approach of a Rayleighian functional. This method leads to equivalent results as in previous studies [69,79] and provides a systematic way of deriving the dynamics of $\theta$.

We will follow a minimisation approach based on Onsager’s variational principle [112,113] in a terminology introduced by Doi [114,115]. We define a function called the Rayleighian, which reads:

$$\mathcal{R}=\mathcal{D}+\dot{\mathcal{F}}$$

(1)

where $\mathcal{D}$ is the dissipation potential and $\mathcal{F}$ is the free energy of the fused assembly. Let us consider the stress tensor of the material ${\sigma }_{\alpha \beta }={\sigma }_{\alpha \beta }^{\mathrm{d}}+{\sigma }_{\alpha \beta }^{\mathrm{el}}$, where ${\sigma }_{\alpha \beta }^{\mathrm{d}}$ and ${\sigma }_{\alpha \beta }^{\mathrm{el}}$ are the dissipative and elastic parts of the stress tensor, respectively. The dissipative contribution will be expressed as:

$$\mathcal{D}={\displaystyle \frac{1}{2}}{\int }_{\Omega }{\sigma }_{\alpha \beta }^{\mathrm{d}}{\dot{\epsilon }}_{\alpha \beta }dV$$

(2)

where ${\epsilon }_{\alpha \beta }={\displaystyle \frac{1}{2}}\left({\partial }_{\alpha }{u}_{\beta }+{\partial }_{\beta }{u}_{\alpha }\right)$ is the symmetric strain tensor, $u_{\alpha }$ is the displacement field, and $\Omega$ is the integration domain of the assembly. The free energy of the assembly will be given by a combination of the surface tension and elastic energies:

$$\mathcal{F}={\displaystyle \frac{1}{2}}{\int }_{\Omega }{\sigma }_{\alpha \beta }^{\mathrm{el}}{\epsilon }_{\alpha \beta }dV+{\int }_{\partial \Omega }\gamma dS$$

(3)

where $\partial \Omega$ is the surface integration domain of the assembly.

To illustrate the method, we will consider the case of a Kelvin-Voigt-type viscoelastic solid such that ${\sigma }_{\alpha \beta }^{\mathrm{d}}=2\eta {\dot{\epsilon }}_{\alpha \beta }$ and ${\sigma }_{\alpha \beta }^{\mathrm{el}}=2\mu {\epsilon }_{\alpha \beta }$. This choice was shown to successfully describe the fusion dynamics of mouse embryonic stem cell aggregates [79]. Notice that the current method is not restricted to linear constitutive equations. For example, one can derive the equations for the case of an upper Maxwell convected model at steady state as in the work of Bellehumeur [111]. We will further assume that the system is incompressible (i.e., ${\dot{\epsilon }}_{\gamma \gamma }=0$) and thus the volume of the assembly is conserved during fusion. In order to simplify the problem, we make the assumption of a homogeneous biaxial extensional flow during fusion [69,111], which implies ${\epsilon }_{\alpha \beta }\simeq ϵ\left(\theta \right)\mathrm{diag}(1/2,1/2,-1)$, with fusion occurring along the z-axis, where we define the strain in the system as $ϵ\left(\theta \right)=1-{\displaystyle \frac{R\left(\theta \right)}{2R_0}}(1+cos\theta )$. The strain rate $\dot{ϵ}$ can be easily computed and reads:

$$\dot{ϵ}\left(\theta \right)={\displaystyle \frac{R\left(\theta \right)}{2R_0}}\left({\displaystyle \frac{sin\theta }{2-cos\theta }}\right)\dot{\theta }$$

(4)

where $R\left(\theta \right)=R_0{2}^{2/3}{(1+cos\theta )}^{-2/3}{(2-cos\theta )}^{-1/3}$. The latter expression is obtained by imposing volume conservation during fusion [72,79,116]. Using the previous definitions, we find that the dissipative contribution reads:

$$\mathcal{D}=4\pi R_0^3\eta {\dot{ϵ}}^2\left(\theta \right)$$

(5)

where we used that the total volume of the assembly is ${\displaystyle \frac{8}{3}}\pi R_0^3$. Similarly, one can compute ${\dot{\mathcal{F}}}_e$, which reads:

$${\dot{\mathcal{F}}}_e=8\pi R_0^3\mu ϵ\left(\theta \right)\dot{ϵ}\left(\theta \right)$$

(6)

Finally, considering that the surface of the assembly can be expressed as $S\left(\theta \right)=4\pi R^2\left(\theta \right)(1+cos\theta )$ [72,79,116], the term $\dot{{\mathcal{F}}_s}$ reads:

$${\dot{\mathcal{F}}}_s=-2\gamma \left({\displaystyle \frac{\pi R^2\left(\theta \right)sin2\theta }{2-cos\theta }}\right)\dot{\theta }$$

(7)

In our problem, $\theta$ plays the role of the generalised coordinate within the Onsager variational framework, and thus we can obtain the dynamics of $\theta$ by minimising the Rayleighian $\mathcal{R}$ with respect to $\dot{\theta }$, i.e., $\partial \mathcal{R}/\partial \dot{\theta }=0$ [114,115]. The contact angle dynamics then reads:

$$\dot{\theta }={\displaystyle \frac{(2-cos\theta )}{\tau sin\theta }}\left[{\displaystyle \frac{1}{\mathrm{Ec}}}\left(1+cos\theta -{\displaystyle \frac{2R_0}{R\left(\theta \right)}}\right)+2cos\theta \right]$$

(8)

where $\mathrm{Ec}={\ell }_e/R_0$ is an elastocapillary number [103,104], and $\tau =R_0/v_c$ is the viscous relaxation rate of fusion [69]. It is worth remarking that the previous parameters are size-dependent [79]. Equation (8) can be shown to be equivalent to the one obtained in Ref. [79]. It is instructive to consider the approximation $R\left(\theta \right)\simeq R_0$ in Equation (8), which is reasonable in many cases due to the small variation in radius during fusion when the cell proliferation rate is much slower than the rate of fusion [69,79]. In this limit, we find that at steady state ($\dot{\theta }=0$), the equilibrium angle of fusion takes a simple expression:

$${\theta }_{\mathrm{eq}}\simeq \mathrm{arcsec}\left(1+2\mathrm{Ec}\right)$$

(9)

In the viscous limit ($\mathrm{Ec}\to \infty$, i.e., $\mu =0$) the equilibrium angle is ${\theta }_{\mathrm{eq}}=\pi /2$, and coalescence is complete. For small angles ($\theta \ll 1$), Equation (8) reduces to the well-known expression $\dot{\theta }\simeq {\displaystyle \frac{2}{\tau \theta }}$ [69].

A major drawback in continuum approaches of tissues is that cells are not considered as individual entities in the model. This fact leads to important limitations when modelling cellular aggregates, for example, not accounting for the deformability of cells and their effect in the mechanics of the 3D aggregate. Most importantly, the phenomenological parameters used in the continuum description cannot be linked to cellular behaviour, which limits the predictive power of these models. Next, we discuss a different modelling approach, considering cells as agents that interact to form a cellular aggregate.

4.2. Agent Based Modelling of Tissue Fusion

With the arrival of microscopy techniques such as light-sheet [117] or two-photon [118] microscopy in the last decades, it is now possible to track the motion of individual cells with high spatial and temporal resolution within 3D tissues [119,120,121]. Notably, recent studies have been able to study single-cell movement and cell shape changes in 3D cell aggregates during tissue fusion [66,97]. The previous studies motivate the necessity of cell-based computational models of 3D tissues, where cell-cell interactions give rise to large-scale tissue movements. Such models consider cell-level processes such as cell adhesion, division, death, differentiation, and migration [122], allowing for the representation of heterogeneous cell populations and the modelling of complex tissue architectures. Furthermore, such models can integrate both local and long-range signalling, that, when coupled with mechanical forces, may lead to intricate emergent behaviours. Thus, one can develop computational models to engineer tissue shapes in silico, followed by experimental validation using in vitro stem cell systems.

Here, we focus on cell-based models for the study of 3D multicellular systems [69,97,99,122,123,124,125,126,127,128,129,130,131,132,133]. The simplest type of models are centre-based models (CBMs) [122,129], where each cell is effectively represented by its centre of mass. More realistic models that take into account cell shape are deformable cell models (DCMs) [99,122,123]. Finally, energy-based approaches also exist such as kinetic Monte Carlo models (KMC) [99,134]. Here we will restrict ourselves to force-based models, which provide a more intuitive way of understanding complex cellular interactions. Next we briefly present each of these models and their applications in the study of tissue fusion (see Figure 3) [69,79,97,99,102].

4.2.1. Centre Based Models—The Overlapping Sphere Model

Despite its simplicity, the overlapping sphere model (OS) [129] allows for a fairly detailed description of cell-cell interactions, making it a suitable method for simulating tissue fusion [79,102]. Each single cell in the aggregate is represented by its centre of mass; therefore, the collection of cells is modelled as a set of points $\{{\mathbf{r}}_1,{\mathbf{r}}_2,\dots ,{\mathbf{r}}_N\}$ in space, where N is the total number of cells in the system. Cells can then interact with each other via forces that act on the cell centres. Usually, inertial effects are neglected due to the low Reynolds number in the system [135], where viscous forces dominate over inertial forces. The equations of motion for each cell i read [122]:

$${\lambda }_s{\mathbf{v}}_i+\sum _{j\in {\mathsf{U}}_i}^{N_i}{\lambda }_{ij}({\mathbf{v}}_i-{\mathbf{v}}_j)=\sum _{j\in {\mathsf{U}}_i}^{N_i}({\mathbf{F}}_{ij}^s+{\mathbf{F}}_{ij}^a)$$

(10)

where ${\mathbf{v}}_i=d{\mathbf{r}}_i/dt$ is the velocity of the cell, ${\mathsf{U}}_i$ is the set of neighbours, $N_i\equiv \left|{\mathsf{U}}_i\right|$ is the total number of neighbours, ${\lambda }_s$ is an effective friction coefficient between cell i and the extracellular matrix, and ${\lambda }_{ij}$ is a cell-cell friction tensor [122]. Friction forces are represented on the left-hand side of the equality and must be balanced by the sum of passive forces (cell-cell adhesion/repulsion) ${\mathbf{F}}_{ij}^s\left(t\right)$ and active forces (e.g., cell-cell protrusions) ${\mathbf{F}}_{ij}^a\left(t\right)$ acting on cell i [79,102]. Equation (10) constitutes a linear problem described by a sparse symmetric matrix that can be efficiently solved by using a conjugate gradient method [122]. A radius $R_i$ is assigned to each cell i, and two cells are considered neighbours if a certain condition is met regarding the distance between their centres ${\mathbf{r}}_{ij}={\mathbf{r}}_i-{\mathbf{r}}_j$, such that $r_{ij}\equiv \left|\right|{\mathbf{r}}_{ij}\left|\right|<r_{max}$, with $r_{max}$ usually being a function of $R_i$. Below this threshold, cells will interact. Typically, ${\mathbf{F}}_{ij}^s=F\left(r_{ij}\right){\mathbf{n}}_{ij}$, where ${\mathbf{n}}_{ij}={\mathbf{r}}_{ij}/r_{ij}$ is the unit vector connecting cells i and j. The function $F\left(r_{ij}\right)$ meets the criteria of having a rest length or equilibrium distance $r_{eq}$ such that $F\left(r_{eq}\right)=0$. Furthermore, it is required that: $F>0$ for $r_{ij}<r_{eq}$ accounting for the limited compressibility of cells. If cell adhesion is desired to be included in the model, then it is additionally required that $F<0$ for $r_{ij}>r_{eq}$. Various force functions exist that meet these requirements, ranging from simple functions such as the linear spring [136,137,138,139] to more complex functions such as forces derived from Morse [140] or Lennard-Jones-like potentials [141,142], as well as Hertz forces [143,144]. On the other hand, the active forces ${\mathbf{F}}_{ij}^a$ are usually regarded as protrusive cellular interactions and are stochastic in nature [79,102,145]. Such forces introduce active cellular rearrangements in the tissue, which critically impact its mechanical state, enabling the system to transition from a solid to a fluid-like behaviour [79,102]. Furthermore, in combination with a continuum description of tissue fusion, simulations allow the study of how cellular behaviour impacts the mechanical parameters at the supracellular scale [79,102].

4.2.2. Deformable Cell Models

The main limitation of CBMs is that they do not account for cell shape. To address this problem, DCMs can be used whereby the cell is modelled as a deformable object. Two main strategies exist: (a) describing each cell as a collection of smaller sub-cellular elements (subcellular elements method (SEM) [108,123,125], or cellular particle dynamics (CPD) [69,70,99,130]), and (b) explicitly describing the cell membrane as a deformable surface (active cellular foams [97,131,133]).

Cellular Particle Dynamics and Subcellular Elements Method

SEM and/or CPD [69,70,99,108,123,125] follow a similar approach and can be thought of as a natural extension of CBMs. These methods take into account the volume and shape of the cells by dividing the cell into smaller elements (particles) that provide a coarse-grain description of the subcellular space. The key difference between the two methods lies in the choice of the force fields. In what follows we will only focus on the CPD method, which was used to simulate tissue fusion [69,99]. In CPD a cell is broken down into smaller sub-units called cellular particles (CPs), that are each represented by their centre of mass. The CPs interact with each other via short-range attractive and repulsive forces, and their dynamics account for changes in cell shape. In a typical fusion scenario where two aggregates come into contact, the system dynamics is determined by the forces acting on the CPs. Since these systems are found in the low Reynolds number regime, inertial terms can be neglected, and the CPs will undergo overdamped Brownian dynamics. The equation of motion of the ith CP found in the nth cell with position ${\mathbf{r}}_{i,n}$ follows the following Langevin dynamics:

$$\lambda {\displaystyle \frac{d{\mathbf{r}}_{i,n}}{dt}}={\mathbf{F}}_{1,i,n}+{\mathbf{F}}_{2,i,n}+{\mathit{\xi }}_{i,n}\left(t\right)$$

(11)

where $\lambda$ is a certain friction coefficient and ${\mathbf{F}}_{1,i,n}$ and ${\mathbf{F}}_{2,i,n}$ account for the intracellular and intercellular forces, respectively, that act on the cellular particle (CP) i of cell n. Intracellular random forces are modelled by using a stochastic force ${\mathit{\xi }}_{i,n}\left(t\right)$. For details on the specific form of the intracellular and intercellular forces; see Refs. [69,99,108]. Tissue fusion of two spherical multicellular aggregates was studied using CPD simulations [69,99,108]. By combining experiments and theory, the authors were able to relate cellular and tissue scale properties. Cell sorting was additionally studied in Ref. [99], showing that CPD can provide a realistic description of multicellular structure formation once correctly calibrated. Finally, the formation of toroidal rings and tubular structures through tissue fusion shown in Ref. [69] highlights the importance of the physical properties of 3D spheroids on bioprinting applications.

Active Cellular Foam Models

Active cellular foam models explicitly take into account cell shape via a discretisation of the cell surface into interconnected nodes forming a triangulated surface. This approach makes the cell surface deformable, allowing for changes in cell shape due to shear and compressive forces between cells [97,122,130,131,133]. While these models allow for a more accurate representation of cells where changes in cell shape can be explicitly accounted for, these benefits come at a significantly higher computational cost. Recently, a novel 3D active foam model was used to study tissue fusion [97]. In this work, cells were described as viscous shells that interacted with each other by means of adhesive tension and random protrusion forces. The authors found that complete fusion is linked to tissue fluidisation, which depends on several parameters such as cell-cell tension, persistence time of the protrusions, cortical relaxation time, and active motility [97]. Interestingly, cell-cell friction promoted jamming in the system, constituting a distinct physical mechanism leading to arrested coalescence.

5. Final Remarks and Future Directions

The phenomenon of arrested coalescence in tissue fusion stands as a simple valuable tool to characterise in a quantitative manner the rheological properties of 3D spheroids, distinguishing fluid-like from solid-like tissues. It is important to remark that the rheological behaviour of cellular systems is generally timescale-dependent [146], and care must be taken when interpreting rheological parameters. Indeed, different rheological methods might provide different values for the parameters if the timescale at which the system is probed differs. In this review, we covered the study of homotypic fusion, i.e., the fusion of identical 3D cellular aggregates. However, it is well known that tissues can have different mechanical affinities [76,105,147], leading to the formation of tissue boundaries, which are essential during tissue development [148]. We believe that theoretical models of heterotypic tissue fusion will be very valuable in the future to infer mechanical parameters describing the affinity betweeen different tissue types. Finally, the ultimate goal is to bridge tissue-scale behaviour to cellular behaviour. In this regard, some first steps have been done in the context of tissue fusion where agent-based computational modelling plays a major role [79,97,99]. We envision that only when tissue rheology can be predicted from cellular behaviour will the value of these technologies be important for clinical purposes [6].