MDCK Cystogenesis Driven by Cell Stabilization within Computational Analogues

The study of epithelial morphogenesis is fundamental to increasing our understanding of organ function and disease. Great progress has been made through study of culture systems such as Madin-Darby canine kidney (MDCK) cells, but many aspects of even simple morphogenesis remain unclear. For example, are specific cell actions tightly coupled to the characteristics of the cell's environment or are they more often cell state dependent? How does the single lumen, single cell layer cyst consistently emerge from a variety of cell actions? To improve insight, we instantiated in silico analogues that used hypothesized cell behavior mechanisms to mimic MDCK cystogenesis. We tested them through in vitro experimentation and quantitative validation. We observed novel growth patterns, including a cell behavior shift that began around day five of growth. We created agent-oriented analogues that used the cellular Potts model along with an Iterative Refinement protocol. Following several refinements, we achieved a degree of validation for two separate mechanisms. Both survived falsification and achieved prespecified measures of similarity to cell culture properties. In silico components and mechanisms mapped to in vitro counterparts. In silico, the axis of cell division significantly affects lumen number without changing cell number or cyst size. Reducing the amount of in silico luminal cell death had limited effect on cystogenesis. Simulations provide an observable theory for cystogenesis based on hypothesized, cell-level operating principles.


Text S1. Supporting Information TIGHT JUNCTION reorganization
TIGHT JUNCTIONS (TJs) were implemented for three reasons: 1) To prevent CELLS from coming into contact with multiple LUMENS. 2) To allow two neighboring LUMENS to merge without violating (1). 3) To allow LUMENS to expand without violating (1). TJs are an accounting mechanism stored as pairs of points. A TJ is defined as pairs of points in adjacent CELLS in contact with the same LUMEN. TJs are counted around a point in two ways ( Figure S12A). 1) If the location is within a LUMEN then the surrounding points are surveyed and any two adjacent points in different CELLS are considered to be TJs and stored. 2) If the location is within a CELL then the surrounding points are surveyed and if any two adjacent points are in a different CELL and a LUMEN, then the point within the other CELL and the current point are considered a TJ and stored.
During index change, any index change that would change the number of TJs surrounding that location will be rejected ( Figure S12B). Also, any index change from a CELL to a LUMEN where there are two TJs before and after the index change will be rejected. The latter rule is included because in rare cases when a DYING CELL is shrinking it is possible for a neighboring CELL that contacts a LUMEN to come into contact with a second LUMEN without the number of TJs changing. During the execution of individual CELL logic, LUMEN merging and expansion through TJ reorganization may occur. All points within TJs are surveyed and the first possible merge or accepted expansion will be executed. Only a single merge or expansion can be executed during a simulation cycle because once one does there are new TJs and the TJ list is incorrect. Rather than update the list and redo it, the process ends and is begun again during the next simulation cycle.
LUMEN merging will always happen if it is allowed. For this to occur a point within a TJ must be neighboring another point within a TJ that is in contact with a different LUMEN (FIGURE S12C). If that occurs two of the points, each contacting a different LUMEN, will each change their index to that of one of the LUMENS, and then the two LUMENS, now in contact with each other, will merge to form a single LUMEN. Essentially, when four CELLS contacting two different LUMENS are in close proximity to each other then the two LUMENS will merge. LUMEN expansion through TJ reorganization only occurs if it is energetically favorable as in the index change step detailed in the text, except that the TJ changing penalties are not assessed. For a point within a TJ to change its index into LUMEN, it first checks that no neighboring points are TJs contacting other LUMENS, that no neighboring points are contacting other LUMENS, and that no neighboring points are MATRIX or UNPOLARIZED CELLS ( Figure S12D). Then the point will check to see if any of the neighboring points satisfy the conditions that will allow this point to change into a LUMEN. At least one neighboring point must be in a POLARIZED or stable CELL that is not already in a TJ and is not the same CELL as the current point. If this is the case the location will calculate the energy change generated by converting from CELL to LUMEN and find the resulting probability p by running the result through an Acceptance Function. If a pseudorandom number r[0,1] is less than p the change will be accepted and the simulation updated accordingly ( Figure S12E). This code allows LUMENS to expand as they would normally without the possibility of CELLS contacting multiple LUMENS. It also allows LUMENS separated by a singlelocation-wide area of CELLS to merge together into a single LUMEN.

Timed shift ISMA
The TS ISMA used an internal clock to determine when a CELL would stabilize. The internal clock was based on the variable shiftCounter, initialized to equal shiftDelay after a CELL POLARIZED. ShiftCounter was decremented at each simulation cycle, and when it reached zero the CELL would change to the stabilized state.

Geometrical mechanism ISMA
The GM ISMA was an earlier version that contained a number of differences in its implementation. The foremost was the method used to determine when CELLS stabilized. In the GM ISMA the following differences existed: • CELLS stabilized when their wedge area was more than twice their actual area instead of relying on LUMEN size.
• CELL clustering was calculated in a different fashion, in which two random number calculations were performed instead of one.
• LUMEN target area was calculated based on lumenGrowthRate, LUMEN perimeter, CELL stretch, and the number of stable CELLS bordering the LUMEN.
• CELLS had an increased likelihood of DYING if they did not contact the MATRIX but did contact stabilized CELLS.

Computational objects
The ISMA consists of a number of computational objects, the most significant of which are listed below: • Point: a grid location

Technical specification
CC3D functions in either 2D or 3D and lets users choose between a square and hexagonal grid. It provides architecture for calculating changes in energy and accepting or rejecting changes in grid locations. The architecture is designed from a system-based perspective. Each simulation cycle, each aspect of the system is executed, from the index change step that selects random points, to the plug-ins that update aspects of the system. The CPM is straightforward and mathematically simple, making it easy to understand and its execution fast. The ISMAs were constructed using CompuCell3D 3.

Potential molecular counterparts to TS and LS ISMA mechanisms
The sirtuin protein SIRT1, a protein deacetylase, has been shown to be involved in cellular senescence, the early phases of which are reversible [2]. This is associated with downregulation of SIRT1 and an increased function of cyclic AMP-regulated kinase (AMPK) [3]. SIRT1 downregulation increases the stability of the AMPK regulator LKB1 and thus inhibits cell proliferation [4]. A possible cell autonomous timing mechanism that could cause cell stabilization might involve downregulation of SIRT1 and upregualtion of LKB1.
Evidence suggests is possible that the tension generated at the luminal membrane is transduced by the subapical F-actin network. It maintains luminal integrity and allows recycling endosomes to aggregate, regulating the protein and lipid composition of the apical plasma membrane. By regulating this F-actin network, cells can control lumen and cyst size.
The Rho family small GTPhase Cdc42 is an ideal candidate for such a regulatory role [5]. Cdc42 is recruited to the apical plasma membrane by the Posphatidylinositol (4,5) bisphosphate (PI(4,5)P2) binding protein Annexin 2, and the loss of Cdc42 disrupts MDCK lumen formation ( [6]. The GTPase exchange factor Tuba activates Cdc42 at the apical membrane, allowing it to control apical exocytosis and thus expand and maintain the apical plasma membrane [7]. Cdc42, in concert with PI(4,5)P2, also helps to polymerize actin by regulating N-Wasp and Arp2/3, thus maintaining the subapical F-actin network [8]. Cdc42 also prevents excessive apical constriction by antagonizing the GTPase RhoA through p190RhoGAP [9,10]. Through the multiple regulatory roles of apical vesicle exocytosis, F-actin scaffold maintenance and regulation of apical constriction, Cdc42 is in a central position to control lumen size.
Cdc42 depletion early in MDCK cystogenesis leads to a loss of central lumen formation, further supporting its central role in this process [6]. Overexpressing a WT Cdc42 in drosophila pupal eye cells increases the apical membrane area [10], possibly reducing luminal tension. This lower tensions and higher surface area allows more fluid to be pumped into the luminal space by transmembrane pumps, increasing lumen size until the apical membrane tension is restored. This molecular mechanism presents one potential explanation for how MDCK cells could sense luminal tension and react accordingly.

Experiments to validate molecular mechanisms behind the TS and LS ISMA
An additional experiment that could validate the LS ISMA would involve disrupting other regulators of luminal tension. Depletion of key regulators of apical tension such as Cdc42 and its partner Protein kinase C zeta/lambda, using RNA interference should decrease lumen size.
Increasing Cdc42 activity, either by overexpression of WT Cdc42 or a constitutively active form (Cdc42-G12V) would be predicted to increase lumen size. In addition, depletion of RhoA and its effectors ROCK I and ROCK II which regulate myosin contractility should also decrease lumen size, since there will be insufficient luminal tension to support lumen expansion. In addition to potentially validating the molecular mechanisms, if disrupting Cdc42 causes cysts to develop multiple lumens and these cysts with smaller lumens fail to stabilize, it will indicate that stabilization is not caused by a timing-dependent mechanism, and is instead caused by a lumen or geometry dependent mechanism.
An ectopic increase in SIRT1 by stable overexpression or knockdown in LKB1 by RNA interference should prevent cysts from stabilizing and promote continued cyst growth. Thus, if cysts grown with cells with overexpression of SIRT1 or LKB1-KD displayed higher cell number at later time points, it would suggest that the molecular mechanism of cell stabilization was linked to SIRT1 and LKB1.