Elucidate growth control mechanisms using reconstructed spatiotemporal distributions of signaling events.

A fundamental biological question is how diverse and complex signaling and patterning is controlled correctly to generate distinct tissues, organs, and body plans, but incorrectly in diseased cells and tissues. Signaling pathways important for growth control have been identified, but to identify the mechanisms their transient and context-dependent interactions encode is more difficult. Currently computational systems biology aims to infer the control mechanisms by investigating quantitative changes of gene expression and protein concentrations, but this inference is difficult in nature. We propose it is desirable to explicitly simulate events and orders of gene regulation and protein interactions, which better elucidate control mechanisms, and report a method and tool with three examples. The Drosophila wing model includes Wnt, PCP, and Hippo pathways and mechanical force, incorporates well-confirmed experimental findings, and generates novel results. The other two examples illustrate the building of three-dimensional and large-scale models. These examples support that reconstructed spatiotemporal distributions of key signaling events help elucidate growth control mechanisms. As biologists pay increasing attention to disordered signaling in diseased cells, to develop new modeling methods and tools for conducting new computational studies is important.

Introduction

How complex and precise developmental signaling and patterning is accurately controlled during embryogenesis in order to generate correct and distinct tissues, organs, and body plans is a fundamental question in biology. The difficulties of examining interactions between genes and proteins under in vivo conditions and unravelling positive and negative feedbacks encoded by these interactions make computational models increasingly used to integrate experimental findings and to unveil control mechanisms. Signaling in a group of connecting cells comprises a series of spatiotemporally ordered gene and protein interactions. Many interactions have been defined in signaling pathways (e.g., Wnt and Notch pathways) [1], but more, including those wrong ones in diseased cells (e.g., cancer cells) [2], remain unclear, because they are context-dependent and show emergent behaviors. It is interesting, important, and challenging to analyze the order and disorder of gene and protein interactions under varied conditions.

Multicellular computational models fall into two classes - lattice models that examine cells with a regular shape in a fixed array and vertex models that inspect cell shape, cell growth, cell division, and tissue growth [3], [4]. Since it is difficult to simulate cell divisions in the lattice model framework and it is inconvenient to solve a large-scale differential equation system in the vertex model framework, modeling platforms such as Chaste have been developed [5], and models built using such platforms may share features of both lattice and vertex models. But, no matter a model is built using a specific method or a platform, inferring the spatiotemporal order of gene and protein interactions from protein concentrations is uneasy and unreliable. It is argued that “the heart of the matter (developmental signaling) is not so much the individual molecules, but more the flow of information and the logic of the system they participate in” [6]; this argument calls for new methods capable of directly exploring the order of gene and protein interactions.

To facilitate describing and exploring the order and disorder of signaling, we have gradually developed a method and tool. The method is to reconstruct signaling events and the tool is a specific programming language. The tool, with cellular automata and object-oriented features (called Cellang++), enables signaling among genes and proteins to be simulated as message passing between objects, and enables cell divisions to be simulated easily in a lattice model framework. Models built using this method and tool indicate that gene and protein interactions follow specific spatiotemporal orders, which help decipher the logic as well as the control mechanisms of signaling [7], [8]. Here this method article reports the method and tool with three new applications.

Methods

The design principle of the method and tool

Cellular automata are a kind of tools for studying complex and dynamic systems [9]. The key feature of such systems (including both physical and biological systems) is that local, simple interactions generate global, complex patterns. It is impressive that simple signaling among cells creates complex tissues and organs. Many researchers assume a natural correspondence between a biological cell and an automata cell and use the latter to model the former. However, the extra features of biological systems demand more flexible and extendable cellular automata-style tools. A tissue or organ contains cells of different types, and these cells can change their type and position and produce new cells. Within each cell, genes and proteins have different attributes and functions, with genes being activated or repressed and proteins being produced or degraded. Thus, local interactions happen within and among cells. Moreover, the interactions within and among cells controlling tissue growth and patterning are not simple; instead, they have rich semantics and show distinct spatiotemporal orders.

To facilitate describing heterogeneous cells and molecules (e.g., genes and proteins), it is advisable to add object-oriented facilities into a cellular automata language, allowing to encapsulate molecules into objects and to describe objects and cells using “local” programs. This also allows interactions between molecules to be simulated as message passing between objects, and allows signaling between cells to be simulated as message passing between cells. By continuously capturing all message passing events in each and every cell, how local interactions generate global and complex growth patterns can be revealed. Thus, we implemented a cellular automata language with object-oriented features to realize the three abstractions: using an automata cell to simulate a biological cell, using an object to simulate a molecule, and using message passing between objects/cells to simulate interactions between molecules/cells. Except several types of sophisticated cells (e.g., neural and muscle cells), many biological cells can be simulated as automata cells if the pattern of tissues and organs, but not the shape of individual cells, is studied. Like the development of many programming languages, we used the C language to realize the cellular automata language, and a model built using the language is first translated into C programs and then complied into an executable file.

Encapsulating genes and molecules into objects

Programmed using this cellular automata language, each model comprises three parts: a Cellang++ program that describes the structure and activities of the cell (the program is shared by all cells) and the molecules within, a text file that describes the cell array and initial input, and two precompiled C files that contain the generic graphic windows. The modeler needs to develop only the first two parts.

The Cellang++ program includes (a) a set of cell fields, (b) a cell program, (c) a message queue (called msgq, implicitly defined by the system), and (4) a set of objects. Each molecule (e.g., a gene or protein) is encapsulated into an object, which includes (a) a set of molecule fields, (b) a molecule program, and (c) a message queue called msgq (implicitly defined by the system) (Fig. 1). Cell fields describe cell attributes (e.g., cell type and age) and molecule fields describe molecule attributes (e.g., protein concentration). Cell and molecule fields are updated in each round of computation. In the cell and molecule programs, variables (which are used to compute cell and molecule fields and to process messages) can be declared. Except the two built-in variables rtime and time (which report the current time and current running step), all variables have only temporary values (they are cleared and re-declared in each round of computation). A set of constants can be declared ahead of cell fields; their values are fixed and most constants are model parameters.

Fig. 1

Using message passing to simulate molecular interactions. (A) An object can receive messages from multiple objects and take actions accordingly, making its activities event-driven. (B) When the Warts concentration (indicated by the field v whose value ranges between 0 and 10) exceeds the parameter para1, Warts sends a Rep message to Yorkie. Warts checks if the message queue msgq contains Act and Rep messages from Yorkie and Dachs. (C) Yorkie checks if msgq contains Rep message from Warts. If Yorkie concentration exceeds the parameter para2, Yorkie sends an Act message to Vestigial. “cell” is used in the messages if the event occurs in the same cell, otherwise a relative address such as [+1,+1] and [-1,+1] is used. Abbreviations: Warts (Wts), Yorkie (Yki), Dachs (D_P, with P indicating the proximal compartment).

In the cell and molecule programs, the sendmsg statement sends messages from one cell/molecule to another cell/molecule, and the if msge $ operation checks messages in an msgq from other cells/molecules. Each message has four components: the target cell, the target molecule, the message, and a parameter that carries quantitative information. For example, if A activates B, A sends an activating message to B, which is captured in B as A_Act_B_0_0 (0_0 indicates B and A are in the same cell) or A_Act_B_p1_n1 (p1_n1 indicates B is in the cell at the relative address [+1,-1]).

Defining message passing between molecules

A message can be defined upon an event (a change of gene state or protein concentration), or upon the occurrence of related events (i.e., the receiving of related messages). If Hill functions are used to describe the nonlinear transcriptional activation/repression and protein activation/repression, the coefficients in Hill functions can be used to define signaling events. When A’s concentration exceeds the half-maximal activation/repression coefficient in the Hill function that describes A nonlinearly activating/repressing B, A sends an activating or repressing message to B.

Displaying signaling events in cells and in the cell space

To display captured signaling events, three sets of graphic windows are developed in the precompiled C files. These windows are opened and closed using hotkeys during simulation. The first set are for debugging models, the second set are for displaying protein concentrations and signaling events in specific cells, and the third set are for displaying protein concentrations and signaling events in the cell space (Supplementary Fig. 1; Fig. 2). The size and names of the first and second sets of windows and the number and names of the third set of windows fit the number and names of fields in the cell and the number and names of fields in molecules automatically, thus these windows are generic. To display these graphic windows, one or multiple color map files should be defined.

Fig. 2

Graphic windows display protein concentrations and signaling events in the 3D cell division model. (A1-A3) The distribution of cells in different planes (the YZ-plane of X = 18 and the XY-plane of Z = 29 and Z = 35) and a time point (T = 158). (B1-C3) Cyclin E and cyclin B concentrations in the same planes and time point (warm and cold colors indicate high and low concentrations). (D-F) The distributions of the E2F1_Act_CycE, E2F1_Act_E2F1, and APCFzr_Ubi_CycA events in the XY-plane of Z = 29 at T = 158 (blue color indicates the occurrence of events). Abbreviations: cyclin E (CycE), cyclin B (CycB), the APC and Fzr complex (APCFzr), ubiquitination (Ubi). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Simulating cell division and cell movement

The copyto([i, j, k], s) statement performs cell division. If [i, j, k] = [0, 0, 0] (i.e., the relative address of the current cell), the daughter cell is inserted into a random position adjacent to the mother cell, otherwise it is inserted into the position specified by the relative address (e.g., [1, 1, 0]). If the specified position is occupied by a cell, the movement of the cell, and possibly also related cells, is performed automatically (Supplementary Fig. 2). In ([i, j, k], s), s = 0 or 1 indicate symmetric or asymmetric cell division (i.e., the daughter cell inherits the same protein concentrations from the mother cell or adopts random protein concentrations). After a cell division, the mother and daughter cells use the current time rtime as their new birth time (this makes the mother cell not immediately dividable). The swapto[i, j, k] statement performs cell movement, moving the current cell at [0, 0, 0] to [i, j, k].

Solving differential equations

If differential equations are defined in objects to compute protein concentrations and gene states, the codes for solving these equations are generated automatically when the Cellang++ program is compiled into the C programs. Partial and ordinary differential equations in all cells are solved simultaneously using the method of lines [10] (in detail, using the Runge-Kutta method with adaptive time steps). In the error control term |err|<|U|·relerr + abserr, the default value for the two error thresholds relerr and abserr is 0.001. As long as |err|>=|U| · relerr + abserr occurs in any equation in any cell, the time step is halved and the current round of computing is repeated in all cells. An initial condition should be defined for each differential equation, and four kinds of boundary conditions (Neumann, Dirichlet, mixed, and periodic) can be defined for partial differential equations. The nonlinearity and stiff of equations, the number of equations, the number of cells, and the value of relerr and abserr together determine a model’s running speed (i.e., steps of numerical solution per second). To run the Drosophila wing model (see Application A) on a personal computer with one CPU (i5-2400, 3.1 GHz), 8G memory, and CentOS 6.0, it takes 1 min to run 627 steps and needs 3–5 h (it takes more time if more graphic windows are opened) to reach T = 210 (T = rtime, and it is non-dimensional).

Built-in functions

Several sets of built-in functions are implemented. First, pos(1), pos(2), and pos(3) return the position of the current cell on the X, Y, and Z axis. These functions are used to define conditions to locate specific cells (a cell that matches the conditions is called the current cell and its relative address is [0, 0, 0]). Second, Hr() and Ha() compute activating and repressive Hill functions. Third, del() computes molecule diffusion in a 2D or 3D space. In addition, all mathematical functions in the C language can be called in the Cellang++ program to perform quantitative computation.

Application A

Background

During embryogenesis cell divisions must be precisely controlled to make tissues and organs reach the correct size and shape [11], [12], [13]. Early studies suggested that growth and patterning is activated by morphogen gradients and that growth stops when these gradients become shallower than a threshold [14], [15]. However, subsequent findings from Drosophila wing and eye development contradicted this simple hypothesis (Supplementary Note 2).

In Drosophila wing development, cell divisions are initiated in the middle of the wing disc along the dorsal–ventral (DV) boundary and growth is mainly along the proximal–distal (PD) axis (Fig. 3A). Key protein interactions controlling wing growth along the PD axis have been characterized (Fig. 3BC). Interactions between these proteins make these proteins form spatial gradients in the wing and biased distributions in each cell. Cell divisions generate a mechanical force, which decreases cytoskeletal tension (especially in the periphery of the wing, i.e., the proximal region, making cells stretched) and increases compression stress (especially in the center of the wing, i.e., the distal region, making cells compressed) [16]. The stress activates Warts (and consequently represses Yorkie) in adjacent cells, downregulating cell divisions [17], [18], [19]. Hence, the core mechanism controlling wing growth along the PD axis includes four elements: Wnt/PCP/Hippo signaling, and mechanical force (Fig. 3C). When and where the four elements interact with each other to initially promote but later stop growth remains poorly understood. Especially, the activation of Warts by the mechanical force created by cell divisions generates a spatiotemporal negative feedback, which should be important for growth arrest. In addition, experiments reveal that the increase of Dachs activates Hippo signaling and that Dachs is degraded by the ubiquitin ligase FbxL7 (Fig. 3C) [20], [21]; but how the increase of Dachs can be reconciled with the degradation of Dachs in the same cells remains unclear. We wish to use a new model to address the two questions.

Fig. 3

Drosophila wing and the wing growth model. (A) Third instar wing disc (left) and ex vivo everted wing (right) (adapted from [50]). Colors show which regions in the wing disc will form which structures in the developed wing. A, P, D, and V indicate anterior, posterior, dorsal, and ventral. The red line in the center of wing pouch but at the distal end of the developed wing indicates the DV boundary. In the developed wing the ventral part of the disc ends up on the back side. (B) A schematic depiction of the Wingless, Vestigial, Four-jointed, and Dachsous gradients along the PD axis. Red arrows indicate Wingless diffusion. (C) A schematic of the model. The a, b, and c shaded areas indicate the Wnt, PCP, and Hippo pathways, respectively. Numbered links indicate experimentally reported protein interactions. Links with an arrow, hammer, or dot indicate activation (e.g., A → B), repression (e.g., A-|B), or binding (e.g., A–B). Dashed links indicate interactions with less experimental support (and not in the primary model). The parentheses around FatDs indicate that the Fat-Dachsous binding is not explicitly modeled. Proteins with a caret (^) have a constant production rate. Proteins that are underlined have a biased (polarized) intracellular distribution along the PD axis, and two differential equations are used to compute its concentration in the distal (shaded blue) and proximal (shaded in pink) compartments. At the two cell junctions D and P indicate the distal and proximal compartments, and big and small fonts of Fat and Dachsous indicate high and low concentrations, respectively. (D) The model’s initial conditions (IC). “#” indicates that Wingless is expressed only in the DV boundary cells. Abbreviations: Wingless (Wg), Frizzled (Fz), Vestigial (Vg), Dachsous (Ds), Four-jointed (Fj), Warts (Wts), Yorkie (Yki), Wg-bound Frizzled (FzB), Dachs-Dachsous (DDs), Fat-Dachsous (FatDs), FbxL7 (F7), Fat-FbxL7 (FatF7), Activating (Act), Repressive (Rep), Stress (Sts), ubiquitination (Ubi). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Lattice models were used to examine Wnt signaling (the diffusion and downstream targets of Wg) and PCP signaling (the gradients of and interactions between Fat, Dachsous, and Four-jointed) in the wing [22], [23], [24]. After it was proposed that cell proliferation is regulated by the mechanical force generated by proliferating cells [25], vertex models focusing on the mechanical force were developed. One model reveals that cells are compressed in the center of Drosophila wing and that growth stops once the compression level in the center reaches a threshold and the compression gradient in the other areas drops below a certain level [26]. After Warts was identified as the molecular target of the mechanical force [19], studies revealed that cell divisions reduce cytoskeletal tension (or increase cell compression) and decrease Yorkie activity [17], and that the spatiotemporal changes of Yorkie activity and mechanical force are correlated [18]. We wish to use reconstructed signaling events to unveil these spatiotemporal changes in more detail, especially, when and where the mechanical force represses Yorkie via Warts.

The model

Using the described method and tool, we built a model that includes the Wnt, PCP, and Hippo pathways and mechanical force. It includes features of lattice and vertex models in a straightforward way: protein diffusion and interaction described by differential equations as in lattice models, and tissue growth simulated by cell divisions as in vertex models. Some simplifications were made to manage the complexity of the model without significantly compromising its power and validity. First, cell division was simulated, but cell growth was not, and cell age was used as a proxy for cell size because cell size does not matter much for control of wing growth. Second, if a protein has a biased (polarized) distribution in a cell (i.e., intracellular distribution) along the PD axis, two equations were used to compute its concentrations in the proximal and distal compartments in each cell; this implicitly divided a cell into two compartments. Third, all proteins were assumed to have the same half-lives.

Model parameters were set based upon experimental findings, constraints among interacting proteins, and tuning the model to produce observed phenotypes. Since Vestigial and Yorkie are required for cell division [27], and cell divisions occur unsynchronized and stochastically [28], four conditions (the concentration of Vestigial, the concentration of Yorkie, the age of cells, and a random number that determines the chance a dividable cell divides. As a cell divides, the daughter cell was inserted proximal to this cell and inherited the protein concentrations from this cell (Supplementary Fig. 2). Experimentally identified protein interactions were realized, and model robustness was tested by changing individual parameters and by performing parameter sampling (Supplementary Table 2, 3).

Results

The model reproduced many experimentally observed cell and protein activities (Supplementary Note 4) [29], [30], [31], [32], [33], [34], [35], [36], [37], including the finding that cell divisions occur more frequently in the proximal regions at late developmental stages. As long as a dividing cell generated a stress (i.e., the mechanical force, hereafter also called cell stress, with the parameter stresspara > 0.0) on its neighbors, growth arrest occurred. A key feature of Drosophila wing growth is that even when the anterior and posterior regions have different cell division rates, a wing of normal size and shape is still generated [38]. We examined whether the model could recapitulate this subtle phenotype by running the model with cells having different cell cycle lengths in the anterior and posterior regions. Under all parameter settings the model eventually generated the same cell numbers in the anterior and posterior regions (Supplementary Fig. 4). If the parameter stresspara = 0, unequal numbers of cells were produced in the two regions. These results provide confidence that the model is accurately reflecting the developmental process.

By examining the spatiotemporal distribution of Yorkie repression, we next explored when and where cell divisions are repressed and obtained the following findings. At the early stage the high Wg and Vestigial concentrations in the center of the wing cause increased cells, activated Warts, and densely distributed Warts_Rep_Yorkie event (Fig. 4A-E). As the wing grows, as experimentally observed [16], activated Vestigial in proximal wing cells regulates Four-jointed and Dachsous expression, causing the sharp slope of the Four-jointed and Dachsous gradients gradually shift proximally (Fig. 3B; Supplementary Fig. 3). The repression of Yorkie by Warts (activated by cell stress in the distal (central) region of the wing) and the activation of Dachs (activated by Fat-Dachsous interaction in the proximal (periphery) region of the wing) make the cell divisions reduced in the center and increased in the periphery. However, as Warts is repressed by Dachs continuously, Warts causes fewer Warts_Rep_Yorkie in the center. Meanwhile, the gradual activation of Warts by increased cell divisions in the periphery increases Warts_Rep_Yorkie. Together, these two regulations make the densely distributed Warts_Rep_Yorkie event shift from the center to the periphery (Fig. 4A’’-E’’). We thus propose that the densely distributed Warts_Rep_Yorkie in the periphery indicates growth arrest. To examine whether this distribution would occur robustly under varied conditions, parameter sampling (to change multiple parameters simultaneously and examine results influenced by these changes) was performed (Supplementary Table 3). We found that whenever wing growth occurs, the densely distributed Warts_Rep_Yorkie in the periphery is generated (Fig. 5), and when in several cases Warts_Rep_Yorkie occurs in the whole wing, Yorkie repression by Warts is still more prominent in the periphery (Fig. 5T2, W2, X2), indicating that Warts_Rep_Yorkie in the periphery is a reasonable indicator of growth control.

Fig. 4

The spatial distributions of protein concentrations and signaling events at different time points during wing growth. In all panels the vertical bar indicates the DV boundary. (A-A’’) show the experimentally observed growth process (from [16] with permission). In (B-C’’) warm and cold colors indicate high and low concentrations of Warts and Yorkie. In (D-D’’) blue and white sites indicate the presence and absence of the Warts_Rep_Yorkie event in cells, and a ring-shape dense distribution of the event in the periphery is apparent at T = 210. (E-E’’) show cell division rates. As experimentally observed, Yorkie concentration varies less significantly than Warts concentration does, and cell division rates are initially high (light blue), then decline (deep blue). The parameter settings were agepara = 10, randpara = 0.99, stresspara = 0.04 for controlling cell age, random cell divisions, and stress. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5

Results of parameter sampling. Panels O2 and O3 (showing the distribution of Warts and Yorkie concentration) indicate that when protein concentrations differ limited it is uneasy to infer signaling events. (C-E, JK, Y) show the distribution of cell divisions, in these cases the parameter changes failed to make cell divisions. In all other cases the Warts_Rep_Yorkie event is densely distributed either globally (in panels L, N, S, T1, W1, and X1) or in the periphery. (T2, W2, and X2) show the distributions of Yorkie concentration, which indicate that Yorkie repression by Warts is more prominent in the periphery. Parameter changes (Supplementary Table 3) and phenotypes are compared with the default ones (Supplementary Table 2; Fig. 4).

Since Dachs indirectly activates Yorkie and Yorkie stimulates cell division (Fig. 3C), it would be rational to observe high Dachs concentration in proliferating cells [39], [40]. Both our simulations and experimental studies generate this observation. However, it was recently reported that Dachs is degraded by FbxL7 in wing cells [20], [21]; we thus examined how the degradation of Dachs can be reconciled with the increase of Dachs. Under all parameter settings, if Dachs has a biased distribution and is degraded only in the proximal compartment in each cell, Dachs concentration is slightly increased in the whole cell (Supplementary Fig. 5). Experimental studies did not reveal whether Fat is involved in the degradation of Dachs and Dachsous [20], [21]. We found that, if we let both FbxL7 and the Fat-FbxL7 complex degrade Dachs, Dachsous, and the Dachs-Dachsous complex (Fig. 3C), the distributions of all these proteins agree with experimental findings, but otherwise not. This is an experimentally testable prediction that both FbxL7 and its complex Fat-FbxL7 degrade Dachs, Dachsous, and the Dachs-Dachsous complex.

Discussion

Upon the experimental findings that cells are compressed in the center, that the stress in compressed cells activates Warts, and that Warts represses cell division by repressing Yorkie, it is hypothesized that growth arrest is initiated in the center of the wing. But this hypothesis cannot explain how growth is repressed later in the periphery [26]. The reconstructed spatiotemporal distribution of the Warts_Rep_Yorkie event instead suggests that growth arrest is initiated in the periphery. We also examined another hypothesis, which is that the global Four-jointed and Dachsous gradients will eventually become too shallow to activate Dachs and Yorkie, thus cause growth arrest globally [41], [42]. The model indicates that the changes of Dachsous and Four-jointed gradients are inadequate to cause growth arrest (Supplementary Note 7).

In a growing tissue or tumor, cell divisions anywhere activated by any growth-promoting pathway can activate growth, but growth arrest is coordinated globally, making it difficult to study growth arrest by examining only an isolated part of the tissue or a subset of the relevant pathways. Thus, the mechanisms controlling growth arrest are less studied and less understood. This wing growth model recapitulates multiple experimental findings [16], [20], [21], [35], [38], [40], [43], and generates two novel findings. First, growth arrest starts in the periphery of the wing, instead of in the center as has been hypothesized [26], and second, Fat is involved in the degradation of Dachs and Dachsous. These results suggest that the distributions of protein concentrations and signaling events together better unveil when and where growth arrest is initiated.

Application B

The growth and patterning of most tissues and organs depends on the orderly control of cell divisions in the 3D space. It remains poorly known how different crosstalk among these evolutionarily conserved signaling pathways generates various signaling cascades and spatiotemporally ordered cell divisions. It is important to model cell divisions in 3D space to address the question.

It is easy to use this tool to build 3D models. To indicate this easiness, here we show how to transform a Drosophila cell cycle model into a 3D cell division model. Most cell cycle models are single-cell models and cell divisions are not simulated, but the one we developed using this tool is a 2D model (Supplementary Fig. 6) [8]. We first extended the cell cycle model into a 2D cell division model by adding the condition that controls cell divisions and the copyto statement that implements cell divisions, then extended the 2D cell division model into a 3D one. The 2D and 3D form of the cell division model has only two differences: the parameter in the dimension definition and the parameter in the copyto statement (Supplementary Fig. 7). By pressing the hotkey “x”, “y” or “z” one can choose to access the YZ-, XZ-, and XY-planes of the 3D model, and by further pressing the hotkey “+” and “-” one can go up or down along the chosen axis to check and capture protein concentrations and signaling events in specific planes (Fig. 2). This cell division model can be easily integrated with different signaling pathways (Supplementary Fig. 8), facilitating the investigation of growth and patterning of tissues, organs, and tumors.

Application C

To simulate cell divisions in more detail during the growth of normal tissues and organs and tumors, it is necessary to incorporate a cell cycle model into a tissue growth model containing multiple signaling pathways. As an example, we used the tool to integrate the 2D cell division model into the 2D wing growth model (Fig. 6A; Supplementary Fig. 6), generating a highly biologically detailed model for examining the orderly control of cell divisions by the Wnt, PCP, and Hippo pathways and the cell cycle pathway in the wing. The two models are bridged by the activation of cyclin E (a conserved protein that activates cell divisions in multicellular organisms) by Yorkie [44], thus we used Yorkie to replace the growth factor in the equation of cyclin E to implement the Yorkie-activated cyclin E expression. The captured protein concentrations and signaling events reveal how the spatiotemporally controlled cell divisions lead to wing growth and patterning (Fig. 6B-F). This model can be further extended to include more genes and pathways.

Fig. 6

The expression of cyclin E and cyclin A and the distribution of the Skp2_Ubi_CycE (Skp2 ubiquitinating cyclin E) event in the developing wing. The high expression of cyclin E (C) and subsequently cyclin A (D) in the periphery of wing disc agrees with the experimentally observed high expression of Yorkie (A) in the late stages (panel A is from [18]).

Concluding remarks

Experimental studies have revealed that the development of distinct tissues, organs, and body plans is controlled by a set of highly conserved genes and pathways. This finding not only supports the argument that the heart of developmental signaling and patterning is more about the flow of information and the logic of the system these genes and pathways participate in [6], but also implies that these genes and pathways critically control the growth of various organoids [45], [46]. The difficulties of continuously detecting transient protein interaction and gene regulation changes under in vivo conditions make computational studies increasingly important. Although multiple modeling tools and platforms (such as Chaste [5]) have been developed, to build complex models is still challenging. During a signaling process, protein concentrations may change insignificantly (especially in non-dimensional models), a gene or protein may be regulated by multiple regulators simultaneously, and it is difficult to differentiate the causes from the consequences [47]. These make it difficult to uncover the order and disorder of signaling using purely quantitative models. As biologists pay increasing attention to disordered signaling in diseased cells, to develop new modeling methods and tools for new computational studies is needed.

Biologically detailed models sharing the features of lattice and vertex models are beneficial; for example, they help examine the role of Wnt diffusion during wing growth (Supplementary Note 4). This method and tool not only makes it easy to build detailed models, but also enables reconstructing signaling events. In addition to revealing the order and logic of normal signaling, reconstructed signaling events in biologically detailed models can help reveal how to rectify wrong signaling in diseased cells (Supplementary Fig. 8) [2], [48]. For researches using organoids to study brain development, tissue repair, and tumor growth [45], [46], reconstructed signaling events are especially useful for identifying mechanisms that can control desired growth and patterning. Even if reconstructed spatiotemporal distributions of signaling events may not be readily validated experimentally, they are valuable for guiding future experiments. If computed protein concentrations agree widely with experimental findings and if subtle phenotypes are reproduced, one can have confidence in the simulated signaling events.

We note that the described tool is not necessarily the sole implementation of the proposed method. The tool is suitable for building multicellular models but not single-cell ones, for building biological detailed models but not concise ones, and for modeling cells without shape and size but not cells with complex shape and different size. Cells with two compartments can be conveniently described, and hexagonal epithelial cells can be described by dividing a cell into six computational units with a specific cell neighborhood [49]. Differential equations are not indispensable, because messages that carry quantitative information can describe complex molecular interactions. Given that quantitative models are computationally costly and qualitative models may be incompetent, models based on message passing with quantitative information in messages provide a choice for the third approach.

CRediT authorship contribution statement

Hao Zhu: Conceptualization, Data curation.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.