TY - JOUR
AU - Qiu, Mingfeng
AB - 1 Introduction Tissues are basic forms of biological materials and constitute the organizational level between individual cells and entire organs. A typical tissue is the epithelium, a layer of cells lining the surfaces of animal organs. Epithelial cells are polarized along the apicobasal axis [1]. They bear locally specific proteins both on the surface and internally. At the basal surface, each cell is decorated with integrins, a family of specialized adhesive trans-membrane proteins. Integrins attach the cell to the basal lamina, a dense polymer substrate, which is a type of extracellular matrix (ECM) [2,3]. The apical surface, which can be associated with an apical ECM, is either in contact with the basal surface of another layer or forms an interface with a lumen or air external to the organ [4]. Confluence of the tissue is mostly ensured through adhesive trans-membrane protein complexes of the cadherin family in the apical region. They form cell-cell junctions such that cells are tightly packed together to form a monolayer. The shape and motion of biological tissues emerge from the underlying cellular structures and subcellular activities, resulting in more complicated dynamics than observed in passive materials [5–8]. Understanding the behavior of cells and tissues, such as how their size and shape are controlled by their subcellular activities and how large-scale events like convergent extension are coordinated, is central to the study of development, homeostasis and disease [9–16]. The vertex model is a successful theoretical framework to study epithelial tissue morphology at cellular resolution. It originates from a geometric description of cellular packings [17–19]. Inspired by a vertex-based study of dynamics of cellular patterns in soap foams [20], the epithelial model was later augmented with a mechanical description [21–23]. The vertex model abstracts an epithelial tissue as a tiling of polyhedrons [24–27]. A common dimensional reduction approach is a representation of the tissue in two dimensions (2D) where cells are reduced to polygons in an apical cross-section [23,27,28]. In the resulting network, edges represent junctions between neighboring cells and vertices correspond to sites where three or more cells meet. The polygonal configuration is commonly governed by a phenomenological Hamiltonian, a function describing the energy of the system dependent on the geometry of individual cells. The energy function in the seminal model by Farhadifar et al. [23] represents three major force contributions in the cell: the contractile forces generated by the cytoskeleton when myosin motor proteins walk on actin filaments, the resistance from the incompressible cytosolic fluid and the adhesion between neighboring cells. In addition, this framework also includes a set of update rules that account for topological changes in the network that are necessary to perform cell division, cell extrusion, and cell neighbor exchanges. Vertex models like the one of Farhadifar et al. [23] offer a theoretical foundation to bridge subcellular activities, in [23] actomyosin contractility and a cell cycle, to tissue-scale behavior. The vertex-based framework provides an intuitive discretization of a confluent tissue on a subcellular level. Vertex models have been successful in exploring the properties of epithelial tissues. They are able to relate cell packing geometry to the balance of physical forces and to predict a rigidity transition of the tissue in response to changes in single-cell properties [23,29–31]. Vertex models have also been applied to model a variety of specific morphogenetic processes, e.g., Drosophila germband extension [32,33], ventral furrow formation [34], dorsal closure [35,36], and wound healing [37]. Two seminal reviews have been dedicated to the vertex model [27,28]: The first gives a detailed overview of the formulation first introduced in [23] and its applications by a wide community of modelers [28]. A second review presents a general form of the model in the language of virtual work and classifies different vertex models based on the spatial dimensions and geometry that they cover [27]. Both reviews have also highlighted extensions to the vertex model, such as incorporating planar cell polarity and the generalization to a three-dimensional (3D) geometry. A recent book has covered the vertex model, focusing on a didactic presentation of the construction of the model and earlier developments in the field [38]. In recent years, extensions and modifications of the vertex model, specifically aiming to capture subcellular structures, dynamics, and feedback mechanisms, are being explored, highlighting the need for a coherent picture of model variants and their respective properties. By now, the vertex model has been established as one of the most popular models to study epithelial tissues from a physical perspective, with an increasingly large and interdisciplinary community. The vertex model introduced by Farhadifar et al. [23], popularized by its capability to capture a solid-to-fluid transition in tissues [29,39], is by far the most cited vertex model study and recently being referred to as the “conventional model” or “standard model” in the realm of vertex models [40–42]; it serves as a reference model in this review. Meanwhile, alternative and extended versions of the model have been developed, partly driven by advances in experimental tools like genetic engineering and imaging techniques. Novel details on the motor proteins and related molecular processes at the nano- and micro-scale are being discovered, which the reference model does not capture. For example, while the reference model only implements uniform contractility along the cell perimeter, during development, actomyosin is first localized in the medial region before later on localizing to the cortex along cell-cell junctions [43]. To account for these experimental observations, a large number of variants have been proposed to extend the vertex model, ranging from stochastic parameters in the energy function to coupling with biochemical networks, e.g., [34,36,44,45]. However, recent studies do not always explore physical properties and compare new model variants to properties of the reference model, and it remains unclear what impact extensions have on the theory of vertex models. Here, we review novel extensions of the vertex model in the context of their motivations and discuss insights on the fundamental properties of vertex models obtained by these extensions. We build on frameworks set by previous works and reviews and categorize model extensions that address different complexity limits of the vertex model. Our objective consists of three specific goals: First, we aim to survey the latest model extensions and to provide a topical reference for vertex model practitioners to make informed modeling choices, when studying specific biological and physical questions, or for developing further extensions. Second, we strive to make this review an accessible and up-to-date introduction for researchers new to this field, including both theoreticians and experimentalists. Finally, we summarize key questions to be studied and debated. We find that a perspective attracting increasing attention is the notion of bridging spatial scales: coupling the vertex model to subcellular activities and studying emergent behavior on the tissue scale. Incorporating active regulation at smaller scales into a scalable model framework provides many opportunities and challenges for future studies. In structuring our review, we start by introducing what we consider our reference model followed by extensions to the framework. In each section, we highlight experimental observations and theoretical considerations supporting a particular extension, describe the mathematical formulation, and give examples of applications. Specifically, we cover non-uniform and adaptive contractility, coupling with cell-cell junctional remodeling mechanisms related to mechanosensitivity, active topological transitions, and tissue rheology. Then we examine how various model extensions preserve or modify tissue properties, the most well-known of which are tissue rigidity and cell packing geometry. We highlight the lack of comparative studies and unifying frameworks in this aspect. Finally, we put the model in the broader context of complementary models and methods to study tissue behavior. As extensions with a subcellular focus have been mostly implemented in 2D apical vertex models so far, we will restrict our introduction of the model formalism to 2D and only explicitly highlight three-dimensional aspects if they are relevant to the implementation or conclusions regarding subcellular aspects. 2 A vertex model for epithelial tissues In preparation for our review, we briefly recall the formulation and properties of the vertex model that the reviewed extensions are based on [23,29,39]. In a confluent cellular tissue such as an epithelial monolayer, cells are bound together through adherens junctions – adhesive molecular complexes at cell-cell boundaries. Tension can be generated by the contractile actomyosin cytoskeleton along the perimeters of the cells and propagated through the adherens junctions. Both the actomyosin ring and the adherens junctions are located close to the apical cortices of the cells. This allows, in many cases, a description of a confluent tissue as a 2D polygonal network, with each edge representing a cell-cell junction consisting of the adherens junctions, the cell-cell boundary, and the actomyosin cytoskeleton on both sides of the boundary (Fig 1A and 1B). Vertices represent junctions where three or more cells meet, and their coordinates, the degrees of freedom of the model, are subject to the physical law of force balance (Fig 1B and 1C). In a vertex model simulation, one usually proceeds by iterating between two steps: (1) update vertex positions based on the forces exerted by the neighboring cells (1C); (2) perform topological transitions in the vertex network that allow for cell neighbor exchanges, cell extrusions, and cell divisions based on a set of rules (Fig 1D). In the following, we will introduce the two steps specifically for the reference model [23,39] and we will briefly discuss the resulting properties of the model [29]. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. Forces, topological transitions, and phase transitions in a vertex model for confluent epithelial tissues. (A) The vertex model incorporates what are considered the main mechanical components in polarized epithelial cells: cytosol filling the volume of a cell; in the apical region: an actomyosin ring along the cell perimeter and adherens junctions connecting two cells, and a cell membrane enclosing the volume. The two-dimensional model is a cross-section close to the apical surface of the tissue. (B) Cells in the tissue are approximated by polygons that form a confluent tiling. In the vertex model, we attribute specific forces to the mechanical components: The cytosol resists compression and expansion forces; the actomyosin ring exerts contractile forces on a cell and the adherens junctions result in attractive forces along the cell-cell boundary. The adhesion forces compete with the cortical tension from the cell membrane and actomyosin ring. (C) Vertices in the network satisfy force balance based on the forces from the mechanical components. Effective forces on a vertex i are along the edge length vector pointing away from reference vertex i (tension) and the normal vector defined as perpendicular to the line connecting the two neighbors of vertex i and pointing towards the outside of the cell (pressure). (D) A common set of update rules accounts for topological transitions in the temporal evolution of a vertex network. T1 transitions are cell neighbor exchanges where an edge between cells a and b falls below a certain threshold and continues to shrink until it vanishes; then, a new edge between cells c and d is inserted perpendicular to the old edge. T2 transitions are cell extrusions where a cell with an area below a threshold value vanishes and leaves a single vertex. Cell division, as part of proliferation, occurs when a cell splits into two. (E) Phase diagram of tissue mechanical states as a function of parameters kP (perimeter contractility) and kl (line tension). When the control parameter is below (white line), the ground state is a regular hexagonal lattice, which becomes unstable once p0 increases past [23,39]. For a disordered, metastable tissue, a solid-to-fluid transition occurs as p0 increases past , characterized by vanishing energy barriers for T1 transitions [29]. The exact threshold of the solid-to-fluid transition can vary depending on cell packing disorder [10] and implementation details of the dynamics [46,47]. The gray area marks where the model breaks down due to vanishing areas. https://doi.org/10.1371/journal.pcbi.1012993.g001 First, we introduce the equations of motion that describe the dynamics of the vertices of a network of epithelial cells. As the considered tissue motion is typically in a regime where viscous forces dominate inertial forces, it is assumed that the system is governed by over-damped dynamics: (1) Here, the left-hand side is a friction force with vertex velocity derived from vertex position . is an effective drag coefficient that accounts for dissipation from external sources. External dissipation arises from interactions with the environment, e.g., the ECM polymer substrate. While external friction is commonly used to model system friction and internal dissipation is neglected, more general formulations are possible. We discuss dissipation and its effect on time scales and how recent models implement this in Sect 4. The friction force balances , the sum of all other forces exerted on vertex i. Two approaches have been used to formulate the force equations driving vertex motion. One proposes explicit functional forms of the forces by considering the membrane tension and the cortical pressure directly [21]. A more popular approach involves specifying a phenomenological Hamiltonian to derive [22,23]. While this phenomenological Hamiltonian is generally not a conserved total energy, it has proven useful in describing tissue mechanics in two ways. First, it is well documented that the ground states of the phenomenological Hamiltonian, i.e., the most relaxed configurations of the network, display fundamental properties also observed in tissues (Fig 1E) [23,29,39], from which the vertex model’s popularity today likely arises. Second, the phenomenological Hamiltonian serves as a mathematical tool to derive mechanical forces, providing an extendable framework for modeling realistic tissue morphologies [22,27,28]. Compared to the force-based formulations of the vertex model, the derivative of the phenomenological Hamiltonian results in qualitatively similar forces but is not strictly equivalent [28]. For our review, we will follow the approach based on the phenomenological Hamiltonian due to its more widespread use. We will refer to the phenomenological Hamiltonian as an energy function, a terminology widely adopted in the field to emphasize the analogy to energy in passive systems [27,28]. Let be the energy of the system dependent on all vertex positions x. Minimizing E with respect to the coordinates of vertex i, we can rewrite the equations of motion as (2) where denotes the gradient with respect to the coordinates of vertex i. It is often assumed that mechanical equilibrium is reached on a time scale much shorter than any driving process of tissue deformation [23]. Under this quasi-static assumption, the equation of motion simplifies to (3) The quasi-static Eq 3 can be solved by directly minimizing the energy function [23,39]. An alternative method, that has been used already in earlier works, is to evolve Eq 2 to steady state [22,48], in which case the friction term can be interpreted as a numerical relaxation method. In the following, we describe the energy function that is commonly used in 2D apical vertex models to phenomenologically describe the intra- and inter-cellular mechanics of epithelial tissues [23]: (4) The three terms correspond to three types of forces typically considered to contribute to the sum of forces acting on vertex i (Eq 1). The first term, summing over all cells a, penalizes the change of the current cell area Aa from a predefined target cell area A0, prefactored by an elasticity coefficient kA. This term, referred to as area elasticity, represents bulk forces resisting area or volume changes of the cell (blue arrows in Fig 1B). Examples include the hydrostatic pressure from the cytosol or osmotic pressure across the cell membrane. The second term is linear with respect to the cell boundary length and represents an effective surface energy between adjacent cells. This includes both cortical tension (linear part of the perimeter elasticity, red arrows in Fig 1B) and adhesion forces by the adherens junctions (yellow arrows in Fig 1B). It sums once over all cell-cell junctions, denoted by , an unordered pair of neighboring cells a and b. kl is the tension coefficient, is the length of the contact surface between cells , and vertices i and j are connected by edge . If cortical tension dominates adhesion, then kl > 0, while if adhesion dominates, kl < 0. The area elasticity and line tension terms were already proposed in an earlier vertex model by Nagai and Honda [22]. The extension of this energy function with a perimeter contractility term was first proposed in the seminal study by Farhadifar et al. [23] to better capture a characteristic feature of epithelial cells: the contractile actomyosin cytoskeleton that forms during embryonic development at the apical cortex of a cell and can be found as a belt-like structure in later stages of development and in mature tissues [23,43]. In the energy function, the perimeter contractility term represents a second-order effect of the cytoskeleton where kP is the contractility coefficient and is the perimeter of cell a. Taking the partial derivative of the energy function with respect to a vertex position , we can write the total force on vertex i (Fig 1C) as: (5) where is the set of cells containing vertex i and is again an unordered pair of neighboring cells. The expression is a decomposition of the force into stresses , and strains , . To simplify the expression, we define and . is the unit vector along the edge between cells a and b, pointing away from vertex i, la,i is the length of the line segment connecting the two adjacent shoulder vertices to vertex i in cell a (dashed line in Fig 1C), and is a unit vector perpendicular to the line segment between the shoulder cells, pointing toward the outside of cell a (Fig 1C). This decomposition can be generalized and is instrumental in the statistical inference of cell and tissue stresses from experimental data [49–51]. Using Eq 4, the mechanical stresses can be interpreted as a cell pressure from area conservation and an edge tension from cell-cell adhesion and cortical tension: (6)(7) In the following, we introduce the second simulation step in the vertex model framework which aims to capture cell rearrangements in a confluent tissue. In the reference model, and many other models, these are implemented as update rules acting on the network’s topology [23,28,39]. Adopting terminology from foam physics, where similar phenomena exist [5], the first two transitions are called T1 and T2 which describe cell neighbor exchanges and cell extrusions, respectively. The third transition is specific to biology: a cell divides into two daughter cells. Here, we outline a set of specific procedures for performing the three common topological transitions in the vertex model (Fig 1D) [23,28,39]: Cell neighbor exchange, also referred to as a T1 transition: If , where is a length threshold below which the bicellular junction is considered disintegrated, delete edge ab and merge attached vertices. Then, split the vertex again and create a new edge with a length perpendicular to the old edge. Cell extrusion, also referred to as a T2 transition: If , where is a minimum cell area, merge all vertices that define cell a and eliminate cell a. Cell division: If , where is a critical cell area beyond which a cell divides, split the cell a into two new cells with a new edge, which intersects the original cell center at random orientation and creates two new vertices. The cell division rule is usually accompanied by a more or less simplistic cell cycle model, e.g., a time-dependent deterministic target cell area A0(t) that is increasing monotonically in time. In such a case, the decision threshold is set, e.g., to double the initial target area . It is evident that the transition rules involve sharp and artificial threshold values such as the minimum edge length for triggering a transition, and the length of the new edge after a T1 rearrangement. For now, we can interpret the thresholds as bound to precision limits in experimental observations. We note that the algorithms mentioned above are only exemplary for the three typical topological transition rules in vertex models and a variety of implementations exists [39,52]. We will discuss nuances of implementation choices in Sect 4. Last, we highlight a physical perspective of the vertex model that has illuminated avenues to characterize and change material properties of a tissue [23,29,39,53,54]. One can rewrite the energy function in Eq 4 as (8) by omitting a constant and defining a target cell perimeter where the factor 2 accounts for the convention to only sum over each edge once [29,39,55]. In this alternative formulation, two important phase transitions are controlled by the target cell shape (Fig 1E): (1) Below , the ground state of the tissue, i.e., the global minimum of the energy function, is a lattice of regular hexagons. When p0 increases past , the regular hexagonal ground state becomes energetically unstable to deformation [23,39,56]. (2) In a disordered multicellular tissue, which is a metastable network in the vertex model, a solid-to-fluid transition occurs as p0 increases past a critical threshold [29]. This transition is characterized by the vanishing energy barriers for cell rearrangements, particularly T1 transitions, allowing the tissue to flow. As the two adjacent hexagons in a four-cell unit become pentagons while undergoing a T1 transition (Fig 1D), the shape index of the solid-to-fluid transition approximately corresponds to that of a regular pentagon [29]. The exact value of can vary depending on the extent of cell packing disorder [10] and implementation details such as topological transition rules and energy minimization schemes [46,47]. Based on the reference model introduced here, we discuss extensions to this reference model in the next section. 3 Subcellular extensions of the vertex model 3.1 Spatiotemporal dependence of model parameters Protein activities in tissues are subject to thermal fluctuations and active regulations. In particular, the actomyosin cortex is subject to dynamic processes giving rise to contractile forces dependent on molecular turnover (Fig 2B) [57,58]. To incorporate temporal and spatial inhomogeneities in the mechanical properties of bicellular junctions, a straightforward generalization of the vertex model is to allow for a time-dependent, edge-specific line tension coefficient in the energy function with the perimeter contractility term (Eq 4) or the target perimeter in the perimeter elasticity formulation (Eq 8) [44]. In the following, we categorize extensions that introduce parameter heterogeneity in the reference model by their spatial or spatiotemporal variation and group them according to their governing method: a stochastic model, a deterministic rule-based model, or a deterministic model-free observation-based pattern (Fig 2A). We summarize these extensions in Table 1, highlighting different model variants, their biological motivations, specifications of relevant parameters and conclusions from specific studies. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. Spatiotemporal dependence of model parameters and alternative cell contractility models. (A) Different models with spatiotemporally dependent coefficients aim to represent observed effects of (B) the complex recruitment dynamics of myosin to actin fibers in the cytoskeleton. (A) Stochastic models can describe spatial static inhomogeneity of myosin distribution. This can be expanded towards temporal stochastic processes (e.g., Ornstein-Uhlenbeck) where myosin abundance on an actin fiber varies over time. Deterministic models impose a specific condition for myosin distribution, e.g., with respect to an angle or a clock. These models are directly motivated to capture observations like convergent extension or active contractility. Some observed patterns are implemented model-free into vertex model parameters. A representative example is the contractile cable along the perimeter of the amnioserosa during dorsal closure of Drosophila embryo development. (C) The actomyosin contractility mode expressed in a vertex model with a perimeter contractility term [23] is sometimes called a purse-string mechanism; it represents a continuous actomyosin belt along the cell perimeter and generates uniform contractility. Apart from that, a series of non-uniform contractility models exist. The edge model assumes a discontinuous actomyosin belt and thus reduces coupling of the contractility structures [80]. The medial cytoskeleton observed in premature tissues has been modeled as isotropic radial fibers and anisotropic parallel fibers [35,81,82]. Parallel fibers carry a load and resist isotropic extension without introducing additional degrees of freedom [82]. https://doi.org/10.1371/journal.pcbi.1012993.g002 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. Overview of vertex model extensions (Part 1). https://doi.org/10.1371/journal.pcbi.1012993.t001 To capture spatial variation of myosin activity at different cell-cell junctions, kl,ab can be drawn from a Gaussian distribution with mean value kl,0 and variance [49,50] or from a uniform distribution [59]. Independent of the specific choice of a distribution, stochastic target perimeters result in increased shear modulus and rigidity of tissues and lead to a heterogeneous solid state with distinct percolation properties which could be relevant for tumor invasion dynamics [44]. To capture temporal fluctuations in tension, for example due to protein binding and unbinding, kl,ab(t) can be described by an Ornstein-Uhlenbeck process (9) where is a Gaussian white noise with zero mean and non-correlation between two edges ab and cd, − t2) where is the relaxation time of the process [7,60–63]. Temporal fluctuations in tensile forces drive fluidization in tissues characterized by loss of correlation between tension and strain [37,64,65]. In other cases, spatial and temporal inhomogeneities are required to be deterministic. A classical example of deterministic spatial inhomogeneities is angle-dependent tension, which can be used to model planar cell polarity – the polarized distribution of protein activities within the plane of the tissue – as well as planar cell chirality, which involves the intrinsic asymmetry and directional movement of cells within the tissue. Anisotropic tensions can explain tissue convergent extension [10,33], the increase in the stability of four- and higher-order rosette vertices [6,66] and asymmetrical cell intercalation and organ rotation [67–69]. Spatial inhomogeneities are also often specified to a desired pattern. For example, linear bundles of actomyosin cytoskeleton, termed cables, span multiple cells. They play critical roles in morphogenesis, often contributing to closure and invagination events [54,70–75]. These supracellular cables can be modeled by imposing a higher tension coefficient along a series of consecutive edges [76–78]. Deterministic temporal programs have already been used in the model by Farhadifar et al. to include a cell cycle via increasing target cell areas before a cell division occurs [23], although the activation of a division cycle is initiated stochastically. Similarly, the decrease of target cell perimeters has been a successful model to capture the material properties of the amnioserosa of Drosophila during dorsal closure. Here, a combination of heterogeneous and decreasing target perimeters arrests fluidization of the tissue [36]. This suggests that stochasticity combined with deterministic growth or shrinkage is a common motif to explain observed tissue properties. A comparison between driven oscillations of target areas A0(t) and target perimeters P0(t) finds that area oscillations that represent the medial activity of actomyosin tend to produce more anisotropic cells, suggesting that different pools of contractility may produce distinct cell morphologies [79]. 3.2 Alternative cell contractility models The continuous actomyosin belt represented in the perimeter contractility model (Eq 4) produces uniform contractility along the cell boundary. In epithelial cells, however, the actomyosin displays distinct subcellular localizations that vary with confluence, functions, and developmental stages. Mature, confluent epithelial tissues predominantly exhibit a continuous, circumferential actomyosin structure along the cell perimeter, associated with adherens junctions [43]. In contrast, during early developmental stages, when confluency may be temporarily disrupted for large-scale tissue reorganization [87], actomyosin structures appear more medial, discontinuous, and heterogeneous—as has been observed in gastrulation of chick embryos [87] and Drosophila embryos [43,85]. Similar actomyosin structures have been observed in cell cultures of human intestinal [88] and Madin-Darby canine kidney (MDCK) cells [89]. These formations, often referred to as the medial pool of actomyosin, are essential for pulsed contractions and apical constriction [43,85]. To capture the aforementioned observations, some models incorporate non-uniform contractility (Fig 2C). Here, we consider two main types of these extensions, namely edge-based contractility and various medial contractile structures (summarized in Table 1). Firstly, tricellular junctions are thought to pose a physical obstacle to the transport of proteins along the cell perimeter, and thus the contractile force may not be uniform along the cell perimeter [90–93]. Experimental studies indicate that each bicellular junction functions as an independent contractile unit, where actin cables are anchored end-on to cadherin complexes binding to tricellulin of the tricellular tight junctions ensuring tissue confluency [94,95]. Vertex models with edge-based contractility terms aim to individually describe bicellular junctions of a cell [32,34,35,55]. In an edge-based contractility model, the last term of the vertex model energy function is modified (10) Here, in the third term, contractility is based on individual edges, with an elastic coefficient . Due to the quadratic expression, this energy function is different from the perimeter contractility model (Eq 4), even when the line tension coefficient kl and the contractility coefficient are constant for all edges. Analogous to the reformulation of the perimeter-based model (Eqs 4 and 8), one can absorb the line tension into the quadratic term in Eq 10 and introduce a target edge length l0. In doing so, the energy function takes a form, equivalent up to a constant: (11) In 2D apical vertex models embedded in 3D space, an anisotropic edge-based contractility induces chiral cell sliding driving asymmetric twisting similar to what is observed in organs like the heart [83,84]. A vertex model with both perimeter contractility and edge contractility can capture mixed effects of contractile actomyosin at the cell border [80]. Further, some researchers incorporate additional contractile structures into the model to represent force-generating mechanisms that are not confined to cell perimeters, such as the star-like formations recently observed on the basal side of epithelial cells [96]. When studying the oscillations of cell areas in the process of dorsal closure during Drosophila development, a series of studies add radial stress fibers, termed spokes, connecting the vertices with the center of each cell [32,35,45,81]. The spokes are implemented as contractile, elastic intracellular edges. In addition to representing the medial pool of actomyosin, they also play a similar role to pressure forces due to area elasticity. These stress fibers are still isotropic in the sense that they are oriented along different directions well distributed in space provided that cell shapes are not strongly elongated. Some other studies explore the effects of anisotropic contractile structures. Parallel stress fibers aligned with imposed unidirectional stretching are employed in an effective force-based vertex model to explore the scaling of the number of fibers with cell apical areas [82]. Similarly, the dynamics of various subcellular stress fiber geometries in a stress-fiber reinforced vertex model reveal a size-dependent response of cells [97]. One comparative study confirms that edge-based contractility models belong to the same sub-class of vertex models as perimeter-based contractility models, exhibiting a crossover from soft to stiff network behaviors [42]. Other studies suggest that different contractility modes and structures might influence cell packing geometry and the threshold of a rigidity transition [79,82]. Another way to incorporate anisotropic contractility from the medial myosin pool is to write an anisotropic continuous stress in the bulk of the cells and derive their effective forces on the cell-cell junctions and vertices [99]. This type of model has been used to investigate how epithelial colonies elongate through collective effects, with tissue elongation resulting from anisotropy in the average cell elongation [99], and how the V-shaped muscle compartments emerge in zebrafish, revealing that spatially modulating the mechanical environment around and within tissues can lead to complex organ shape formation []. In these two examples, the anisotropic bulk stress is crucial for symmetry breaking during tissue growth and sculpting tissue shapes, respectively. However, recent comparative work has shown that the two different formulations that have been employed to derive the forces on vertices [99], lead to distinct tissue dynamics due to differences in their treatment of non-affine deformations [100]. Moreover, neither description of anisotropic bulk stresses could fully represent equivalent forces arising from cell-cell junctions in the reference vertex model [100]. Overall, these results suggest that different modeling choices could lead to distinct predictions, and it is important to align model design with experimental observations of contractile structures, which vary across biological contexts. 3.3 Subcellular regulation of cell contractility and mechanosensitive response The reference model [23] describes cell-cell junctions as elastic springs driven to equilibrium by forces derived from an energy function Eq 4. In other words, the mechanical forces along the edges in the model are determined by the deformation of the polygonal tiling through a positive linear relationship. This simplified representation of the mechanical behavior of cell-cell junctions is hardly sufficient considering that the junctions are tightly regulated by complex, interacting subcellular elements such as cytoskeleton dynamics, mechanotransductive and other signaling pathways, and gene expressions. Two early seminal studies, using a lateral 2D vertex model, explore how deformation-dependent edge target length could lead to apical constriction and epithelial invagination [48,101]. The showcased mechanical self-organization suggests that tissues might leverage mechanical signals such as force or deformation to actively regulate their contractile forces. This type of feedback loop can result in spontaneous tissue configuration changes, a behavior loosely referred to as mechanosensitivity. With the advance in experimental and imaging techniques, knowledge of how cells actively control their mechanics on the subcellular scale provides increasingly solid foundations to the mechanosensitive hypothesis (Fig 3B) [11,82,91,102–110]. Together with the accumulation in experimental evidence, theoretical interest in mechanical self-organization has also grown rapidly [111]. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. The effective relation between forces and cell shapes in the vertex model can be changed by introducing submodules featuring active regulation and mechanosensitive feedback describing activities at the subcellular scale. (A) The mechanics in the reference vertex model [23] are governed by an energy function (a phenomenological Hamiltonian) that dictates a linear relation between forces and cell shape (observed e.g., in the ground states or a proliferation-free simulation). To be able to capture realistic tissue morphologies, vertex models are commonly extended to a full simulation framework with components refining the subcellular scales, resulting in effectively nonlinear relations between forces and cell geometry [65]. The subcellular activities include the activation and turnover of cytoskeleton machinery and adhesion proteins, the relaying of mechanical signals through mechanotransduction molecules, the modulation of mechanical forces by the spatiotemporal dynamics of signaling pathways, and the activation and deactivation of genetic programs instructing protein synthesis. (B) Examples of submodules that actively regulate quantities are strain-dependent remodeling of the target edge length, edge length-dependent active tension, or strain-dependent remodeling of the density of a biomolecule, e.g., myosin. In systems with mechanosensitive feedback, these remodeled quantities are usually interdependent and formulated as coupled differential equations. https://doi.org/10.1371/journal.pcbi.1012993.g003 The vertex model is a convenient tool to incorporate junction regulatory effects (Fig 3A). Cell division, as implemented in the reference model [23], is already an active event that reflects the cell cycle, an effect beyond the passive linear junctional tension dictated by the energy function, although it does not form an explicit feedback loop. One of the earlier studies combining vertex modeling and experiments on how mechanosensitivity contributes to morphogenesis encodes a feedback loop by letting the cell target parameters depend on the strain and suggests that strain-enhanced apical constriction is critical in the invagination of an optic-cup organoid [112]. Notably, the authors already deploy a 3D vertex model. This approach of implementing feedback loops in mechanics-dependent model parameters is a generalization of our first category of extensions with spatiotemporal dependence of model parameters. Meanwhile, it can be seen as a long-time limit of other implementations using differential equations to describe the evolution of model parameters in time, representing subcellular scale activities to study their effects at the tissue level. In more recent work, one approach, using the edge-based elasticity model (Eq 11), is to introduce remodeling of the target edge length l0,ab, described by the following generic formulation (12) where f is a general edge regulation function. We use the subscript ab for all variables to emphasize that they can be edge-specific. On the left-hand side of the equation is the absolute rate of change of the rest length [113], or alternatively, can also be the relative rate of change [34]. On the right-hand side, a variety of regulation mechanisms can be incorporated by adopting different forms of f. For example, the actomyosin cytoskeleton is well-known to be viscoelastic [114]. Thus, describes linear Maxwell viscoelasticity with viscosity . More importantly, the target edge length can be remodeled by the strain or strain rate [65,113], which is an active regulation mechanism. This type of feedback is proposed to explain the mechanical ratcheting of cell-cell junctions [113]. A recent variant of this model shows that strain-regulated tension, together with mechanical memory dissipation implemented through the fluctuation of tension coefficients and the relaxation of tension and strain, leads to transiently stable higher-fold vertices (rosettes) [65]. Notably, in contrast to other mechanosensitive models, this model incorporates negative feedback, where tension decreases under edge extension and increases under edge contraction. The negative feedback is found to stabilize rosettes [65]. Biochemistry can be more explicitly described by considering the dependence of f on the quantity or activity of a certain molecule, generically written as mab. In this case, an additional evolution equation is required to govern the dependent variable mab, which can be written in the form (13) In one example, mab denotes the average edge density of myosin motors, with Eq 12 describing the active walking of myosin motors contracting actin filaments under low tension and the slipping of motors under high tension, while the myosin density mab, in turn, depends on the strain rate of the edge through Eq 13 [34]. In this case, tissue viscoelasticity and cell size variability emerge from the feedback between the strain rate of cell-cell junctions and myosin recruitment, coupled with junctional viscoelastic relaxation [34]. In a 3D vertex model, the effect of remodeling several mechanical parameters – the surface area elasticity, the target surface area, and the surface tension – on apicobasal polarity is investigated [115]. By combining a vertex model with the biochemistry of key proteins, the study shows that the mechano-polarity feedback loop can drive spontaneous tissue folding. Besides remodeling the target edge length or the target cell perimeter, an additional or alternative scheme is to augment the edge tension with an active force. With the decomposition into cell pressure and edge tension (Eq 5), one can extend the edge tension as (14) where the additional active force is regulated by phenomenological laws describing subcellular processes. Similar to f in the edge target length remodeling formulation, can depend on an array of factors. For example, an additional contractile force can be added, which becomes active only when the cell area exceeds a threshold. This feedback enhances tissue strength and resistance to rupture [116]. Others assume that the active force is dependent on additional quantities, such as a myosin activity level (e.g., [35]). Persistent cell neighbor exchanges and convergent extension could arise from this type of mechanochemical feedback [32,45,117]. The active force term can also depend on a cell polarity vector. The coupling results in a chiral patterning when the tissue is confined on a spherical surface, demonstrating how mechanosensitivity interacts with topology [118]. In these cases, again new rate equations need to be introduced to govern the additional variables, either itself or the biochemical species mab. In general, both positive and negative feedback loops can be easily incorporated into the outlined framework. Some researchers consider further complexities in the regulation of cell mechanics. For example, the coupling of a lateral inhibition model with a vertex model reveals how mechanics and biochemistry combine to instruct the alternating cellular pattern in the inner ear [119]. With further progress in imaging technology and discovering molecular processes underlying the transmission and interpretation of mechanical signals by biological systems, we believe that with these extensions (summarized in Table 2) we are only seeing the beginning of what might become a major part of vertex model studies. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 2. Overview of vertex model extensions (Part 2). https://doi.org/10.1371/journal.pcbi.1012993.t002 3.1 Spatiotemporal dependence of model parameters Protein activities in tissues are subject to thermal fluctuations and active regulations. In particular, the actomyosin cortex is subject to dynamic processes giving rise to contractile forces dependent on molecular turnover (Fig 2B) [57,58]. To incorporate temporal and spatial inhomogeneities in the mechanical properties of bicellular junctions, a straightforward generalization of the vertex model is to allow for a time-dependent, edge-specific line tension coefficient in the energy function with the perimeter contractility term (Eq 4) or the target perimeter in the perimeter elasticity formulation (Eq 8) [44]. In the following, we categorize extensions that introduce parameter heterogeneity in the reference model by their spatial or spatiotemporal variation and group them according to their governing method: a stochastic model, a deterministic rule-based model, or a deterministic model-free observation-based pattern (Fig 2A). We summarize these extensions in Table 1, highlighting different model variants, their biological motivations, specifications of relevant parameters and conclusions from specific studies. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. Spatiotemporal dependence of model parameters and alternative cell contractility models. (A) Different models with spatiotemporally dependent coefficients aim to represent observed effects of (B) the complex recruitment dynamics of myosin to actin fibers in the cytoskeleton. (A) Stochastic models can describe spatial static inhomogeneity of myosin distribution. This can be expanded towards temporal stochastic processes (e.g., Ornstein-Uhlenbeck) where myosin abundance on an actin fiber varies over time. Deterministic models impose a specific condition for myosin distribution, e.g., with respect to an angle or a clock. These models are directly motivated to capture observations like convergent extension or active contractility. Some observed patterns are implemented model-free into vertex model parameters. A representative example is the contractile cable along the perimeter of the amnioserosa during dorsal closure of Drosophila embryo development. (C) The actomyosin contractility mode expressed in a vertex model with a perimeter contractility term [23] is sometimes called a purse-string mechanism; it represents a continuous actomyosin belt along the cell perimeter and generates uniform contractility. Apart from that, a series of non-uniform contractility models exist. The edge model assumes a discontinuous actomyosin belt and thus reduces coupling of the contractility structures [80]. The medial cytoskeleton observed in premature tissues has been modeled as isotropic radial fibers and anisotropic parallel fibers [35,81,82]. Parallel fibers carry a load and resist isotropic extension without introducing additional degrees of freedom [82]. https://doi.org/10.1371/journal.pcbi.1012993.g002 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. Overview of vertex model extensions (Part 1). https://doi.org/10.1371/journal.pcbi.1012993.t001 To capture spatial variation of myosin activity at different cell-cell junctions, kl,ab can be drawn from a Gaussian distribution with mean value kl,0 and variance [49,50] or from a uniform distribution [59]. Independent of the specific choice of a distribution, stochastic target perimeters result in increased shear modulus and rigidity of tissues and lead to a heterogeneous solid state with distinct percolation properties which could be relevant for tumor invasion dynamics [44]. To capture temporal fluctuations in tension, for example due to protein binding and unbinding, kl,ab(t) can be described by an Ornstein-Uhlenbeck process (9) where is a Gaussian white noise with zero mean and non-correlation between two edges ab and cd, − t2) where is the relaxation time of the process [7,60–63]. Temporal fluctuations in tensile forces drive fluidization in tissues characterized by loss of correlation between tension and strain [37,64,65]. In other cases, spatial and temporal inhomogeneities are required to be deterministic. A classical example of deterministic spatial inhomogeneities is angle-dependent tension, which can be used to model planar cell polarity – the polarized distribution of protein activities within the plane of the tissue – as well as planar cell chirality, which involves the intrinsic asymmetry and directional movement of cells within the tissue. Anisotropic tensions can explain tissue convergent extension [10,33], the increase in the stability of four- and higher-order rosette vertices [6,66] and asymmetrical cell intercalation and organ rotation [67–69]. Spatial inhomogeneities are also often specified to a desired pattern. For example, linear bundles of actomyosin cytoskeleton, termed cables, span multiple cells. They play critical roles in morphogenesis, often contributing to closure and invagination events [54,70–75]. These supracellular cables can be modeled by imposing a higher tension coefficient along a series of consecutive edges [76–78]. Deterministic temporal programs have already been used in the model by Farhadifar et al. to include a cell cycle via increasing target cell areas before a cell division occurs [23], although the activation of a division cycle is initiated stochastically. Similarly, the decrease of target cell perimeters has been a successful model to capture the material properties of the amnioserosa of Drosophila during dorsal closure. Here, a combination of heterogeneous and decreasing target perimeters arrests fluidization of the tissue [36]. This suggests that stochasticity combined with deterministic growth or shrinkage is a common motif to explain observed tissue properties. A comparison between driven oscillations of target areas A0(t) and target perimeters P0(t) finds that area oscillations that represent the medial activity of actomyosin tend to produce more anisotropic cells, suggesting that different pools of contractility may produce distinct cell morphologies [79]. 3.2 Alternative cell contractility models The continuous actomyosin belt represented in the perimeter contractility model (Eq 4) produces uniform contractility along the cell boundary. In epithelial cells, however, the actomyosin displays distinct subcellular localizations that vary with confluence, functions, and developmental stages. Mature, confluent epithelial tissues predominantly exhibit a continuous, circumferential actomyosin structure along the cell perimeter, associated with adherens junctions [43]. In contrast, during early developmental stages, when confluency may be temporarily disrupted for large-scale tissue reorganization [87], actomyosin structures appear more medial, discontinuous, and heterogeneous—as has been observed in gastrulation of chick embryos [87] and Drosophila embryos [43,85]. Similar actomyosin structures have been observed in cell cultures of human intestinal [88] and Madin-Darby canine kidney (MDCK) cells [89]. These formations, often referred to as the medial pool of actomyosin, are essential for pulsed contractions and apical constriction [43,85]. To capture the aforementioned observations, some models incorporate non-uniform contractility (Fig 2C). Here, we consider two main types of these extensions, namely edge-based contractility and various medial contractile structures (summarized in Table 1). Firstly, tricellular junctions are thought to pose a physical obstacle to the transport of proteins along the cell perimeter, and thus the contractile force may not be uniform along the cell perimeter [90–93]. Experimental studies indicate that each bicellular junction functions as an independent contractile unit, where actin cables are anchored end-on to cadherin complexes binding to tricellulin of the tricellular tight junctions ensuring tissue confluency [94,95]. Vertex models with edge-based contractility terms aim to individually describe bicellular junctions of a cell [32,34,35,55]. In an edge-based contractility model, the last term of the vertex model energy function is modified (10) Here, in the third term, contractility is based on individual edges, with an elastic coefficient . Due to the quadratic expression, this energy function is different from the perimeter contractility model (Eq 4), even when the line tension coefficient kl and the contractility coefficient are constant for all edges. Analogous to the reformulation of the perimeter-based model (Eqs 4 and 8), one can absorb the line tension into the quadratic term in Eq 10 and introduce a target edge length l0. In doing so, the energy function takes a form, equivalent up to a constant: (11) In 2D apical vertex models embedded in 3D space, an anisotropic edge-based contractility induces chiral cell sliding driving asymmetric twisting similar to what is observed in organs like the heart [83,84]. A vertex model with both perimeter contractility and edge contractility can capture mixed effects of contractile actomyosin at the cell border [80]. Further, some researchers incorporate additional contractile structures into the model to represent force-generating mechanisms that are not confined to cell perimeters, such as the star-like formations recently observed on the basal side of epithelial cells [96]. When studying the oscillations of cell areas in the process of dorsal closure during Drosophila development, a series of studies add radial stress fibers, termed spokes, connecting the vertices with the center of each cell [32,35,45,81]. The spokes are implemented as contractile, elastic intracellular edges. In addition to representing the medial pool of actomyosin, they also play a similar role to pressure forces due to area elasticity. These stress fibers are still isotropic in the sense that they are oriented along different directions well distributed in space provided that cell shapes are not strongly elongated. Some other studies explore the effects of anisotropic contractile structures. Parallel stress fibers aligned with imposed unidirectional stretching are employed in an effective force-based vertex model to explore the scaling of the number of fibers with cell apical areas [82]. Similarly, the dynamics of various subcellular stress fiber geometries in a stress-fiber reinforced vertex model reveal a size-dependent response of cells [97]. One comparative study confirms that edge-based contractility models belong to the same sub-class of vertex models as perimeter-based contractility models, exhibiting a crossover from soft to stiff network behaviors [42]. Other studies suggest that different contractility modes and structures might influence cell packing geometry and the threshold of a rigidity transition [79,82]. Another way to incorporate anisotropic contractility from the medial myosin pool is to write an anisotropic continuous stress in the bulk of the cells and derive their effective forces on the cell-cell junctions and vertices [99]. This type of model has been used to investigate how epithelial colonies elongate through collective effects, with tissue elongation resulting from anisotropy in the average cell elongation [99], and how the V-shaped muscle compartments emerge in zebrafish, revealing that spatially modulating the mechanical environment around and within tissues can lead to complex organ shape formation []. In these two examples, the anisotropic bulk stress is crucial for symmetry breaking during tissue growth and sculpting tissue shapes, respectively. However, recent comparative work has shown that the two different formulations that have been employed to derive the forces on vertices [99], lead to distinct tissue dynamics due to differences in their treatment of non-affine deformations [100]. Moreover, neither description of anisotropic bulk stresses could fully represent equivalent forces arising from cell-cell junctions in the reference vertex model [100]. Overall, these results suggest that different modeling choices could lead to distinct predictions, and it is important to align model design with experimental observations of contractile structures, which vary across biological contexts. 3.3 Subcellular regulation of cell contractility and mechanosensitive response The reference model [23] describes cell-cell junctions as elastic springs driven to equilibrium by forces derived from an energy function Eq 4. In other words, the mechanical forces along the edges in the model are determined by the deformation of the polygonal tiling through a positive linear relationship. This simplified representation of the mechanical behavior of cell-cell junctions is hardly sufficient considering that the junctions are tightly regulated by complex, interacting subcellular elements such as cytoskeleton dynamics, mechanotransductive and other signaling pathways, and gene expressions. Two early seminal studies, using a lateral 2D vertex model, explore how deformation-dependent edge target length could lead to apical constriction and epithelial invagination [48,101]. The showcased mechanical self-organization suggests that tissues might leverage mechanical signals such as force or deformation to actively regulate their contractile forces. This type of feedback loop can result in spontaneous tissue configuration changes, a behavior loosely referred to as mechanosensitivity. With the advance in experimental and imaging techniques, knowledge of how cells actively control their mechanics on the subcellular scale provides increasingly solid foundations to the mechanosensitive hypothesis (Fig 3B) [11,82,91,102–110]. Together with the accumulation in experimental evidence, theoretical interest in mechanical self-organization has also grown rapidly [111]. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. The effective relation between forces and cell shapes in the vertex model can be changed by introducing submodules featuring active regulation and mechanosensitive feedback describing activities at the subcellular scale. (A) The mechanics in the reference vertex model [23] are governed by an energy function (a phenomenological Hamiltonian) that dictates a linear relation between forces and cell shape (observed e.g., in the ground states or a proliferation-free simulation). To be able to capture realistic tissue morphologies, vertex models are commonly extended to a full simulation framework with components refining the subcellular scales, resulting in effectively nonlinear relations between forces and cell geometry [65]. The subcellular activities include the activation and turnover of cytoskeleton machinery and adhesion proteins, the relaying of mechanical signals through mechanotransduction molecules, the modulation of mechanical forces by the spatiotemporal dynamics of signaling pathways, and the activation and deactivation of genetic programs instructing protein synthesis. (B) Examples of submodules that actively regulate quantities are strain-dependent remodeling of the target edge length, edge length-dependent active tension, or strain-dependent remodeling of the density of a biomolecule, e.g., myosin. In systems with mechanosensitive feedback, these remodeled quantities are usually interdependent and formulated as coupled differential equations. https://doi.org/10.1371/journal.pcbi.1012993.g003 The vertex model is a convenient tool to incorporate junction regulatory effects (Fig 3A). Cell division, as implemented in the reference model [23], is already an active event that reflects the cell cycle, an effect beyond the passive linear junctional tension dictated by the energy function, although it does not form an explicit feedback loop. One of the earlier studies combining vertex modeling and experiments on how mechanosensitivity contributes to morphogenesis encodes a feedback loop by letting the cell target parameters depend on the strain and suggests that strain-enhanced apical constriction is critical in the invagination of an optic-cup organoid [112]. Notably, the authors already deploy a 3D vertex model. This approach of implementing feedback loops in mechanics-dependent model parameters is a generalization of our first category of extensions with spatiotemporal dependence of model parameters. Meanwhile, it can be seen as a long-time limit of other implementations using differential equations to describe the evolution of model parameters in time, representing subcellular scale activities to study their effects at the tissue level. In more recent work, one approach, using the edge-based elasticity model (Eq 11), is to introduce remodeling of the target edge length l0,ab, described by the following generic formulation (12) where f is a general edge regulation function. We use the subscript ab for all variables to emphasize that they can be edge-specific. On the left-hand side of the equation is the absolute rate of change of the rest length [113], or alternatively, can also be the relative rate of change [34]. On the right-hand side, a variety of regulation mechanisms can be incorporated by adopting different forms of f. For example, the actomyosin cytoskeleton is well-known to be viscoelastic [114]. Thus, describes linear Maxwell viscoelasticity with viscosity . More importantly, the target edge length can be remodeled by the strain or strain rate [65,113], which is an active regulation mechanism. This type of feedback is proposed to explain the mechanical ratcheting of cell-cell junctions [113]. A recent variant of this model shows that strain-regulated tension, together with mechanical memory dissipation implemented through the fluctuation of tension coefficients and the relaxation of tension and strain, leads to transiently stable higher-fold vertices (rosettes) [65]. Notably, in contrast to other mechanosensitive models, this model incorporates negative feedback, where tension decreases under edge extension and increases under edge contraction. The negative feedback is found to stabilize rosettes [65]. Biochemistry can be more explicitly described by considering the dependence of f on the quantity or activity of a certain molecule, generically written as mab. In this case, an additional evolution equation is required to govern the dependent variable mab, which can be written in the form (13) In one example, mab denotes the average edge density of myosin motors, with Eq 12 describing the active walking of myosin motors contracting actin filaments under low tension and the slipping of motors under high tension, while the myosin density mab, in turn, depends on the strain rate of the edge through Eq 13 [34]. In this case, tissue viscoelasticity and cell size variability emerge from the feedback between the strain rate of cell-cell junctions and myosin recruitment, coupled with junctional viscoelastic relaxation [34]. In a 3D vertex model, the effect of remodeling several mechanical parameters – the surface area elasticity, the target surface area, and the surface tension – on apicobasal polarity is investigated [115]. By combining a vertex model with the biochemistry of key proteins, the study shows that the mechano-polarity feedback loop can drive spontaneous tissue folding. Besides remodeling the target edge length or the target cell perimeter, an additional or alternative scheme is to augment the edge tension with an active force. With the decomposition into cell pressure and edge tension (Eq 5), one can extend the edge tension as (14) where the additional active force is regulated by phenomenological laws describing subcellular processes. Similar to f in the edge target length remodeling formulation, can depend on an array of factors. For example, an additional contractile force can be added, which becomes active only when the cell area exceeds a threshold. This feedback enhances tissue strength and resistance to rupture [116]. Others assume that the active force is dependent on additional quantities, such as a myosin activity level (e.g., [35]). Persistent cell neighbor exchanges and convergent extension could arise from this type of mechanochemical feedback [32,45,117]. The active force term can also depend on a cell polarity vector. The coupling results in a chiral patterning when the tissue is confined on a spherical surface, demonstrating how mechanosensitivity interacts with topology [118]. In these cases, again new rate equations need to be introduced to govern the additional variables, either itself or the biochemical species mab. In general, both positive and negative feedback loops can be easily incorporated into the outlined framework. Some researchers consider further complexities in the regulation of cell mechanics. For example, the coupling of a lateral inhibition model with a vertex model reveals how mechanics and biochemistry combine to instruct the alternating cellular pattern in the inner ear [119]. With further progress in imaging technology and discovering molecular processes underlying the transmission and interpretation of mechanical signals by biological systems, we believe that with these extensions (summarized in Table 2) we are only seeing the beginning of what might become a major part of vertex model studies. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 2. Overview of vertex model extensions (Part 2). https://doi.org/10.1371/journal.pcbi.1012993.t002 4 Properties of extended vertex models With emerging model variations, an important question is how different models preserve certain physical properties compared to the reference model [23]. In this section, we discuss extended vertex models, focusing on tissue properties, including extensions that alter these properties through means other than subcellular dynamics, structures, and feedback. One property of significance is tissue rigidity. As reviewed in Sect 2, it is well known that the reference model displays a rigidity transition between fluid and solid characteristics with changing line tension and perimeter contractility (kl and kP in Eq 4, respectively) [23,29]. This transition has been found to occur at a critical target cell shape (Fig 1E) [29,53]. Imaging data from tissues without directed stresses reveal that the cell rearrangement rate undergoes a sharp transition at a critical average cell shape . This threshold corresponds to the critical target cell shape in the vertex model [10]. In tissues with a source of anisotropy, however, the cell shape alone does not suffice to describe or control rigidity. Under anisotropic conditions it is possible to obtain ordered, rigid states with elongated hexagonal cells beyond the threshold where a regular hexagonal lattice in an isotropic tissue loses order and rigidity. As a result, rigidity must be evaluated taking into account other quantities, e.g., cell shape alignment [10]. Several other observables and control parameters have a profound impact on the solid-to-fluid phase transition previously characterized by the target cell shape alone, resulting in higher-order transitions (Fig 4A,4B and 4D) [10,118,120]. Additionally, 2D vertex sheets restricted on spherical surfaces, as is the case in many realistic biological scenarios such as embryos and organoids, lead to tissue fluidity becoming a function of Gaussian curvature (Fig 4C) [46,121,122]. A similar effect is seen in a vertex model where edges take the shape of circular arcs based on the Young-Laplace law [123]. Here, the solid-to-fluid transition already occurs at a less elongated target cell shape of p0 = 3.73. A recent review highlights different physical mechanisms that cause the collective arrest of cell motion in tissues [30]. While vertex models have been primarily devoted to studying tension-driven rigidity so far, “freezing” a tissue by reducing parameter fluctuations and cell motility is another possible mechanism of rigidification in systems featuring a source of activity, e.g., a cell cycle or traction forces [63,118]. Disentangling contributions to rigidity is not trivial, especially in active vertex models where both mechanisms are at play [30]. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Dynamical phase diagrams describing a solid to fluid transition in various extended vertex models. (A) Average vertex coordination number (from to ) against target cell shape p0 in a model with active T1 transitions [120]; (B) alignment of cell shape against average measured cell shape in an anisotropic tension model [10]; (C) Gaussian curvature against target cell shape in a vertex model restricted to a sphere (positive curvature) or hyperbolic disk (negative curvature) [121,122]; (D) alignment rate of cell polarity against traction forces in a traction force model on a spherical surface [118]. Asterisks indicate the critical target cell shape or average measured cell shape ) in the reference model [29,39]. https://doi.org/10.1371/journal.pcbi.1012993.g004 The reference vertex model [23] and extended vertex models are used to explore cell packing geometry and reveal mechanisms for the characteristic stability of rosettes – vertices which belong to five or more cells [28]. Packings and rosettes are observable in experiments and are closely related to tissue rigidity [29,120]. For cell packing geometry in general, one way to generate disorder in the vertex model is to enable programmed growth and division on randomly selected cells [23]. Recent studies demonstrate that tissues with realistic cell packing can emerge from models without stochastic cell proliferation [36,65]. Instead, these models incorporate stochastic tension coefficients along with either cell shrinkage or mechanosensitive feedback between edge tension and target edge length. This indicates that gradual deterministic changes with a source of stochasticity appear to be a common motif in causing tissue disorder. On the other hand, there is evidence that different contractility models can lead to distinct cell morphologies [79]. At the same time, cell packing disorder is closely related to tissue rigidity by modifying its transition threshold [10,44]. Therefore, it remains to be explored how specific model extensions that have emerged recently influence cell packing and their implications for tissue rigidity. Among various cell packing features, four-fold vertices and higher-fold rosette vertices are of particular interest. The average vertex coordination number , summed over all V vertices, measures the density of rosettes, and in general cell packing disorder. The presence of higher-fold vertices, similar to isotropic cell shapes, is a good indicator for a rigid tissue [10,120]. A hexagonal lattice has , and the higher is, the higher the critical target cell shape beyond which the tissue becomes fluid (Fig 4A) [120]. Rosettes can be identified from experimental images relatively easily, thus providing another observable for tissue rigidity. In regular T1 transitions in the vertex model, the four-fold vertex state is energetically unstable and thus short-lived (Sect 2 and Fig 1D). In epithelial tissues, however, four-fold vertices and rosettes are frequently observed. The unexpected stability of higher-fold vertices is a subject of study where multiple of the reviewed extensions have explanation power. High rosette stability can be achieved by low amplitude, high-frequency fluctuations of edge tensions such that edges are virtually trapped inside a T1 transition [37,64]. Similarly, delayed T1 transitions, a consequence of the finite time scales required by molecular complexes to undergo rearrangements, increase the frequency and stability of rosettes and make tissues more elastic [66]. Increased stability of rosettes can also be linked to tension fluctuations arising from mechanosensitive feedback between tension and target edge length, independent of implementation details of T1 transitions [65]. These studies offer insights into how the frequency and stability of rosettes synergize with model extensions including fluctuations and T1 delay time. They constitute the beginning steps for investigating the relationship between higher-fold vertices and more sophisticated extensions such as non-uniform contractility and mechanosensitive models. The developments centered around tissue rheology and geometry imply that extensions of the vertex model can profoundly change predicted tissue properties, highlighting the need for comparative studies and unified views. Here, we identify three major lines of research. The first one continues along the subject of rigidity. A major result is the classification of the rigidity transition in the reference model [23,29] as one particular case of a broad class of rigidity arising from geometric incompatibilities that can be described by a minimal length [47]. Whether different vertex model variants display a rigidity transition has only been investigated recently by comparison of shear moduli in the quasi-static limit of vanished shear [42]. The authors have concluded that (i) the confluent foam model without the quadratic contractility term (see Sect 2) is unstable in the regime where the vertex model is in a fluid state; (ii) vertex models with edge contractility (Eq 10) [32,34,35,45,81] behave similarly to perimeter-based models (Eq 4) for small values of the area compressibility kA in displaying a transition from fluid-like to solid-like when l0 is decreased. They have also suggested the uncertainty of a quasi-static limit for a model with thresholded yielding and tension remodeling [113], which likely requires a stationary state that is constructed differently. A recent study generalizes the theory of elastic properties of athermal, under-constrained systems [47] to include thermal fluctuations in an edge-based contractility model without area elasticity (kA = 0). In thermal, under-constrained systems, it shows that energetic rigidity (also present in athermal systems) and entropic rigidity (the limit of infinitely stiff springs) “add up” like being connected in series [124]. This novel view unifies the physics of vertex models with that of polymer networks. The second line, which we believe holds potential for further exploration, is the influence of topological transition rules and more generally the effects of sources of activity. The effects of sharp thresholds and their potential effects and artifacts on tissue dynamics have been addressed in a small number of studies [39,52,66]. Different algorithms of topological transition have been compared in [39]. Here, a simple T1 transition algorithm (detailed in Sect 2) is compared with more complex T1 procedures, where transitions are selected based on the most energetically favorable changes in the network. While similar in polygon distributions and cell area variations, energy-based algorithms are less numerically stable because they depend on very fine changes in energy landscapes. Additionally, the study reveals that omitting quasistatic cell growth prior to applying a cell division rule (outlined in Sect 2) does not affect results significantly. A study comparing a large range of implementation choices finds that simulation outcomes are robust to variations in T1 edge lengths, T2 transition thresholds, and simulation step order [52]. However, the study also reveals that the cell cycle duration and the interplay of growth, division, and rearrangements influence the tissue size and mechanics. Further, the time step size influences rearrangement dynamics and higher-order methods do not improve convergence. As topological transitions affect rosette stability and tissue rigidity [37,64,120], simulation timing choices can significantly impact tissue behavior [52,66]. From a broader perspective, topological update rules belong to a class of elements that drive tissues out of thermal equilibrium, thereby making them active. Examples of such elements include a growth and division rule, biochemical circuits, and mechanical feedback loops. In the long term, it will be instructive to categorize active elements and clarify their effects on model behavior. So far, some authors have emphasized “active” components, e.g., an active vertex model with cell self-propulsion [125] or actively oriented cell divisions [126], while others have not stated so explicitly even though their models include similar elements such as cell division which in itself can be considered an active process [23]. A clear and unified terminology distinguishing active and passive elements could pave the way toward a better understanding of the potential effects of different model components. The third line, where further research can bring a significant understanding of tissue dynamics, is how different time scales interact. Many studies mentioned above already explore temporal effects [64,66,79]. In one recent comparative study, different friction coefficients used in the vertex model lead to distinct results in cell growth [127]. This is because the time scale of viscous relaxation arising from the friction with the ECM polymer substrate competes with intrinsic cellular time scales such as the cell cycle length and duration of growth. For models incorporating more complicated dissipation mechanisms such as tensor friction coefficients [128–130], there could be multiple viscous time scales to consider. Other models include internal dissipation from, e.g., area, edge and perimeter dissipation [128], neighbor-neighbor dissipation [130], and vertex-center viscosity [131]. However, as highlighted in an earlier review [27], dissipation processes in the vertex model and tissue mechanics remain poorly understood, and for vertex models employing over-damped motion as a numerical tool to solve energy minimization problems, caution is needed to ensure artificial viscous time scales do not distort results [27,52,127]. For models without a viscous time, other time scales, particularly those introduced by mechanosensitive feedback (Eqs 12 and 13), may be relevant. Future investigations explicitly addressing the interaction of such time scales could provide key insights into tissue dynamics. Examining the properties of tissues raises the question of whether a unified theoretical framework for the vertex model can be constructed – one that could encompass simplified versions under different conditions and clarify how these simplified versions correspond to existing model parameters. Could a universal theory for vertex models emerge, akin to the ideal gas law in thermodynamics or the Navier-Stokes equations for fluid dynamics? The phenomenological Hamiltonian provides a promising foundation and has already yielded significant insights, particularly in understanding tissue rigidity. However, systematic, comparative studies remain limited, though they are essential to furthering both the practical applications of vertex models and our understanding of tissue mechanics. 5 Inferring mechanical stresses and properties in cells and tissues Inference methods using statistical models to estimate mechanical stress in tissues are an orthogonal approach to mechanistic vertex models. The vertex model and other models discussed later in this review, e.g., the Voronoi model, can be classified as mechanistic models. Mechanistic models are designed to capture hypothesized mechanisms that are able to explain experimental observations. Statistical methods, on the other hand, aim for a statistical model that is able to capture the structure of observational data. In a confluent network of cells, this translates to a statistical method that, based on cell geometries extracted from cell boundary segmentations from experimental data, infers effective cell boundary tensions and cell pressures (Eq 5). This method assumes local force balance in the bicellular junction network, similar to the vertex model (Eq 1) [49,50,132–134]. This makes it possible to compare predictions by vertex models based on specific phenomenological energy functions (e.g., Eq 4B) to estimates of effective forces (Fig 1C). A recent study includes inference methods into their model design process by first estimating how stresses in a tissue relate to strains [41]. By an information criterion, they find that in Drosophila pupal wing, pupal notum, and embryo an anisotropic edge-based contractility model (Eq 10) describes the relation between forces and cell shape better than the conventional perimeter-based model (Eq 4). They find that positive feedback between junction tension and shrinkage increases at specific developmental stages. Future studies will profit from including inference methods into their model design process by first estimating stresses in a tissue to characterize nonlinearities in the relation between stress and strain (Fig 3). Moreover, because stress inference is non-invasive, it can be used to estimate local stresses from a time-lapse dataset [49,132]. Mapping stress estimates in space and time enables comparison to vertex model simulations based on effective forces on a vertex, namely the changing forces that a specific tricellular junction is exposed to over a course of minutes or hours. In general, estimates of stresses during topological transitions could be helpful to profile energy barriers that are otherwise intractable with commonly used destructive laser ablation measurements. Details of stress inference methods and applications have been discussed in a previous review [51]. Furthermore, a novel approach aims to infer flow properties of tissues from their geometry and dynamics, based on vertex model simulations [40]. The method relies on the ratio of the strain rate increment and the stress increment to distinguish between regimes. Comparing results from conventional vertex model simulations in the solid phase with Drosophila pupal wing data, the study finds that the distribution of bond lengths can serve as a proxy for local distances to yield stress and thus reveal the regime in which flow is occurring [40]. This represents a promising direction of research on inference of tissue properties by leveraging a deeper understanding of the vertex model. 6 Complementary models for epithelial tissues Here, we provide a short overview of several other types of models that are related to and complement the vertex model by providing alternative approaches to tissue mechanics at different spatial resolutions, degrees of freedom, and levels of details in the mechanics. Comparison with related models provides valuable insights into the strengths and limitations of the vertex model in resolving large system sizes, adhesion dynamics, and cell shapes. 6.1 Continuum models capturing vertex model properties At macroscopic scales, continuum models, coarse-grained from vertex models, provide a possible means to bridge cellular-level mechanics with large-scale tissue behavior, without resolution of single cells. They formulate partial differential equations governing the evolution of several field functions, e.g., displacements, stress, strain, and strain rate, and are thus amenable to a variety of analytical and computational tools. In this way, tissues at macroscopic scales consisting of far more cells than in a vertex model can be studied. There have been many attempts at deriving the continuum limits of various elements of the vertex model. An energy function was formulated in terms of the strain tensor, and a solid-solid transition due to geometric frustration has been uncovered [135]. A normal mode analysis can be applied to the vertex model to compute the linear loss and storage moduli, revealing complex linear viscoelastic rheology of a 2D tissue [129,136]. Full dynamical theories could be derived by directly postulating the equations of elasticity or hydrodynamics and fitting at least part of the parameters to vertex model simulations or experiments [61,137,138]. They can also be constructed by homogenization techniques [77] or thermodynamic formalisms including Onsager’s approach [139], Poisson brackets [140] and a kinetic theory formulating a Fokker-Planck equation [141]. Such procedures can often give relations between macroscopic material parameters and vertex model parameters. A detailed survey of these studies and their predictions compared to each other and those from the vertex model are beyond the scope of this review. We point out, however, that these models focus on coarse-graining the reference model [23], and little is known about the continuum limits of extended vertex models like the ones discussed in Sect 3. 6.2 Vertex-like models with different resolutions of cell-cell junctions Several models can be considered “vertex-like” models in that they use a vertex model energy or strongly focus on resolving cell-cell junctions and topological transitions (Fig 5). The Voronoi model is a coarse-grained version of the vertex model with cell centers as degrees of freedom. It relies on the duality between the Delaunay triangulation of cell centers and the Voronoi tesselation, which resembles a membrane network, to update the equations of motion of the cell centers based on a vertex model energy function (Eq 8). It offers lower implementation costs due to the implicit treatment of topological transitions at similar computational cost ( with N degrees of freedom) [125]. However, in the 2D case, the Voronoi model lacks a solid-to-fluid transition due to constraints imposed on cell areas and perimeters [125,142–144]. Similarly, a comparative study between 3D vertex models and 3D Voronoi models suggests broader distributions of cell shapes in vertex models, highlighting the power of vertex models in capturing heterogeneity in tissue properties [144]. Although this restricts the benefit of Voronoi models to mostly solid 3D systems, in fact many systems of interest in biology belong to this category [27,118,145–147]. At finer scales, models like the N-vertices models [7,13,148–150], the deformable polygon model (DPM) [151], and the apposed cortex adhesion model (ACAM) [152] allow for curvature, slippage and delamination of cortices, offering a more realistic description of cell-cell contacts. Recent experimental findings support the notion of a compression-dominated solid, which is observed in DPM simulations as a tension-to-compression transition between two solid regimes [151,153]. In general, further investigation is needed to distinguish between different phase transitions, described by various terms, such as epithelial-to-mesenchymal transition (EMT), jamming, rigidity, or percolation transition. Many of the observed phase transitions in epithelial tissues are not exclusively driven by only one of the commonly discussed mechanisms: cortical tension, cell-cell adhesion, tissue density or cell motility. Disentangling multiple mechanisms is of particular interest in, e.g., tissues undergoing embryonic development or cancer progression [30,153–155]. Multi-phase field models and continuum approximations of cell membranes provide avenues for exploring topological transitions, energy landscapes and continuous cell shape changes in epithelial tissues [156–160]. Moving forward, comparing and coupling model extensions and variants at different length scales and benchmarking them against material property studies in other models, such as the cellular Potts model, can offer deeper insights into tissue mechanics and behavior [161–164]. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 5. Vertex-like models by resolution of cell-cell junction dynamics. Increasing resolution of cell-cell junctions from left to right: The Voronoi model is a coarse-grained version of the vertex model with cell centers as degrees of freedom and implicit treatment of T1 transitions [125,142,143]; the vertex model as reference [28]; the deformable polygon or N-vertices model resolves the cell cortex by expanding the definition of vertices to a subjunctional flexible chain with explicit adhesion between cell membranes [7,13,148–151]; the apposed cortex adhesion model focuses on adhesion dynamics of cell-cell junctions during topological transitions with viscoelastic continuum cortices adhering to each other by an explicit continuous adhesion energy [152]. https://doi.org/10.1371/journal.pcbi.1012993.g005 6.1 Continuum models capturing vertex model properties At macroscopic scales, continuum models, coarse-grained from vertex models, provide a possible means to bridge cellular-level mechanics with large-scale tissue behavior, without resolution of single cells. They formulate partial differential equations governing the evolution of several field functions, e.g., displacements, stress, strain, and strain rate, and are thus amenable to a variety of analytical and computational tools. In this way, tissues at macroscopic scales consisting of far more cells than in a vertex model can be studied. There have been many attempts at deriving the continuum limits of various elements of the vertex model. An energy function was formulated in terms of the strain tensor, and a solid-solid transition due to geometric frustration has been uncovered [135]. A normal mode analysis can be applied to the vertex model to compute the linear loss and storage moduli, revealing complex linear viscoelastic rheology of a 2D tissue [129,136]. Full dynamical theories could be derived by directly postulating the equations of elasticity or hydrodynamics and fitting at least part of the parameters to vertex model simulations or experiments [61,137,138]. They can also be constructed by homogenization techniques [77] or thermodynamic formalisms including Onsager’s approach [139], Poisson brackets [140] and a kinetic theory formulating a Fokker-Planck equation [141]. Such procedures can often give relations between macroscopic material parameters and vertex model parameters. A detailed survey of these studies and their predictions compared to each other and those from the vertex model are beyond the scope of this review. We point out, however, that these models focus on coarse-graining the reference model [23], and little is known about the continuum limits of extended vertex models like the ones discussed in Sect 3. 6.2 Vertex-like models with different resolutions of cell-cell junctions Several models can be considered “vertex-like” models in that they use a vertex model energy or strongly focus on resolving cell-cell junctions and topological transitions (Fig 5). The Voronoi model is a coarse-grained version of the vertex model with cell centers as degrees of freedom. It relies on the duality between the Delaunay triangulation of cell centers and the Voronoi tesselation, which resembles a membrane network, to update the equations of motion of the cell centers based on a vertex model energy function (Eq 8). It offers lower implementation costs due to the implicit treatment of topological transitions at similar computational cost ( with N degrees of freedom) [125]. However, in the 2D case, the Voronoi model lacks a solid-to-fluid transition due to constraints imposed on cell areas and perimeters [125,142–144]. Similarly, a comparative study between 3D vertex models and 3D Voronoi models suggests broader distributions of cell shapes in vertex models, highlighting the power of vertex models in capturing heterogeneity in tissue properties [144]. Although this restricts the benefit of Voronoi models to mostly solid 3D systems, in fact many systems of interest in biology belong to this category [27,118,145–147]. At finer scales, models like the N-vertices models [7,13,148–150], the deformable polygon model (DPM) [151], and the apposed cortex adhesion model (ACAM) [152] allow for curvature, slippage and delamination of cortices, offering a more realistic description of cell-cell contacts. Recent experimental findings support the notion of a compression-dominated solid, which is observed in DPM simulations as a tension-to-compression transition between two solid regimes [151,153]. In general, further investigation is needed to distinguish between different phase transitions, described by various terms, such as epithelial-to-mesenchymal transition (EMT), jamming, rigidity, or percolation transition. Many of the observed phase transitions in epithelial tissues are not exclusively driven by only one of the commonly discussed mechanisms: cortical tension, cell-cell adhesion, tissue density or cell motility. Disentangling multiple mechanisms is of particular interest in, e.g., tissues undergoing embryonic development or cancer progression [30,153–155]. Multi-phase field models and continuum approximations of cell membranes provide avenues for exploring topological transitions, energy landscapes and continuous cell shape changes in epithelial tissues [156–160]. Moving forward, comparing and coupling model extensions and variants at different length scales and benchmarking them against material property studies in other models, such as the cellular Potts model, can offer deeper insights into tissue mechanics and behavior [161–164]. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 5. Vertex-like models by resolution of cell-cell junction dynamics. Increasing resolution of cell-cell junctions from left to right: The Voronoi model is a coarse-grained version of the vertex model with cell centers as degrees of freedom and implicit treatment of T1 transitions [125,142,143]; the vertex model as reference [28]; the deformable polygon or N-vertices model resolves the cell cortex by expanding the definition of vertices to a subjunctional flexible chain with explicit adhesion between cell membranes [7,13,148–151]; the apposed cortex adhesion model focuses on adhesion dynamics of cell-cell junctions during topological transitions with viscoelastic continuum cortices adhering to each other by an explicit continuous adhesion energy [152]. https://doi.org/10.1371/journal.pcbi.1012993.g005 7 Software for vertex and vertex-like models As vertex models become increasingly complex and more popular with practitioners across disciplines – specifically mathematics, physics, engineering, biology, and health sciences – the development of accurate, extensible, scalable, and accessible computational libraries becomes crucial. A variety of software packages caters to different needs in a diverse community of modelers (summarized in Table 3). Surface Evolver is a standard solver that computes the geometry of any given surface by minimizing an associated energy function using gradient descent, accommodating arbitrary topology, volume constraints, boundary constraints, and boundary contact angles [165]. It continues to be used to solve the quasi-static equations of motion, especially for models based on the reference model (eqs 2 and 8) [29,39]. Cell-based Chaste is a submodule of the Cancer Heart and Tissue Environment, designed for comprehensive cell population modeling [166,167]. It integrates cell cycle, cell motion, and biomolecule transport models, supporting various frameworks, including lattice-free models (vertex, Voronoi, overlapping spheres, and immersed boundary models), and on-lattice models (cellular automata and cellular Potts models). SAMoS is a Delaunay-based Voronoi model leveraging Delaunay-Voronoi duality for efficiency, as it avoids recalculating the entire Voronoi tessellation at each time step [125]. cellGPU implements a GPU-accelerated vertex model that is approximately ten times faster than its CPU counterpart and a hybrid GPU/CPU Voronoi model offering a speed increase of up to three orders of magnitudes for tissues in the topological transition-free solid regime [168]. TiFoSi is a simulation tool for planar epithelial cellular dynamics, integrating mechanosensitive feedback, gene expression, and interactions between cell populations. It allows customizable cell cycles, cleavage properties, protein partitioning during division, and cell communication [169]. Tyssue is a Python package for 2D and 3D vertex models, offering a flexible API for rapid model development and a range of energy function terms [170]. It also supports 2D apical surfaces embedded in 3D as well as full 3D modeling. The Vertex Model Solver of Tissue Forge enables event-based modeling and integrates vertex models with particle-based and subcellular methods, offering real-time simulation visualization [171]. neoVM is an implementation of a graph vertex model that utilizes knowledge graphs and a non-redundant database to streamline topological transitions [172]. It executes graph transformations through predefined queries, simplifying implementation and data handling. Most of the available software packages are implemented under an object-oriented paradigm and provide generalized representations of energy functions and numerical solvers in a modular fashion that make it convenient to couple to stochastic models and ODEs that govern a vertex model parameter. The diversity of these libraries manifests the popularity of vertex models in studies of tissue mechanics. We expect their capabilities to be continuously enhanced and expanded. Comparative studies of vertex models and their fundamental properties, as envisioned in earlier sections, could contribute to more comprehensive and generalized software packages. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 3. Overview of software packages for vertex and vertex-like models. https://doi.org/10.1371/journal.pcbi.1012993.t003 8 Conclusion In this review, we have surveyed recent advancements in the vertex model. We find that extensions to a seminal vertex model [23] have largely aimed at coupling tissue mechanics with model constructs describing subcellular activities. These extensions are made possible through several types of implementations: local and temporal dependence of model para-meters, stochastic fluctuations, non-uniform or non-junctional contractility structures, and coupled biochemical and mechanical regulatory mechanisms. Vertex models have been successfully applied to obtain insight in both tissue properties in general and the mechanics of biological processes, specifically tissue dynamics, e.g., dorsal closure [35,36,76,78], convergent extension [6,10,32,33,35,45,66], and tumor invasion [36,44]. In synthesizing the predictions on tissue properties, we point out that different models could lead to similar results in certain aspects, such as the generation of disorder by either growth or shrinkage coupled to stochasticity [36,65]. In other aspects, however, distinct extensions can lead to significant differences in properties, such as cell morphologies and tissue rigidity transitions. This leads to questions regarding their comparability and the quest for a unified theoretical framework, which have just started to be addressed [6,41,42,47]. We also call for attention to several model components whose effects have not been sufficiently examined, including different time scales in the model and topological transitions as a source of activity. Statistical inference methods, based on a cellular network geometry and a force balance assumption, provide an orthogonal approach to estimating mechanical parameters from experimentally measured cell geometries, enabling data-driven model selection, force estimates, and model comparisons in epithelial tissues [41,49,51,134]. Furthermore, we present a comparison of several other models to the vertex model. Continuum models, coarse-grained from vertex models, enable the study of large-scale tissue behavior by bridging cellular mechanics and macroscopic properties using analytical and computational tools, but their applicability to extended vertex models remains largely unexplored. Models similar to the vertex model in resolving cell shapes include the Voronoi model, the N-vertices model, the deformable particle model, and the apposed cortex adhesion model [144,151,152]. These models differ from vertex models in how they resolve cell-cell junction dynamics, which affects their tissue-scale behavior and emphasizes the multiscale nature of tissue dynamics. As vertex models gain popularity, the development of accurate, extensible, scalable, and accessible computational libraries like Surface Evolver, Cell-based Chaste, SAMoS, CellGPU, TiFoSi, Tyssue, Tissue Forge, and NeoVM has become crucial to facilitate simulations of biological processes across disciplines. With further collaborative effort across experimental, computational, and theoretical communities, we expect the vertex model to become a well-formalized and robust framework to study how the physical properties of tissues are influenced by kinetics and dynamics at the subcellular scale. It is evident that in recent years, there has been significant progress in extending vertex models and studying their properties. We expect that research in vertex models will continue advancing in multiple directions. First, the diversity of extensions will continue growing to incorporate additional mechanical and biochemical processes now accessible by experimental approaches such as optogenetics [15]. Of particular interest could be adapting vertex models to further study different developmental stages [43] and more complex cell topologies which require sophisticated 3D models [145,147,172]. Second, comparative studies are essential to clarify the differences in vertex models equipped with different extensions. Through these studies, it could become clear which physical properties of tissues are conserved to form a well-defined class of vertex models and how properties are potentially lost [42]. One established aspect is tissue rigidity. On the other hand, topological transition rules and their effects remain underexplored and could benefit from systematic works similar to [6,65,66]. Third, the development of a multiscale approach in many recent vertex models presents opportunities to unravel the dynamics of epithelial tissues comprehensively and bridge molecular processes to the scale of organs. In particular, integrating submodules at smaller scales effectively supplies parameters at the cellular scale, e.g., tension coefficients [34,35]. Four, more extensive comparisons between model predictions and experimental data would enhance model predictability and fidelity, allowing for deeper insights into tissue dynamics, and potentially illuminating how tissue properties could be engineered for medical applications. We believe the vertex model framework holds promise for advancing our understanding of epithelial tissue mechanisms, morphodynamics, and properties by integrating subcellular structure, dynamics, and feedback. Acknowledgments MQ wishes to thank Francis Corson and Vincent Hakim for stimulating discussions. The authors would like to thank Clinton Durney for feedback on both the scientific content and the clarity of the manuscript and Karin Lange, Lange Design Intelligence for offering technical expertise and advising on the final artwork of the figures and tables.
TI - Vertex models capturing subcellular scales in epithelial tissues
JF - PLoS Computational Biology
DO - 10.1371/journal.pcbi.1012993
DA - 2025-05-21
UR - https://www.deepdyve.com/lp/public-library-of-science-plos-journal/vertex-models-capturing-subcellular-scales-in-epithelial-tissues-CQFBPjdoOf
SP - e1012993
VL - 21
IS - 5
DP - DeepDyve
ER -