Mechanisms of Regulating Cell Topology in Proliferating Epithelia: Impact of Division Plane, Mechanical Forces, and Cell Memory

Regulation of cell growth and cell division has a fundamental role in tissue formation, organ development, and cancer progression. Remarkable similarities in the topological distributions were found in a variety of proliferating epithelia in both animals and plants. At the same time, there are species with significantly varied frequency of hexagonal cells. Moreover, local topology has been shown to be disturbed on the boundary between proliferating and quiescent cells, where cells have fewer sides than natural proliferating epithelia. The mechanisms of regulating these topological changes remain poorly understood. In this study, we use a mechanical model to examine the effects of orientation of division plane, differential proliferation, and mechanical forces on animal epithelial cells. We find that regardless of orientation of division plane, our model can reproduce the commonly observed topological distributions of cells in natural proliferating animal epithelia with the consideration of cell rearrangements. In addition, with different schemes of division plane, we are able to generate different frequency of hexagonal cells, which is consistent with experimental observations. In proliferating cells interfacing quiescent cells, our results show that differential proliferation alone is insufficient to reproduce the local changes in cell topology. Rather, increased tension on the boundary, in conjunction with differential proliferation, can reproduce the observed topological changes. We conclude that both division plane orientation and mechanical forces play important roles in cell topology in animal proliferating epithelia. Moreover, cell memory is also essential for generating specific topological distributions.


Geometric Model of Cell
The underlying physics of our two-dimensional cellular model is that of a well-studied topic of surface energy minimization, also called bubble formation (1; 2; 3). In this model, each cell minimizes its surface energy, as prominently observed in the developing retina of drosophila, where differential expression of N-cadherin leads to the formation of an overall shape that minimizes their surface contact with surrounding cells (4). Formally, a biological cell is represented by the combination of three types of geometric elements ( Fig S1). First, a geometric cell c i is a spatial region representing the volume of cell i (Fig S1a). It is a disk when in isolation, but can be a disk segment when the cell is contacting other cell(s). It can take the shape of the union of a polygon and one or more disk segment(s), when there are multiple contacting cells in the surroundings, and at the same time there is one or more free boundaries (Fig S1b-c). When completely surrounded by other cells, it takes the form of a polygon (Fig S1d). Cells can have different sizes.
Second, edges represent the boundaries of a cell. There are two types of edges: outer edges e i for cell i and inner edges e i,j between cell i and cell j (Fig S1b). An outer edge e i is an arc or a circle, and represents the boundary between cell c i and the outside medium (denoted as c 0 ). e i is a full circle when the cell exists in isolation, but becomes one or more arcs if the medium does not fully surround the cell. An inner edge e i,j occurs when a cell c i is in contact with another cell c j . Their shared boundary is a face with constant surface curvature (1), but is modeled as a straight line segment here. An inner edge appears twice, once for each of the neighboring cell, with the order of the two indices reversed. This reflects the fact that each cell has a separate wall.
Third, vertex is the junction point of three edges. When two cells c i and c j make contact (Fig S1b), their outer edges (arcs) e i and e j intersect at two vertices, v i,0,j (indices in clockwise direction) and v j,0,i , which are also the two end-points of the inner edge e i,j (Fig S1b). When three cells c i , c j and c k intersect (Fig S1c), they form a single vertex v i,j,k . In our two-dimensional model, we assume no more than three cells can intersect as seen in soap Figure S 1: Two-dimensional cell model. a) An isolated cell is modeled as a disk. b) A cell is modeled as a disk segment when contacting other cell(s). An outer edge e i is an arc or a circle, representing the boundary between cell c i and the outside medium (denoted as c 0 ). An inner edge e i,j occurs when a cell c i is in contact with another cell c j . Their shared boundary is modeled as a straight line segment. When two cells c i and c j make contact, their outer edges (arcs) e i and e j intersect at two vertices v i,0,j and v j,0,i , which are also the two end-points of the inner edge e i,j . c) When three cells c i , c j and c k intersect, they form a vertex v i,j,k . d) A cell completely surrounded by other cells is represented as a polygon.
bubbles (1; 2; 3). That is, no more than three edges can meet at a vertex.
For a tissue consisting of n cells, we denote the set of cell centers as Z = {z 1 , · · · , z n }, where z i ∈ R 2 is the coordinates of the center of cell i. The set of edges is denoted as E ≡ {e i } ∪ {e i,j }, and the set of vertices is denoted as V = {v i,j,k }. The overall state S of a tissue with n cells is defined as: S = (Z, E, V ). It fully determines the geometric pattern of the tissue formed by these n cells.

Physical Model of Cell and Cell Growth
Stationary Model We model cells with the assumption that they are in a stationary state, in which changes are slow, and all forces in the system at every moment are balanced out by each other.
We use discrete time steps to model incremental changes of cell volume, which is dictated by the underlying biology, e.g. cell birth, cell death, cell growth and shrinkage, and changes of cell wall properties. In our model, we assume the energy of the cells reaches a minimum at the state S: Forces exerting on the system of cells at state S have a zero net sum. We model the mechanical forces using only the vertices. The geometry of the whole system, including the edges and cells, and the forces in each cell, all will follow once the vertex set V is specified. We have: Mechanical Forces A number of physical forces exist inside a cell. Cytoskeletal microfilament (9; 10; 11), intermediate filaments (12), and cell membrane all exert compression forces on a cell. In addition, there exists adhesion or alternatively repulsion force between cells. These forces are summed up and modeled as a tension force that exerts along the direction e i,j of inner edge (interior cell boundary), or along the tangent direction of outer edge e i (free cell boundary).
There are also expansion forces in a cell. These include those from microtubles (13; 10; 14; 11) and the extracellular matrix (ECM) (15). We model these forces as a pressure force that acts along the direction normal to the edge. Pressure force only exists for an inner edge and not for arc/outer edge. Tension force exists at vertices due to both inner and outer edges.
In our model, cells minimize their energy and take upon the appearance and behavior of soap bubbles (4). According to the soap bubble model, cell walls take the shape of constant mean curvature (CMC) surfaces under fixed volume and pressure conditions (1). The most notable examples are spherical shape and biconcave erythrocyte shape (16). We therefore model cells as intersecting circular disks. Physical forces are modeled to act at vertices (Fig S2). For inner edges, although the physical cell wall will adopt a curved surface under the bubble model, the curvature is small and we simplify it as a straight line segment. The physical forces originally tangent to the curved surface are now decomposed into two components for the straight line segment. They are the tension force T (e i,j ) and the pressure force P (e i,j ). Tension force acts in the direction of shortening the edge, and pressure force follows the direction of the difference of pressure in two cells sharing the edge.
For outer edges, we take the curved surface into account and model it as an arc. In this case, the physical force tangent to the curved surface is modeled directly as the tension force T (e i ) for this cell, in the direction of shortening the arc. The pressure inside the boundary cell containing this outer edge determines the curvature of the outer edge, and therefore determines the tangential direction of the tension force.
Formally, the forces applied to a vertex v in V can be decomposed as which sums over all edges e's ending in the vertex v. Here T(e) and P(e) are the forces acting on edge e through cell wall tension and intracellular pressure, respectively. For the edge e i,j between cells i and j, the tension force is always tangential to the edge e i,j : T where η is the tension coefficient, which may depend on the cell types of both cells, and e i,j is the edge vector. We assume e i,j is in the direction of shortening e i,j , otherwise, we add a coefficient "−1" in front of this formula. Pressure induced force is in the direction normal to the inner edge: where P i and P j are the pressures in two cells, |e i,j | is the length of the inner edge, and n (i,j) is the unit vector normal to the edge in the direction from the cell with higher pressure to the cell with lower pressure. Although pressure force exerts on the whole inner edge, we decompose it equivalently to the two end vertices, each distributed with 1 2 of the total pressure force P (e i,j ).
For an outer edge e i of cell i, the tension force acting on a vertex v is always tangential to the arc e i , in the direction t v of shortening e i : Here 0 denotes the outside medium, η(i, 0) is the tension coefficient of the outer edge of cell i. The internal pressure is compensated by pressure due to the curvature of the arc of the cell boundary. Therefore, pressure induced force is zero. The value of the pressure inside an outer cell i itself is determined by the curvature of the edge: Here r i is the radius of the arc.
In general, forces at the two vertices of an edge can be in different direction, which can result in displacement of the edge. The movement of an edge is the result of the volume change of the cell. The new position of the edge forms an irregular quadrilateral with the edge before the movement. We distribute this volume change to the two vertices, each with a triangle. This volume change happens in each time integration step (Fig S3). there is no tension on an inner edge, and it can be regarded as an imaginary cell wall; b) when η(i, j) = 0.5η(i, 0), there is a strong adhesion force between the two cells; c) when η(i, j) = η(i, 0) = η(j, 0), the two cells behave as if physically they have the same wall; d) when η(i, j) = 1.5η(i, 0), there is a weak adhesion force between the two cells; e) when η(i, j) ≥ 2η(i, 0), The two cells have no adhesion and behave like soccer balls. Adding an inner wall would be more costly, as it is equivalent to adding two outer walls. In this case, the overall energy of the two cells is not reduced.
Cell Shape In our model, the cell walls take the shape of constant mean curvature surfaces under fixed volume and pressure conditions. Physically, each cell has its own wall, and the surface tension η(i, j) at an inner edge depends on the properties of both cell walls. The final shape of a cell depends on the ratio of tension coefficient η(i, j) for inner edge and η(i, 0) for the outer edge. When η(i, j) = 0, there is no tension on an inner edge, and it can be regarded as an imaginary cell wall. When η(i, j) = 0.5η(i, 0), there is a strong adhesion force between the two cells. When η(i, j) = η(i, 0) = η(j, 0), the two cells behave as if physically they have the same wall. When η(i, j) ≥ 2η(i, 0), adding an inner wall would be more costly, as it is equivalent to adding two outer walls. In this case, the overall energy of the two cells is not reduced. The two cells therefore have no adhesion and behave like soccer balls.
Calculating Forces For cell i that experiences cell volume change ∆V i at time step t, the net force at each of its vertices can be calculated based on the assumption of stationary state: |F v(e) × e|. (2) Here the coefficient of 1/2 represents the change in volume of one of the two triangles formed by the division of the irregular quadrilateral (Fig S3). The summation is over all edges e of cell i and over both end vertices {v(e)} of each edge, "×" represents vector cross product, and σ is the constant for the integration step. Eqn (2) is solved to obtain the forces at each vertex.
Updating vertex position. At time t + 1, the location of a vertex v i after the volume change is updated to: where v i (t) and v i (t + 1) are locations of vertex i before and after the time step, respectively. Time t is an integer representing the number of time steps since the initial time, σ is a constant that controls the convergence rate towards stationary state, and F v i is the net force exerting at vertex location v i .

Algorithm for Calculating Stationary State
We assume all cells exist in stationary state at the end of each time step of simulation. Cells can grow or shrink during a time step. The amount of volume changes are assigned from models of underlying biological process. The altered cell volume leads to movement of cell boundaries. In addition, cell wall properties such as the surface tension coefficients may also be different at different time steps. All these changes are introduced in increments of small fractions. For each increment, we solve Eqn (2) to obtain the updated forces. We then move the vertices using Eqn (3) to their new locations. After the final increment of volume or cell property change was introduced, we continue iterations with constant volume and constant cell properties. Vertices are further moved until the system relaxes and reaches stationary state, and a balance of the forces is established (Eqn 2). This is the same as applying a gradient search method to find a local minimum of system energy of the cells (17). We then take the geometric patterns of the cells at this state as that of time step t + 1.
The procedure for computing the stationary state of the cell pattern after one time step is illustrated in Algorithm 1. Here F v i are the forces acting on vertex i; V (t) = (v 1 (t), · · · , v m (t)) is the vector of coordinates of all of the m vertices at time t; ∆V(t) = (∆V 1 , · · · , ∆V m ) is the vector of desired volume changes associated with the vertices for all cells at time t, ∆η is the vector representing desired changes in the cell properties (e.g., cell tension coefficients, cell color) for all cells at time step t. The output is the new coordinates of the vertices V (t + 1) at time step t + 1.
The overall simulation of cellular pattern formation is carried out by repeatedly applying this algorithm to model different biological phenomena, with pre-defined time-dependent volume changes and cell properties changes assigned as input. Stochasticity and other physical factors can be incorporated in schemes that assign these changes.

Topological Changes of Cellular Pattern
An important ingredient in modeling dynamic changes of cells is an accurate account of all topological changes. We discuss these changes below.
Cell Birth, Cell Division, and Cell Death Topological changes occur during cell birth, cell division, and cell death. In our model, a new cell is generated at cell birth. We model this by inserting a new disk. A new cell is also formed if an existing cell divides. For cell division, we add an edge inside the existing dividing cell.
When a cell dies gradually, its volume decreases, and it eventually disappears, making new empty room in the space. A cell can also disappear suddenly. In both cases of cell death, we carry out two primitive operations. First, the outer edge of the dying cell is removed at the moment when cell dies; Second, the inner edges of all the cells contacting the dying cell are replaced with outer edges.
Cell Contact Changes In addition to cell birth and cell death, there are three additional types of topological changes when cells grow or shrink and their boundaries move, resulting in cell rearrangement. We use three primitives to model these topological changes, which occur when the same space would be occupied by more than one cell: Edge Insertion When two cells grow, they may come in contact with each other. When this happens, we add an inner edge to represent the newly formed intersection plane.
Void Removal When three cells are grown together, new inner edges are introduced between two contacting cells. At the moment when three cells meet at a common vertex, we insertion. When two isolated cells come into contact, we add an edge to represent the intersection plane of the boundary between these two cells. (middle) Void removal. When three cells come into contact, the curve triangular empty space is replaced by a vertex where the inner edges meet. (bottom) Edge flip. When two previously disconnected cells expand to meet while pushing away two previously connected cells, we replace one inner edge with a new inner edge.
need to replace the curved triangular empty space (void) with a new vertex where the three inner edges meet.
Edge Flip When two originally disconnected cells expand and come into contact, they may squeeze away two previously contacting cells. In this case, we remove the inner edge between two cells originally in contact, and add a new inner edge between the two cells that now come into contact.
Together with the three topological changes of inserting a cell due to cell birth, inserting an inner edge due to cell division, and deleting a cell due to cell death, we exhaust all possible topological changes of cellular patterns modeled in two dimensional space. In our model, these topological changes can occur at any discrete time step during the simulation.