Advertisement

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Graph topological transformations in space-filling cell aggregates

  • Tanmoy Sarkar,

    Roles Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Writing – original draft, Writing – review & editing

    Affiliation Department of Theoretical Physics, Jožef Stefan Institute, Ljubljana, Slovenia

    ⨯
  • Matej Krajnc

    Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Methodology, Project administration, Software, Supervision, Validation, Writing – original draft, Writing – review & editing

    matej.krajnc@ijs.si

    Affiliation Department of Theoretical Physics, Jožef Stefan Institute, Ljubljana, Slovenia

    ⨯

Abstract

Cell rearrangements are fundamental mechanisms driving large-scale deformations of living tissues. In three-dimensional (3D) space-filling cell aggregates, cells rearrange through local topological transitions of the network of cell-cell interfaces, which is most conveniently described by the vertex model. Since these transitions are not yet mathematically properly formulated, the 3D vertex model is generally difficult to implement. The few existing implementations rely on highly customized and complex software-engineering solutions, which cannot be transparently delineated and are thus mostly non-reproducible. To solve this outstanding problem, we propose a reformulation of the vertex model. Our approach, called Graph Vertex Model (GVM), is based on storing the topology of the cell network into a knowledge graph with a particular data structure that allows performing cell-rearrangement events by simple graph transformations. Importantly, when these same transformations are applied to a two-dimensional (2D) polygonal cell aggregate, they reduce to a well-known T1 transition, thereby generalizing cell-rearrangements in 2D and 3D space-filling packings. This result suggests that the GVM’s graph data structure may be the most natural representation of cell aggregates and tissues. We also develop a Python package that implements GVM, relying on a graph-database-management framework Neo4j. We use this package to characterize an order-disorder transition in 3D cell aggregates, driven by active noise and we find aggregates undergoing efficient ordering close to the transition point. In all, our work showcases knowledge graphs as particularly suitable data models for structured storage, analysis, and manipulation of tissue data.

Author summary

Space-filling polygonal and polyhedral packings have been studied as physical models for foams and living tissues for decades. One of the main challenges in the field is to mathematically describe complex topological transformations of the network of cell-cell interfaces that are present during cell rearrangements, accompanying plastic deformations and large-scale cellular flows. Our work addresses this challenge by storing the topology of the network of cell-cell interfaces into a knowledge graph with a specific data structure, uniquely defined by a metagraph. It turns out that this graph technology, also used by tech giants such as Google and Amazon, allows representing topological transformations as graph transformations, that are intuitive, easy to visualize, and straight-forward to implement computationally.

Citation: Sarkar T, Krajnc M (2024) Graph topological transformations in space-filling cell aggregates. PLoS Comput Biol 20(5): e1012089. https://doi.org/10.1371/journal.pcbi.1012089

Editor: Changbong Hyeon, Korea Institute for Advanced Study, Republic of Korea

Received: November 24, 2023; Accepted: April 19, 2024; Published: May 14, 2024

Copyright: © 2024 Sarkar, Krajnc. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the paper and its Supporting information files. The source code of neoVM package can be found at https://gitlab.com/ijskrajncgroup1/neovm without access restrictions.

Funding: The authors acknowledge financial support from the Slovenian Research Agency [research project No. J1-3009 (TS and MK) and research core funding No. P1-0055 (MK)]. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

The mechanical interactions between individual cells and between the cells and their environment play a crucial role in determining the macroscopic properties and behaviors of animal tissues on a large scale [1–7]. During embryonic development, for example, forces generated within cellular cortices drive precise and highly orchestrated active deformations and collective cellular flows [8–13]. The mechanical forces are transmitted across the tissue through cell-cell contacts, which form a complex spatial network with a dynamically changing topology [14–25]. Apart from better understanding biological processes such as development, regeneration and disease [14, 26–28], studying the mechanical properties of cell aggregates also facilitates the development of biomimetic materials and devices for applications in medicine, engineering, and materials science [6, 29].

Tissue-scale mechanics have been addressed by computational models that represent cells as discrete entities with certain physical properties that phenomenologically describe the mechanics at smaller length scales [30–39]. Recent studies compare various widely used computational models for simulating cell assemblies including cellular automaton, cellular Potts model, phase-field model, particle-based model, deformable-particle model, Voronoi model, and vertex model [40–42]. While many of these approaches naturally incorporate cell rearrangements–processes that are crucial for capturing a realistic rheological tissue response–handling these processes computationally in the vertex model is quite challenging. In particular, the vertex model describes individual cells as polygons (2D) or polyhedra (3D) and parametrizes their shapes by vertex positions [43–46]. Cell-cell contacts are represented by edges in 2D and polygons in 3D and comprise an intricate network whose connectivity needs to be computationally altered upon every cell-rearrangement event. Current state-of-the-art computational tools for vertex-model simulations offer certain solutions to simplify the implementation of these network-reconnection events, however, these methods are not easily generalizable to 3D [25, 39].

Generally, vertex models store the topology and geometry of the tissue in a tabular form within arrays and perform topological transformations of the cell-interface network by updating these arrays according to specified rules [47–50]. Although programming the routines that perform these dynamic array updates is still relatively manageable for planar polygonal cell packings and even for 3D surface packings involving polyhedral prism-like cells [51–55], developing computer codes for simulations of 3D bulk cell aggregates poses a significantly greater challenge. Indeed, since the pioneering work by Honda et al. [43], who first introduced a vertex model of 3D cell aggregates, there have only been a few recent works reporting successful attempts of coding a full 3D vertex model with dynamic cell rearrangements [56–58]. The difficulty of implementing these rearrangements with the conventional data model raises questions about its suitability and challenges our basic understanding of rearrangements in space-filling packings.

Nevertheless, the vertex model often offers a profound mechanistic understanding of tissue behaviors, which is to a large extent facilitated precisely by the explicit presence of cell boundaries and their ability to remodel. Therefore, despite the challenges associated with cell-boundary tracking during cell rearrangements, the vertex model may still be in many respects advantageous over other models, in which cell boundaries need not be explicitly managed. For instance, the vertex model excels in realistically modeling large-scale three-dimensional tissue deformations and cellular flows, that may be driven by local active forces [43, 46, 59, 60] while, at the same time, capturing certain small-scale structural features, including single-cell shapes and their sidedeness within the aggregate.

To address the challenges associated with the implementation of cell rearrangements in the vertex model, we introduce a reformulation of the vertex model called Graph Vertex Model (GVM). We discover a particular graph-data model, which allows formulating topological transformations of cell networks as simple graph transformations. The blueprint outlining the relationships among the components is specified through a metagraph which is designed in a manner that topological transformations are themselves represented by graphs. This design not only enhances their intuitive and visual understanding but also simplifies their implementation, making it accessible even to researchers with limited programming expertise. Moreover, we demonstrate that within the GVM’s data representation, a T1 transition–the basic rearrangement event in 2D cellular tilings–emerges as a subgraph of graph transformations representing cell rearrangements in 3D packings. This allows developing generalized computational codes that are applicable to both 2D and 3D, suggessting that GVM’s data model may be the most natural representation of these systems. As a proof of concept, we develop an open-source Python package neoVM, which implements GVM over a graph database, managed in Neo4j [61, 62]. We use our new approach to study order-disorder transition in 3D cell aggregates. We characterize the transition and find aggregates undergoing most efficient ordering in the vicinity of the transition point.

Vertex model

The vertex model represents a cell aggregate by a three-dimensional packing of space-filling polyhedral cells (Fig 1A). Cell shapes are parametrized by positions of vertices in the (x, y, z) space: (1) Here i = 1, …, Nv, where Nv is the total number of vertices.

thumbnail
Download:
Fig 1. Geometry and topology of cell aggregates.

A Cell aggregate is modeled by a space-filling packing of polyhedral cells. Cell shapes are parametrized by vertex positions ri, which move according to mechanical forces Fi. Topology of the cell network is specified by lists ej, pk, and cl [Eqs (2)–(4)], which store head and tail vertices of edges, oriented edges within polygons, and oriented polygons within cells, respectively. B EV transition merges vertices of an edge into a single vertex, whereas VT transition resolves a vertex into a triangle. VE and TV transitions are inverse transitions of EV and VT transitions, respectively. C Polygons involved in topological transitions from panel B. D Decomposed schematic highlighting from 3 different perspectives of polygons, edges, and vertices, in topological transitions from panel B.

https://doi.org/10.1371/journal.pcbi.1012089.g001

Pairs of vertices are connected by oriented edges, defined by the indices of the constituent vertices: (2) Here and denote indices of head and tail vertices of edge j, respectively (Fig 1A), and j = 1, …, Ne, with Ne being the total number of edges. The signs of the vertex indices, denote head and tail vertices of the edge. While at this point these signs seem redundant, since the order in which the indices appear in ej itself indicates the head/tail role of the corresponding vertices, they will become important later on.

Polygonal cell sides are defined by oriented lists of indices of their constituent edges as (3) where k = 1, …, Np, with Np being the total number of polygons. The edges in the list pk are listed sequentially and denotes the orientation of the μ-th edge (with index ) within the polygon relative to a chosen positive direction (Fig 1A).

Similarly, cells are defined by oriented lists of indices of the constituent polygons as (4) where l = 1, …, Nc, with Nc being the total number of cells within the aggregate. Polygon orientations , where −1 and +1 correspond to the polygon’s normal vector pointing towards the cell center and away from the cell center, respectively (Fig 1A). The direction of the polygon’s normal vector is defined by the right-hand-screw rule, where the fingers curl along the chosen polygon’s positive orientation of the bounding edges (Fig 1A). In contrast to the case of polygons, pk, polygon indices in cl need not be listed in any specific order.

Exemplary structures of cell aggregates with Nc = 128 cells are given in the online repository of neoVM [61] and their generation is described in Methods.

The dynamics of cell-shape changes are described by simulating movements of the vertices, driven by mechanical forces. In a model that neglects inertial effects and only considers friction of vertices with a static background, vertices follow first-order dynamics described by (5) Here, η is the friction coefficient, W is the potential energy of the aggregate, and is a system-specific active-force contribution. The potential energy is typically calculated from geometric properties of cells such as surface areas of cell-cell contacts A, cell volumes V, etc. These quantities are calculated from the vertex positions as sums over geometric elements, i.e., vertices, edges, polygons (Methods).

Topological transitions

In addition to changing their shapes, cells also change their relative organization within the tissue by exchanging their neighbors [43]. These reorganization events alter the topology of the network of cell-cell contacts, thereby also affecting the interaction between the vertices. In confluent cell aggregates, cells exchange their neighbors by (i) merging vertex pairs of vanishingly short edges and (ii) resolving these vertices into new edges [63–65]. The topology of the edge network after these transformations locally differs from the initial one. In particular, cells that were initially separated might become neighbors, whereas pairs of initially neigboring cells may separate.

To model cell rearrangements in 3D cell aggregates, the following elementary local topological trasformations need to be considered (Fig 1B–1D): (i) Edge-to-vertex (EV) transition merges vertices of a vanishingly short edge into a single 6-fold vertex, (ii) vertex-to-triangle (VT) transition resolves a 6-fold vertex into a new triangle, (iii) triangle-to-vertex (TV) transition merges all three vertices of a triangular polygon into a single 6-fold vertex, and (iv) vertex-to-edge (VE) transition resolves a 6-fold vertex into a new edge.

An EV followed by a VT completes an ET transition, which transforms a vanishing edge into a triangle, formed in the perpendicular direction to the shrinking edge (Fig 1B and 1C). After an ET transition, initially separated cells become neighbors by sharing the new triangle. Similarly, a TV followed by a VE, completes a TE transition, which transforms a vanishing triangle into an edge, formed in the perpendicular direction to the shrinking triangle (Fig 1B and 1C). After a TE transition, the neighboring cells initially sharing the triangle become separated. We note that in a special case of an epithelial monolayer, where cells only attach to their neighbors laterally but not apically and basally, EV transitions either on basal or apical edges give rise to scutoidal cells [66, 67].

Importantly, as illustrated in Fig 1D, the basic topological transformations in fact resemble multiple edge-to-vertex and vertex-to-edge transitions, known in 2D polygonal networks as T1 transitions. In a T1 transition, a pair of initially neigboring polygons becomes separated, whereas polygons from the remaining (initially separated) polygon pair become neighbors. This similarity between 2D and 3D transformations suggests that it may be possible to unify topological transformations in 2D and 3D space-filling packings, provided that the vertex model’s core architecture be properly reformulated.

The main issue with the conventional implementation of the vertex model is that vertices, edges, polygons, and cells are stored in separate arrays (ri, ej, pk, and cl, respectively), which do not directly encode any interconnections or relationships among their respective elements. In particular, elements that might be spatially and topologically related are generally not stored together in the database and accessing any high-level topology data (e.g., finding cells that share a common vertex) requires inefficient searches over all the elements of the cell network. To avoid these inefficient searches, the conventional vertex models store data in a highly redundant form, where higher-level information about the topology of the cell network are stored in addition to ej, pk, and cl (even though these higher-level information may be calculable from ej, pk, and cl). For instance, to efficiently search for cells that share a certain vertex, lists of cells sharing a common vertex need to be stored for all vertices. Indeed, retreiving this information from ej, pk, and cl on the fly would require highly inefficient looping over all the cells. Due to this data redundancy, algorithms that manipulate the data arrays upon topological transformations in a self-consistent manner are difficult to program.

Results

Knowledge graph

To overcome challenges related with implementing topological transformations in the vertex model, we propose a new approach, based on storing the topology of the cell network into a knowledge graph, which uses a graph- rather than a tabular data model. By construction, the elements that are topologically related are also connected in the knowledge graph and therefore, any high-level information about the topology of the cell network is readily retrievable by querying over the relevant part of the database with no need of storing any redundant data.

Knowledge graph is a graph data structure, which represents a network of real-world entities and relationships between them [68]. These data are stored in a graph database where entities and relationships are represented by nodes and links, respectively, and can, additionally, carry multiple properties. For example, a movie database can be stored as a knowledge graph, in which the data about actors and directors for a given movie are represented by nodes labeled Person and Movie and relationships labeled ACTED_IN and DIRECTED. In such a knowledge graph, the information that Cillian Murphy acted in the movie Oppenheimer, directed by Christopher Nolan, can be stored as (p1:Person {name: “Cillian Murphy”})-[:ACTED_IN]->(m:Movie {title: “Oppenheimer”})<-[:DIRECTED]-(p2:Person {name: “Christopher Nolan”}); here nodes labeled Person and Movie carry the name and the title properties, respectively.

The notation used here and in the following sections follows the syntax of the Cypher graph query language [69]. It is important to note that GVM relies on general principles of discrete mathematics and does not depend neither on the choice of the query language nor on the choice of the database-management framework. Computationally implementing GVM, however, does require some basic knowledge of graph databases and query languages.

Metagraph

Real-world entities in GVM are vertices, edges, polygons, and cells and they are stored in a knowledge graph in a hierarchical manner as nodes labeled Vertex, Edge, Polygon, and Cell, respectively. The topology of the cell network is encoded through relationships labeled IS_PART_OF. These relationships are directed and relate source-target node pairs, where the entity represented by the source node is always hierarchically one level below the entity represented by the target node.

For instance, if a specific polygon p contains a specific edge e, the nodes representing these two entities, i.e., (p) and (e), are connected as (e)-[:IS_PART_OF]->(p). Connecting equally labeled nodes (e.g., a pair of edges) is not allowed and neither is connecting nodes carrying labels that do not follow one another hierarchically (e.g., a vertex-polygon pair). For example, say that one of the polygon p’s vertices is vertex v. Rather than encoding the information that v is part of p directly by a (v)-[:IS_PART_OF]->(p) connection, this information is retrieved hierarchically from the connectivity of v and p through edges. In particular, if v is part of p, it is also necessarily shared by two edges that are both also part of p (say e1 and e2): (v)-[:IS_PART_OF]->(e1)-[:IS_PART_OF]->(p) and (v)-[:IS_PART_OF]->(e2)-[:IS_PART_OF]->(p). An extra connection (v)-[:IS_PART_OF]->(p) would be redundant and is therefore forbidden in GVM, since these two subgraphs already imply that v is one of polygon p’s vertices.

The above rules for the construction of the GVM’s knowledge graph can be conveniently represented by a graph, called metagraph. Much like metalanguage is a language that describes another language, metagraph is a graph that describes another graph and can be viewed as a blueprint for generating actual (valid) manifestations of that graph. From the above definitions, it is obvious that the metagraph of GVM is (:Vertex)-[:IS_PART_OF]->(:Edge)-[:IS_PART_OF]->(:Polygon)-[:IS_PART_OF]-> ->(:Cell) (Fig 2A).

thumbnail
Download:
Fig 2. Graph vertex model.

A The hierarchical structure of GVM vertex→edge→polygon→cell is defined by a metagraph. Relationships between these nodes are labeled IS_PART_OF and carry a sign property denoted by s, σ, and Σ for vertex→edge, edge→polygon, and polygon→cell relationships, respectively. B A subgraph representing a particular local cell state is obtained by pattern matching. This subgraph is then transformed by a graph transformation. C Metagraph of graph-transformation graph connects Vertex, Edge, Polygon, and Cell nodes with green and red relationships, indicating creation and deletion of IS_PART_OF relationships in the GVM’s knowledge graph, respectively.

https://doi.org/10.1371/journal.pcbi.1012089.g002

Additionally, both nodes and relationships carry properties that encode additional information about nodes and relationships. While nodes carry a property, id, which represents the identification numbers i, j, k, and l of vertices, edges, polygons, and cells, respectively, the relationships are prescribed a property sign, whose value is either −1 or +1. This is a contextual property, that puts the relationship’s source node into the context of the target node. In particular, in the subgraphs of type (v:Vertex)-[r:IS_PART_OF]->(e:Edge), the value of r.sign denotes whether vertex v is a head vertex (r.sign=+1) or a tail vertex (r.sign=-1) of edge e [i.e., parameter s in Eq (2)]. In the subgraphs of types (e:Edge)-[r:IS_PART_OF]->(p:Polygon) and (p:Polygon)-[r:IS_PART_OF]->(c:Cell) the same property specifies the orientation of edge e in the context of polygon p and polygon p in context of cell c, respectively [i.e., parameters σ and Σ in Eqs (3) and (4), respectively].

Pattern matching

Transforming the graph database of GVM upon cell rearrangements requires only two steps: (i) Data retrieval, accomplished through pattern matching and (ii) Graph transformation.

In step (i), a suitable (meta)graph pattern is utilized to query the database and identify the nodes relevant to the transformation at hand. After the graph is traversed, instances of the specified graph pattern are returned. These instances are further filtered, using various conditional statements.

The goal of step (i) is to retrieve from the whole graph of GVM a small subgraph comprising solely of the nodes representing the objects (vertices, edges, polygons, and cells) that take part in the specific topological transformation being performed (Fig 2B). Given the unambiguous definition of the graph data structure by the GVM’s metagraph (Fig 2A), the routines that perform this step are easily reproducible. We implement these routines in Cypher and find that each topological transition requires ∼10 distinct short queries, similar to the query shown in Eq (18) to retrieve the relevant data [61].

Graph transformations

In step (ii), the subgraph matched during the pattern-matching step undergoes a transformation based on the rules of the specific cell-rearrangement event being performed.

Like the matched subgraph that is being transformed, the graph transformation itself is represented by a graph. This graph contains exactly the same nodes as the matched GVM subgraph, however with much fewer relationships. In particular, the relationships in the transformation graph are of two types: (i) Green and (ii) red, indicating creations and deletions of :IS_PART_OF relationships in the actual GVM subgraph, respectively (Fig 2C).

Compared to the convoluted codes that perform topological transformations in the conventional vertex model, the task of programming the routines that perform graph topological transformations in GVM is much less challenging. Indeed, our implementation of graph-transformation routines in Cypher comprises of successive calls of ∼5 distinct short queries, similar to examples shown in Eqs (19) and (20), which merely delete and create relationships.

Transformation of contextual properties.

Values of contextual properties s, σ, and Σ to be assigned to the newly created relationships are specified in the transformation graph through the relationship property of green relationships (s′, σ′, and Σ′ for vertex→edge, edge→polygon, and polygon→cell connections, respectively). Unlike the green relationships, the red relationships do not carry any additional properties (Fig 2C).

Transformed values of certain contextual properties can be arbitrarily chosen, however, in most cases they need to be figured out from the known values of other contextual properties. Conveniently, the rules for prescribing these new values can be summarized in two compact formulas. In particular, for a new vertex→edge relationship between nodes representing vertex vi and edge ej, is calculated as (6) where k denotes vertex vk, which was part of edge ej prior to the transformation (i.e., sk,j denotes contextual property of the relationship between nodes representing vertex vk and edge ej prior to the transformation).

For a new edge→polygon relationship between nodes representing edge ei and polygon pj, is calculated as (7) Here, vertex vk is a vertex shared by edges ei and el, which are both part of polygon pj. Additionally, among the two vertices of edge ei, vertex vk is the one that was not part of edge el prior to the transformation (i.e., Eq (7) contains and not sk,l). An analogous equation to Eq (7) also holds for assigning to a newly created polygon → cell relationship.

Eqs (6) and (7) are demonstrated in the subsequent section for the case of a T1 transformation and their meaning is described in the corresponding figure caption.

T1 graph transformation

To demonstrate all steps required to perform a topological transformation in a space-filling cell aggregate, we first turn to a 2D polygonal cellular tiling, where cells rearrange through T1 transitions. Specifically, during a T1 transition, a vanishingly short edge merges into a single vertex, i.e., a four-way junction, which subsequently resolves into a new edge oriented roughly in the perpendicular direction compared to the orientation of the initial edge. As any other topological transformation, GVM performs a T1 transition in two steps as follows (Fig 3).

thumbnail
Download:
Fig 3. Graph transformation for a T1 transition.

A Graphs in the left and right columns correspond to the initial and final cell configurations, respectively. Gray arrows represent relationships labeled IS_PART_OF. The graph in the middle column shows the graph transformation, which includes green and red relationships, indicating relationship creations and deletions, respectively. Additionally, the graph transformation specifies property values of the newly created relationships. B Schematic of determining the value of contextual property for a new vertex→edge relationship generated between vertex v1 and edge e3 (). After the transformation, vertex v1 assumes the same role in the context of edge e3 as the role of v2 in the context of e3 before the transformation (s2,3). This occurs because the edge e3 merely replaces v2 (blue) with v1 (red). C Schematic of determining the value of contextual property for a new edge→polygon relationship generated between edge e5 and polygon p1 (). The calculation of (red) relies on one of the vertices of e5, i.e., v2 (green), and the edge linked to both v2 and p1, i.e., e2 (also depicted in green). The assignment of is determined based on the contextual properties of: (i) edge e2 in the context of polygon p1, σ2,1, (blue) and (ii) vertex v2 in contexts of edges e2, , and e5, , (green). In short, the contextual property aligns or opposes that of σ2,1 depending on the similarity or dissimilarity of and .

https://doi.org/10.1371/journal.pcbi.1012089.g003

(i) Pattern matching performs a series of graph-database queries so as to find nodes representing elements (i.e., vertices, edges, polygons and cells) that participate in the transformation (the initial state in Fig 3A). These queries first identify the edge that undergoes a T1 and label it e5. Subsequently, they identify two vertices connected to e5, randomly labeling one v1 and the other v2. Following this, they locate two polygons that share e5, labeling them p3 and p4. Continuing, edges e1 and e2 are identified as edges connected to v1 (excluding e5) and are also part of polygon nodes p3 and p4, respectively. Similarly, e3 and e4 are determined using the same procedure, where the role of v1 from the previous step is now adapted by v2. Subsequent steps involve finding polygons p1, containing both e1 and e2, followed by identifying p2, containing e3 and e4.

Note that nodes representing cells (polyhedra) are missing in graphs describing a T1 transition (Fig 3). This is because polyhedra do not exist in the 2D representation of the vertex model, where polygons themselves are interpreted as cells.

(ii) After pattern matching, graph transformations convert the initial sub-graph to the final sub-graph, illustrated in the upper-right area of Fig 3A, employing a few transformative operations. The detailed depiction of these graph transformations in the middle panel of Fig 3A reveals several deletions and creations of new relationships. Specifically, at the vertex and edge level, v1 → e2 and v2 → e3 relationships are eliminated, while v1 → e3 and v2 → e2 relationships are established. Furthermore, at edges and polygons level, only e5 → p3 and e5 → p4 relationships are removed, while new e5 → p1 and e5 → p2 relationships are created. What relationships need to be deleted and created is visually represented in the graph-transformations graph by red and green arrows, respectively.

Finally, the newly created relationships need to be prescribed contextual properties using Eqs (6) and (7), (Fig 3B and 3C). For instance, , i.e., the sign of vertex v1 in the context of the edge e3, is assigned a value identical to s2,3, i.e., the sign of vertex v2 in the context of e3 before transformation. In turn, determining the sign of a new edge in the context of a polygon, e.g., e5 in the context of p1 (), relies on contextual properties , , and σ2,1 [Eq (7)]. For a more comprehensive understanding of the origin of Eqs (6) and (7), refer to Fig 3’s caption.

Note that all the graphs shown in Fig 3 are unique in a sense that their connectivity does not depend on the labeling of vertices, edges, and polygons. Of course, relabeling these elements or repositioning the corresponding nodes would affect the visual representation of the graphs, but the graphs themselves (i.e., their connectivities) would remain the same.

ET and TE graph transformations

Graph transformations describing ET and TE transitions as well as EV, VT, TV, and VE transitions are obtained following the exact same procedure as in the case of T1 transition (previous section). Fig 4 shows graph transformations for ET and TE transformations as well as the matched subgraphs representing initial and final cell configurations. Graph transformations for EV, VT, TV, and VE transformations are given in S1 and S2 Figs. For clarity, the graphs representing all three states (E, T, and V) are additionally specified in a more explicit (non-pictorial) form in Methods.

thumbnail
Download:
Fig 4. Graph transformations in a 3D vertex model of polyhedral packings.

A Graph transformation of an ET transition. B Graph transformation of an TE transition. In both panels, graphs in the left and the right column correspond to the initial and the final cell configuration, respectively. Gray arrows represent relationships labeled IS_PART_OF. The graphs in the middle column show graph transformations, which include green and red relationships, indicating relationship creations and deletions, respectively. Additionally, the graph transformation specifies property values of the newly created relationships. In each graph, the node indices increase from left to right in unit steps.

https://doi.org/10.1371/journal.pcbi.1012089.g004

Generalization of topological transformations

As shown in Fig 1D, the transformation patterns during ET and TE transformations are both geometrically as well as topologically quite similar to the more simple T1 transition. This raises a question whether topological transformations in 2D and 3D can be generalized and described using the same, generalized, graph transformation.

Surprisingly, as depicted in Fig 5, our Graph vertex model readily resolves this question. Indeed, when either ET or TE transformation (Fig 5 only shows the case of ET transformation) are applied on a 2D GVM of polygonal cell aggregates, only a part of the initial subgraph is matched, whereas the rest (shown in transparent in Fig 5) corresponds to elements (vertices, edges, polygons, and cells) that do not exist in the 2D model due to the reduced dimensionality. In turn, only the matched subgraph gets transformed and by relabelling and repositioning the nodes, we prove that the remaining graph transformation exactly corresponds to the graph transformation (Fig 3A), describing a T1 transition. In short, T1 graph transformation can be viewed as a subgraph of both ET and TE graph transformations.

thumbnail
Download:
Fig 5. 3D graph transformations reduce to a T1 transition when applied to a 2D vertex model.

A ET transformation when applied to a four-polygon neighborhood described by a 2D vertex model. Parts of graphs shown in transparent are not matched because the corresponding elements do not exist in 2D. After relabeling the nodes (panel B) and repositioning them (panel C), we observe that the matched graph transformation exactly corresponds to , i.e., the graph transformation that performs a T1 transition (Fig 3).

https://doi.org/10.1371/journal.pcbi.1012089.g005

This result generalizes topological transformations in 2D and 3D vertex models, suggesting that it should be possible to develop a generalized computational implementation of GVM, capable of simulating both 2D and 3D space-filling packings. Our implementation of GVM within the neoVM package confirms this hypothesis.

Order-disorder transition in active tissues

As a proof of concept, we develop a custom Python package called neoVM, which manages the GVM’s knowledge graph and its transformations in a graph database management framework Neo4j (Methods).

We use neoVM, to study an order-disorder transition of cell aggregates, driven by active tension fluctuations. In particular, we consider an aggregate of Nc = 128 cells with identical (normalized) volumes (V0 = 1 for all cells), enclosed within a simulation box with periodic boundary conditions. The vertex dynamics are described by Eq (5), assuming the potential energy given by Eq (10).

In addition to the conservative and friction forces, we also include active force dipoles acting along cell edges to induce active cell rearrangements. The total active force on vertex i is a sum of forces acting along edges (i.e., tricellular junctions) sharing that same vertex: (8) Here the indices lmn denote a tricellular junction (i.e., edge), shared by cells l, m, and n; Llmn is the edge length.

The magnitudes of active force dipoles are dynamic quantities that fluctuate with time. In particular γlmn(t) obeys Ornstein-Uhlenbeck dynamics described by (9) where τm is the relaxation time scale associated with turnover dynamics of molecular motor Myosin, γ0 is a baseline tension, whereas ξlmn(t) is Gaussian white noise with properties 〈ξlmn(t)〉 = 0 and 〈ξlmn(t)ξopr(t′)〉 = (2σ2/τm)δloδmpδnrδ(t − t′); σ2 is the long-time variance of the tension fluctuations.

We simulate the above active dynamics at different magnitudes of active noise σ, starting with a Kelvin structure–a crystalline cell arrangement made up of truncated octahedra with 14 facets (8 regular hexagons and 6 squares). The active noise distorts the geometry of the aggregates, which are no longer perfect crystals. In particular, for σ > 0 the average cell shape, quantified by the shape factor (Sl and Vl are cell surface area and volume, respectively) deviates from the Kelvin’s truncated octahedron (Fig 6A). The distorted geometry is also seen in the width of the distribution of edge lengths, which increases with an increasing σ–to a point where vanishingly short edges appear (Fig 6B). These edges undergo ET and TE topological transformations, which in turn triggers cell rearrangements–a signature of a transition from a solid-like to a fluid-like behavior. This transition manifests in disordering of the aggregate, seen in the dependence of an order parameter 1 − f14, describing the fraction of non-14-sided polyhedra (f14 = N14/Nc), on the control parameter σ (Fig 6C). From these results, we obtain an estimate for the transition point σ* ≈ 0.17. Fig 6D–6H show cell configurations at σ = 0, 0.2, 0.3, 0.4, and 0.5.

thumbnail
Download:
Fig 6. Order-disorder transition in 3D cell aggregates.

A Shape parameter q for “thermalized” aggregates versus active noise σ. B Distribution of edge lengths in thermalized aggregates for σ = 0.025, 0.2, and 0.5 (red, green, and blue curves, respectively). C Order parameter 1 − f14 and normalized number of cell-rearrangement events n = (NET + NTE)/Nmax versus active noise σ. D-H Thermalized aggreages at σ = 0, 0.2, 0.3, 0.4, and 0.5. I Order parameter 1 − f14 versus time for σ = 0.01, 0.17 and 0.4 (red, green, and blue curves, respectively) during aggregate ordering. J The initial (t = 0.001; orange curve) and final (t = 1000; magenta curve) distributions of edge lengths at σ = 0.17. K Shape parameter q versus time for σ = 0.01, 0.17 and 0.4 (red, green, and blue curves, respectively) during aggregate ordering.

https://doi.org/10.1371/journal.pcbi.1012089.g006

Next, we are interested in whether the active noise can drive the opposite effect, i.e., ordering. To study this, we start with a disordered cell packing, prepared in advance by packing spheres in the simulation box, using random sequential addition and then constructing Voronoi partitions around sphere centers. This procedure yields a sample with the initial fraction of 14-sided polyhedra f14 = 0.234. Again, we simulate the active dynamics at different magnitudes of the active noise σ. We find that high-σ values keep the aggregate disordered due to frequent cell-neighbor exchanges. In contrast, a sufficiently small level of active noise drives active tissue ordering, which is seen in decreasing 1 − f14 and narrowing of the edge-length distribution over time (Fig 6I and 6J, respectively). At moderate σ ≈ σ* values, the ordering is more efficient compared to small σ values, where the active-noise level is not sufficiently high to allow overcoming local energy barriers for cell rearrangements. Despite a higher degree of disorder, the low-σ states consist of cells whose shapes are closer to the Kelvin’s regular truncated octahedron compared to cells from the σ ≈ σ* case, where cell shapes are perturbed due to active noise (Fig 6K).

The rate of topological-transition events in the simulations of active cell aggregates reaches as high as ∼200/time unit, which in total amounts to ∼105 events per simulation (Fig 6C); note that many of these events are reversible transitions that occur multiple times while the manipulated vertices are still located very close to one another and have not yet properly resolved in the geometric sense. Despite this large number of reconnections, the aggregate does not develop any nonphysical topology, e.g., an edge connecting more than two vertices or a polygon with one missing edge, etc. This demonstrates that topological transformations in space-filling 3D packings can indeed be implemented as graph transformations defined in Fig 4. Importantly, due to GVM’s unambiguous data structure, these transformations are relatively straight-forward to implement and are therefore readily reproducible, which is clearly demonstrated by the implementation of GVM within our neoVM package [61].

Discussion

We reformulated the vertex model of cell aggregates. The new formulation, called Graph Vertex Model (GVM), is based on storing the topology of the cell network into a knowledge graph. We discovered a particular graph data model, uniquely defined by a metagraph, which allows formulating topological transformations of the cell network as simple and mathematically properly defined graph transformations (Fig 2). These transformations are themselves represented by graphs and consist of only the most elementary graph operations, i.e., relationship deletions and creations.

We designed graph transformations for all topological transitions which are required to describe cell rearrangements–including edge-to-triangle (ET) and triangle-to-edge (TE) transformations (Fig 4). Importantly, ET and TE transformations can be both applied to a 2D system, where they reduce to the well-known T1 transition. We showed that this happens because the transformation graph describing a T1 transition is in fact a subgraph of both ET and TE transformation graphs (Figs 3 and 5). Thus, when ET or TE transformations are applied onto a 2D polygonal tiling, only the T1 subgraph is matched and the corresponding transformation executed. This result generalizes topological transformations in 2D and 3D space-filling packings, suggesting that our proposed graph-data structure may be the most natural representation of the topology of space-filling packings. We used GVM to study active cell aggregates, whose cell junctions are subject to fluctuating active tensions (Fig 6). In particular, we characterized the order-disorder transition and found initially disordered aggregates undergoing ordering, which is most efficient for active noise close to the transition point.

Even though the basic GVM’s data model presented here only encodes information on the topology of the network of cell-cell interfaces, GVM already represents an important technological and conceptual step forward in computational models of tissue mechanics. This advancement lays the foundation for the development of knowledge graphs capable of structurally storing live-imaging data, such as that obtained from developing embryos, integrating data on geometry, topology, mechanics, and biochemistry. With this aim, our ongoing work uses GVM as a starting point to develop a comprehensive knowledge-graph database of the early fly development. Making this database interactive and accessible online will allow collaborative research with the aim to progressively expand our collective knowledge base about the mechanics of the embryonic development. Additionally, its graph data structure may even be readily complemented with graph-compatible methods of artificial intelligence (e.g., Graph Neural Networks).

Graph vertex model can be readily extended to describe other cell-scale events that change the local topology of the network of cell-cell interfaces. Indeed, integrating cell division and apoptosis into the model would allow detailed mechanistic studies of spheroid-like cell aggregates as models of tumors or even embryos during the early stages of development, where cells may still be packed within a three-dimensional aggregate [70, 71]. In the context of tumors, for instance, investigating how the stability of the overall tumor shape depends on smaller-scale biomechanical processes such as the cells’ effective surface tension, inhomogeneous cell proliferation, and active noise at cell-cell interfaces, would allow better understanding the mechanical basis of tumors’ transition to malignancy.

For more efficient simulations, capable of dealing with hundreds or even thousands of cells, a considerable effort will also have to be devoted to developing technologies that will improve computational efficiency of the vertex-model simulations over a graph database. While our current implementation of GVM, neoVM [61], primarily serves as a proof of concept, it falls short on the efficiency. The reason for this is mostly twofold: (i) The time integration of the dynamical system is performed in Python, which is generally slower than some low-level programming languages, and (ii) Performing operations on the graph database managed by Neo4j necessitates reading from and writing to the local hard drive during cell rearrangements, where the database is stored. That this indeed limits the performance of neoVM, is also seen in S4 Fig, which shows the simulation time depending on the system size. While the dependence itself is linear, tissue activity shifts this linear dependence towards a higher offset value following a larger number of cell-rearrangement events. This shows that cell rearrangements consume most of the simulation time of neoVM, suggesting that the efficiency of neoVM could be significanlty improved by relying on memory storage instead of accessing the local hard drive.

Materials and methods

Initial configurations

Initial configurations (.vt3d files) used in simulations (Fig 6 and S4 Fig) are available in the git repository of neoVM. Among these initial files, the Kelvin initial configuration is generated by taking a unit cell from Surface evolver [72] and replicating them in all directions until the number of cells in the aggregate reaches 128. On the other hand, the disordered initial conditions are generated by putting spheres in the simulation box using random sequential addition and, subsequently, constructing Voronoi cells around these spheres using the Voro++ package [73].

Calculation of forces

We neglect the inertial effects so that the friction force needs to counterbalance the sum of conservative and active forces, and , respectively. This implies , where the conservative force , with W being the potential energy of the system, whereas . Here, only friction with a static background is considered, η being the associated friction coefficient. The active force can describe different system-specific active mechanisms, e.g., active contractions of the cell membrane due to the activity of the underlying cell cortex or traction forces [19, 20, 36, 74]. This model yields a system of first-order dynamic equations for vertex positions given by Eq (5).

We consider a model, in which cell-cell interfaces are prescribed by effective surface tensions Γlm, which include contributions of the cell cortical tension and cell-cell adhesion [75–77]; the notation lm denotes index of a polygon shared by cells l and m. In this model, the total potential energy of the cell aggregate reads (10) where the first sum goes over all pairs of neighboring cells l and m and the second sum goes over all the cells, κV and V0 being the cell-incompressibility constant and the preferred cell volume, respectively. Note that V0 need not be equal for all cells. In heterogeneous cell aggregates such as tumors, a similar same energy functional could be used, but V0’s would need to be considered distributed.

By definition, the conservative force acting on vertex i is calculated as (11) which further requires calculating gradients of interfacial surface areas and cell volumes as described below.

Surface area of polygonal side k is calculated as a sum of surface areas of triangular surface elements, ‖alm,μ‖, defined by pairs of consecutive polygon vertices rμ and rμ+1, and the polygon’s center of mass (12)

Like in the previous section, Greek indices here do not denote the real vertex identification numbers, but their sequential indices within individual polygons. The surface area of polygon lm reads (13) and its gradient (14) (15) where δij is the Kronecker delta.

Cell volume l is calculated as a sum of volumes of tetrahedra, defined by triangular surface elements (rμ, rμ+1, clm) and with the fourth vertex at the origin (0, 0, 0), as (16)

Its gradient is calculated as (17)

Topological transformations

None of the topological transitions is allowed if the resulting cell configuration breaks any of the following topological rules [51, 57]: (i) Edge pairs may not share more than one vertex, (ii) polygon pairs may not share more than one edge and (iii) cell pairs may not share more than one polygon. In our implementation of GVM within the neoVM packing, these conditions seem to be rigorously imposed, since we never observe them being violated despite conducting numerous simulations involving a large number of rearrangement events (Fig 6). If such violation were to happen, the transformation causing it would need to be immediately reversed.

Cypher queries for pattern matching and graph transformations

The following Cypher query retrieves nodes representing the polygons that share common edge i (18)

The following Cypher query creates a new relationship IS_PART_OF between nodes (v1) and (e3) and assigns property value s23. (19)

The following Cypher query deletes relationship r between nodes (e5) and (p4). (20)

List of relationships

For clarity, we here explicitly list relationships in graphs representing states E, T, and V (Figs 4 and 5, S1 and S2 Figs). All relationships are of type IS_PART_OF.

  • [state E]: v1 → e1, v1 → e2, v1 → e5, v1 → e7, v2 → e3, v2 → e4, v2 → e6, v2 → e7, e1 → p1, e1 → p4, e1 → p8, e2 → p2, e2 → p4, e2 → p6, e3 → p1, e3 → p5, e3 → p9, e4 → p2, e4 → p5, e4 → p7, e5 → p3, e5 → p6, e5 → p8, e6 → p3, e6 → p7, e6 → p9, e7 → p1, e7 → p2, e7 → p3, p1 → c1, p1 → c2, p2 → c1, p2 → c3, p3 → c2, p3 → c3, p4 → c1, p4 → c4, p5 → c1, p5 → c5, p6 → c3, p6 → c4, p7 → c3, p7 → c5, p8 → c2, p8 → c4, p9 → c2, p9 → c5
  • [state T]: v1 → e1, v1 → e3, v1 → e7, v1 → e9, v2 → e2, v2 → e4, v2 → e7, v2 → e8, v3 → e5, v3 → e6, v3 → e8, v3 → e9, e1 → p1, e1 → p4, e1 → p8, e2 → p2, e2 → p4, e2 → p6, e3 → p1, e3 → p5, e3 → p9, e4 → p2, e4 → p5, e4 → p7, e5 → p3, e5 → p6, e5 → p8, e6 → p3, e6 → p7, e6 → p9, e7 → p4, e7 → p5, e7 → p10, e8 → p6, e8 → p7, e8 → p10, e9 → p8, e9 → p9, e9 → p10, p1 → c1, p1 → c2, p2 → c1, p2 → c3, p3 → c2, p3 → c3, p4 → c1, p4 → c4, p5 → c1, p5 → c5, p6 → c3, p6 → c4, p7 → c3, p7 → c5, p8 → c2, p8 → c4, p9 → c2, p9 → c5, p10 → c4, p10 → c5
  • [state V]: v1 → e1, v1 → e2, v1 → e3, v1 → e4, v1 → e5, v1 → e6, e1 → p1, e1 → p4, e1 → p8, e2 → p2, e2 → p4, e2 → p6, e3 → p1, e3 → p5, e3 → p9, e4 → p2, e4 → p5, e4 → p7, e5 → p3, e5 → p6, e5 → p8, e6 → p3, e6 → p7, e6 → p9, p1 → c1, p1 → c2, p2 → c1, p2 → c3, p3 → c2, p3 → c3, p4 → c1, p4 → c4, p5 → c1, p5 → c5, p6 → c3, p6 → c4, p7 → c3, p7 → c5, p8 → c2, p8 → c4, p9 → c2, p9 → c5

Implementation in Python and Neo4j

As a proof of concept, we set up the GVM’s knowledge-graph database in a graph database management framework Neo4j [62]. The core program of the vertex model is implemented in Python and communicates with Neo4j using Py2neo client library. The time integration of the dynamical system is performed in Python, whereas all topological transformations are performed in Neo4j through pattern matching and graph transformations implemented as Cypher queries [69]. Our implementation of GVM is available as an open-source Python package, called neoVM, and is available online [61].

S3 Fig shows the schematic of neoVM’s architecture. The program is initialized by reading the initial geometry and topology of the cell network from an input .vt3d file and storing them into an object t of class tissue. In particular, this object stores lists of vertex, edge, polygon, and cell objects, which encode ri, ej, pk, and cl, respectively, in a tabular form. This is followed by generating an object db of class database, which connects to an empty Neo4j database and fills it with the tissue data according to the rules of the GVM’s metagraph, using function setup_DB(). The initialization is followed by a time loop, which propagates the system forward in time by time steps Δt. At each time step, the dynamical system [Eq (5)] is integrated between t and t + Δt [function solver()] and the program checks whether any of the edges in the cell network meets contitions for topological transitions [function topological_transitions()]. In particular, edge j undergoes an ET transition if it is shorter than a threshold length lth and the rate of change of its length is negative (), i.e., the edge is contracting. If the edge happens to be part of a triangle, a TE transition is performed on the triangle if, additionally, all edges of the triangle are shorter than lth and the area of the triangle is decreasing (dAk/dt < 0).

For every edge/triangle, subject to a topological transition, the core program sends a sequence of Cypher queries to the Neo4j graph database. These queries (i) perform pattern matching to isolate a subgraph relevant for the particular transformation being performed, and (ii) perform graph transformation on that subgraph. After graph transformations, vertex positions are displaced such that the lengths of the newly created edges (edges e7, e8, and e9 in Fig 1) are on the order of lnew ≪ 1. Finally, the local structure of the arrays, encoding ri, ej, pk, and cl, (stored in object t) are updated according to the applied transformations. This is done by converting the altered part of the knowledge graph back into the array format using function update().

Supporting information

S1 Fig. Graph transformations of EV and VT transitions.

(panels A and B, respectively). Graphs in the left and the right column correspond to the initial and the final cell configuration, respectively. Gray arrows represent relationships labeled IS_PART_OF. The graphs in the middle column show graph transformations, which include green and red relationships, indicating relationship creations and deletions, respectively. Additionally, the graph transformation specifies property values of the newly created relationships. In each graph, the node indices increase from left to right in unit steps.

https://doi.org/10.1371/journal.pcbi.1012089.s001

(TIF)

S2 Fig. Graph transformations of TV and VE transitions.

(panels A and B, respectively). Graphs in the left and the right column correspond to the initial and the final cell configuration, respectively. Gray arrows represent relationships labeled IS_PART_OF. The graphs in the middle column show graph transformations, which include green and red relationships, indicating relationship creations and deletions, respectively. Additionally, the graph transformation specifies property values of the newly created relationships. In each graph, the node indices increase from left to right in unit steps.

https://doi.org/10.1371/journal.pcbi.1012089.s002

(TIF)

S3 Fig. The algorithm underlying neoVM.

A 3D cell aggregate is set up from an input .vt3d file by function setup_from_vt3d() and then converted into a graph database, set up in Neo4j by function setupDB(). These initialization steps then followed by a time loop, which iterates between solver() and topological_transitions() functions. The function solver() calculates geometric properties of cells and the associated gradients (i.e., conservative forces) and propagates the system forward in time. The function topological_transitions() loops over all cell edges to find those that meet criteria for topological transitions. For edges that meet these criteria, topological transitions are performed through pattern matching and graph transformations (functions pattern_matching() and graph_transformation(), respectively), applied directly to the graph database. Finally, the tissue is updated accordingly by function update() and positions of vertices involved in the topological transformation are corrected by function correct_vertex_positions().

https://doi.org/10.1371/journal.pcbi.1012089.s003

(TIF)

S4 Fig. Performance of neoVM.

Runtime duration per unit simulation time as a function of cell count. Here, two separate data sets represent two strengths of active fluctuations, i.e., σ = 0.2 and σ = 0.4, shown in red and green dots, respectively. In both cases, simulation runtime shows a linear dependence (in a blue dashed lines) on the size of cell aggregates.

https://doi.org/10.1371/journal.pcbi.1012089.s004

(TIF)

Acknowledgments

We thank Urban Železnik for generating the disordered initial configurations for benchmarking. We thank Tomer Stern, Domen Vaupotič, Staš Adam, and all members of the Theoretical Biophysics Group at Jožef Stefan Institute for fruitful discussions. The main idea for this project surfaced from a very fruitful collaboration with Wisdom Labs. Thanks to all Wisdom Labs folks, especially Andraž Tuš, Kristjan Pečanac, Luka Stopar, Marko Zadravec, Jasna Pečanac, and Jan Lampič.

References

  1. 1. Lecuit T, Lenne PF. Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis. Nature Reviews Molecular Cell Biology. 2007; 8(8):633–644. pmid:17643125
  2. 2. Angelini TE, Hannezo E, Trepat X, Marquez M, Fredberg JJ, Weitz DA. Glass-like dynamics of collective cell migration. Proceedings of the National Academy of Sciences. 2011; 108(12):4714–4719. pmid:21321233
  3. 3. Park JA, Kim JH, Bi D, Mitchel JA, Qazvini NT, Tantisira K, et al. Unjamming and cell shape in the asthmatic airway epithelium. Nature Materials. 2015; 14(10):1040–1048. pmid:26237129
  4. 4. Mongera A, Rowghanian P, Gustafson HJ, Shelton E, Kealhofer DA, Carn EK, et al. A fluid-to-solid jamming transition underlies vertebrate body axis elongation. Nature. 2018; 561(7723):401–405. pmid:30185907
  5. 5. Etournay R, Popović M, Merkel M, Nandi A, Blasse C, Aigouy B, et al. Interplay of cell dynamics and epithelial tension during morphogenesis of the Drosophila pupal wing. eLife. 2015; 4:e07090. pmid:26102528
  6. 6. Discher DE, Janmey P, Wang YL. Tissue cells feel and respond to the stiffness of their substrate. Science. 2005; 310(5751):1139–1143. pmid:16293750
  7. 7. Fiore VF, Krajnc M, Quiroz FG, Levorse J, Pasolli HA, Shvartsman SY, et al. Mechanics of a multilayer epithelium instruct tumour architecture and function. Nature. 2020; 585(7825):433–439. pmid:32879493
  8. 8. Blankenship JT, Backovic ST, Sanny JSP, Weitz O, Zallen JA. Multicellular rosette formation links planar cell polarity to tissue morphogenesis. Developmental Cell. 2006; 11(4):459–470. pmid:17011486
  9. 9. Martin AC, Kaschube M, Wieschaus EF. Pulsed contractions of an actin–myosin network drive apical constriction. Nature. 2009; 457(7228):495–499. pmid:19029882
  10. 10. Rauzi M, Krzic U, Saunders TE, Krajnc M, Ziherl P, Hufnagel L, et al. Embryo-scale tissue mechanics during Drosophila gastrulation movements. Nature Communications. 2015; 6(1):8677. pmid:26497898
  11. 11. Streichan SJ, Lefebvre MF, Noll N, Wieschaus EF, Shraiman BI. Global morphogenetic flow is accurately predicted by the spatial distribution of myosin motors. eLife. 2018; 7:e27454. pmid:29424685
  12. 12. Stern T, Shvartsman SY, Wieschaus EF. Template-based mapping of dynamic motifs in tissue morphogenesis. PLOS Computational Biology. 2020; 16(8):e1008049. pmid:32822341
  13. 13. Stern T, Shvartsman SY, Wieschaus EF. Deconstructing gastrulation at single-cell resolution. Current Biology. 2022; 32(8):1861–1868. pmid:35290798
  14. 14. Bertet C, Sulak L, Lecuit T. Myosin-dependent junction remodelling controls planar cell intercalation and axis elongation. Nature. 2004; 429(6992):667–671. pmid:15190355
  15. 15. Rauzi M, Verant P, Lecuit T, Lenne PF. Nature and anisotropy of cortical forces orienting Drosophila tissue morphogenesis. Nature Cell Biology. 2008; 10(12):1401–1410. pmid:18978783
  16. 16. Rauzi M, Lenne PF, Lecuit T. Planar polarized actomyosin contractile flows control epithelial junction remodelling. Nature. 2010; 468(7327):1110–1114. pmid:21068726
  17. 17. Lecuit T, Lenne PF, Munro E. Force generation, transmission, and integration during cell and tissue morphogenesis. Annual Review of Cell and Developmental Biology. 2011; 27:157–184. pmid:21740231
  18. 18. Lecuit T, Yap AS. E-cadherin junctions as active mechanical integrators in tissue dynamics. Nature Cell Biology. 2015; 17(5):533–539. pmid:25925582
  19. 19. Curran S, Strandkvist C, Bathmann J, de Gennes M, Kabla A, Salbreux G, et al. Myosin II controls junction fluctuations to guide epithelial tissue ordering. Developmental Cell. 2017; 43(4):480–492. pmid:29107560
  20. 20. Krajnc M. Solid–fluid transition and cell sorting in epithelia with junctional tension fluctuations. Soft Matter. 2020; 16(13):3209–3215. pmid:32159536
  21. 21. Krajnc M, Dasgupta S, Ziherl P, Prost J. Fluidization of epithelial sheets by active cell rearrangements. Physical Review E. 2018; 98(2):022409. pmid:30253464
  22. 22. Krajnc M, Stern T, Zankoc C. Active instability and nonlinear dynamics of cell-cell junctions. Physical Review Letters. 2021; 127(19):198103. pmid:34797151
  23. 23. Staddon MF, Cavanaugh KE, Munro EM, Gardel ML, Banerjee S. Mechanosensitive junction remodeling promotes robust epithelial morphogenesis. Biophysical Journal. 2019; 117(9):1739–1750. pmid:31635790
  24. 24. Cavanaugh KE, Staddon MF, Chmiel TA, Harmon R, Budnar S, Yap AS, et al. Force-dependent intercellular adhesion strengthening underlies asymmetric adherens junction contraction. Current Biology. 2022; 32(9):1986–2000. pmid:35381185
  25. 25. Sknepnek R, Djafer-Cherif I, Chuai M, Weijer C, Henkes S. Generating active T1 transitions through mechanochemical feedback. eLife. 2023; 12:e79862. pmid:37039463
  26. 26. Keller R. Cell migration during gastrulation. Current Opinion in Cell Biology. 2005; 17(5):533–541. pmid:16099638
  27. 27. Ajeti V, Tabatabai AP, Fleszar AJ, Staddon MF, Seara DS, Suarez C, et al. Wound healing coordinates actin architectures to regulate mechanical work. Nature Physics. 2019; 15(7):696–705. pmid:31897085
  28. 28. Janiszewska M, Primi MC, Izard T. Cell adhesion in cancer: Beyond the migration of single cells. Journal of Biological Chemistry. 2020; 295(8):2495–2505. pmid:31937589
  29. 29. Huebsch N, Mooney DJ. Inspiration and application in the evolution of biomaterials. Nature. 2009; 462(7272):426–432. pmid:19940912
  30. 30. Alert R, Trepat X. Physical models of collective cell migration. Annual Review of Condensed Matter Physics. 2020; 11:77–101.
  31. 31. Graner F, Glazier JA. Simulation of biological cell sorting using a two-dimensional extended Potts model. Physical Review Letters. 1992; 69(13):2013–2016. pmid:10046374
  32. 32. Glazier JA, Graner F. Simulation of the differential adhesion driven rearrangement of biological cells. Physical Review E. 1993; 47(3):2128–2154. pmid:9960234
  33. 33. Swat MH, Thomas GL, Belmonte JM, Shirinifard A, Hmeljak D, Glazier JA. Multi-scale modeling of tissues using CompuCell3D. Methods in Cell Biology. 2012; 110:325–366. pmid:22482955
  34. 34. Bi D, Lopez JH, Schwarz JM, Manning ML. A density-independent glass transition in biological tissues. Nature Physics. 2015; 11:1074–1079.
  35. 35. Misra M, Audoly B, Kevrekidis IG, Shvartsman SY. Shape transformations of epithelial shells. Biophysical Journal. 2016; 110(7):1670–1678. pmid:27074691
  36. 36. Bi D, Yang X, Marchetti MC, Manning ML. Motility-driven glass and jamming transitions in biological tissues. Physical Review X. 2016; 6(2):021011. pmid:28966874
  37. 37. Mueller R, Yeomans JM, Doostmohammadi A. Emergence of active nematic behavior in monolayers of isotropic cells. Physical Review Letters. 2019; 122(4):048004. pmid:30768306
  38. 38. Kim S, Pochitaloff M, Stooke-Vaughan GA, Campàs O. Embryonic tissues as active foams. Nature Physics. 2021; 17(7):859–866. pmid:34367313
  39. 39. Barton DL, Henkes S, Weijer CJ, Sknepnek R. Active vertex model for cell-resolution description of epithelial tissue mechanics. PLOS Computational Biology. 2017; 13(6):e1005569. pmid:28665934
  40. 40. Osborne JM, Fletcher AG, Pitt-Francis JM, Maini PK, Gavaghan DJ. Comparing individual-based approaches to modelling the self-organization of multicellular tissues. PLOS Computational Biology. 2017; 13(2):e1005387. pmid:28192427
  41. 41. Beatrici C, Kirch C, Henkes S, Graner F, Brunnet L. Comparing individual-based models of collective cell motion in a benchmark flow geometry. Soft Matter. 2023; 19(29):5583–5601. pmid:37439121
  42. 42. Buttenschön A, Edelstein-Keshet L. Bridging from single to collective cell migration: A review of models and links to experiments. PLOS Computational Biology. 2020; 16(12):e1008411. pmid:33301528
  43. 43. Honda H, Tanemura M, Nagai T. A three-dimensional vertex dynamics cell model of space-filling polyhedra simulating cell behavior in a cell aggregate. Journal of Theoretical Biology. 2004; 226(4):439–453. pmid:14759650
  44. 44. Farhadifar R, Röper JC, Aigouy B, Eaton S, Jülicher F. The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing. Current Biology. 2007; 17(24):2095–2104. pmid:18082406
  45. 45. Fletcher AG, Osterfield M, Baker RE, Shvartsman SY. Vertex models of epithelial morphogenesis. Biophysical Journal. 2014; 106(11):2291–2304. pmid:24896108
  46. 46. Alt S, Ganguly P, Salbreux G. Vertex models: from cell mechanics to tissue morphogenesis. Philosophical Transactions of the Royal Society B: Biological Sciences. 2017; 372(1720):20150520. pmid:28348254
  47. 47. Theis S, Suzanne M, Gay G. Tyssue: an epithelium simulation library. Journal of Open Source Software. 2021; 6(62):2973.
  48. 48. Mirams GR, Arthurs CJ, Bernabeu MO, Bordas R, Cooper J, Corrias A, et al. Chaste: an open source C++ library for computational physiology and biology. PLOS Computational Biology. 2013; 9(3):e1002970. pmid:23516352
  49. 49. Sussman DM. cellGPU: Massively parallel simulations of dynamic vertex models. Computer Physics Communications. 2017; 219:400–406.
  50. 50. Zhang T, Schwarz JM. TVM: https://github.com/ZhangTao-SJTU/tvm; 2023.
    • 51. Okuda S, Inoue Y, Eiraku M, Sasai Y, Adachi T. Reversible network reconnection model for simulating large deformation in dynamic tissue morphogenesis. Biomechanics and Modeling in Mechanobiology. 2013; 12:627–644. pmid:22941051
    • 52. Okuda S, Miura T, Inoue Y, Adachi T, Eiraku M. Combining Turing and 3D vertex models reproduces autonomous multicellular morphogenesis with undulation, tubulation, and branching. Scientific Reports. 2018; 8(1):2386. pmid:29402913
    • 53. Sui L, Alt S, Weigert M, Dye N, Eaton S, Jug F, et al. Differential lateral and basal tension drive folding of Drosophila wing discs through two distinct mechanisms. Nature Communications. 2018; 9(1):4620. pmid:30397306
    • 54. Rozman J, Krajnc M, Ziherl P. Collective cell mechanics of epithelial shells with organoid-like morphologies. Nature Communications. 2020; 11(1):3805. pmid:32732886
    • 55. Rozman J, Krajnc M, Ziherl P. Morphologies of compressed active epithelial monolayers. The European Physical Journal E. 2021; 44(7):99. pmid:34287727
    • 56. Okuda S, Sato K. Polarized interfacial tension induces collective migration of cells, as a cluster, in a 3D tissue. Biophysical Journal. 2022; 121(10):1856–1867. pmid:35525240
    • 57. Zhang T, Schwarz JM. Topologically-protected interior for three-dimensional confluent cellular collectives. Physical Review Research. 2022; 4(4):043148.
    • 58. Lawson-Keister E, Zhang T, Nazari F, Fagotto F, Manning ML. Differences in boundary behavior in the 3D vertex and Voronoi models. PLOS Computational Biology. 2024; 20(1):e1011724. pmid:38181065
    • 59. Honda H, Nagai T. Cell models lead to understanding of multi-cellular morphogenesis consisting of successive self-construction of cells. The Journal of Biochemistry. 2015; 157(3):129–136. pmid:25552548
    • 60. Honda H, Eguchi G. How much does the cell boundary contract in a monolayered cell sheet? Journal of Theoretical Biology. 1980; 84(3):575–588.
    • 61. Sarkar T, Krajnc M. neoVM: https://gitlab.com/ijskrajncgroup1/neovm; 2023.
      • 62. Available from: https://neo4j.com/product.
        • 63. Weaire D, Fortes MA. Stress and strain in liquid and solid foams. Advances in Physics. 1994; 43(6):685–738.
        • 64. Reinelt DA, Kraynik AM. Simple shearing flow of a dry Kelvin soap foam. Journal of Fluid Mechanics. 1996; 311:327–343.
        • 65. Schwarz HW. Rearrangements in polyhedric foam. Recueil des Travaux Chimiques des Pays‐Bas. 1965;84(6):771–781.
        • 66. Gómez-Gálvez P, Vicente-Munuera P, Tagua A, Forja C, Castro AM, Letrán M, et al. Scutoids are a geometrical solution to three-dimensional packing of epithelia. Nature Communications. 2018; 9(1):2960. pmid:30054479
        • 67. Rupprecht JF, Ong KH, Yin J, Huang A, Dinh HHQ, Singh AP, et al. Geometric constraints alter cell arrangements within curved epithelial tissues. Molecular Biology of the Cell. 2017; 28(25):3582–3594. pmid:28978739
        • 68. Paulheim H. Knowledge graph refinement: A survey of approaches and evaluation methods. Semantic Web. 2017; 8(3):489–508.
        • 69. Available from: https://neo4j.com/docs/cypher-manual/current/introduction/cypher_overview/.
          • 70. Liu S, Lemaire P, Munro E, Mani M. A mechanical atlas for Ascidian gastrulation. bioRxiv. 2023.
          • 71. Ichbiah S, Delbary F, McDougall A, Dumollard R, Turlier H. Embryo mechanics cartography: inference of 3D force atlases from fluorescence microscopy. Nature Methods. 2023; 20(12):1989–1999. pmid:38057527
          • 72. Brakke KA. The surface evolver. Experimental Mathematics. 1992; 1(2):141–165.
          • 73. Rycroft CH. VORO++: A three-dimensional Voronoi cell library in C++. Chaos. 2009; 19:041111. pmid:20059195
          • 74. Trepat X, Wasserman MR, Angelini TE, Millet E, Weitz DA, Butler JP, et al. Physical forces during collective cell migration. Nature Physics. 2009; 5(6):426–430.
          • 75. Derganc J, Svetina S, Žekš B. Equilibrium mechanics of monolayered epithelium. Journal of Theoretical Biology. 2009; 260(3):333–339. pmid:19576229
          • 76. Manning ML, Foty RA, Steinberg MS, Schoetz EM. Coaction of intercellular adhesion and cortical tension specifies tissue surface tension. Proceedings of the National Academy of Sciences. 2010; 107(28):12517–12522. pmid:20616053
          • 77. Hannezo E, Prost J, Joanny JF. Theory of epithelial sheet morphology in three dimensions. Proceedings of the National Academy of Sciences. 2014; 111(1):27–32. pmid:24367079