Developmental Engineering: A New Paradigm for the Design and Manufacturing of Cell-Based Products. Part II. From Genes to Networks: Tissue Engineering from the Viewpoint of Systems Biology and Network Science

Structure of the article

In section 2, we introduce systems biology. The need for the use of its concepts and methods arises from the fact that cell, tissue, organ, or organism functions can not only be attributed solely to single genes or linear cascades of gene activation but also require extensive crosstalk and feedback of signaling pathways that form networks of interacting genes. Because of the high complexity of the interactions in comparison with linear cascades, it is not possible to describe the network dynamics only by intuition, because the changes caused by the input (e.g., addition of growth factor) are spread throughout the whole network omnidirectionally and iteratively until the gene expression stabilizes. Therefore, the concepts and mathematical methods used in systems biology to decipher what the gene/protein network does and to what cell/tissue function its operation corresponds are of great relevance to tissue engineering.

In section 3, we present an example from the literature of developmental biology (segment polarity pattern formation in Drosophila), in which systems biology methods were used to determine the gene/protein network responsible for the formation of the segment polarity pattern, a macroscopic developmental module at the tissue level. This example not only clarifies how systems biology is used in development for the determination of modularity of tissue forms but also indicates how it should be used to design in vitro biomimetic processes that simulate development in vivo. The important issue here is that the macroscopic modularity of tissue forms is attributed to autonomous (no external factors involved) operation of a gene/protein network that is equally modular. In other words, the macroscopic modularity of the engineered tissue intermediate results from corresponding modularity of a set of genes that define that state. This set of genes is relatively isolated from the rest of the genes/proteins of the overall network operating in the cell, but the connections between its member genes/proteins are strong and not easily affected by perturbation of its parameters and initial conditions.

In section 4, we use the methodology of systems biology presented in section 3 to design a biomimetic process for growth plate development, which we presented briefly in Part I. The purpose of the design is to determine the initial conditions and parameters of the process, so that the gene/protein modular network responsible for the macroscopic modularity of the columnar pattern of the growth plate will be activated, thereby establishing the spatially differential gene expression pattern observed in the growth plate. We will confine discussion of mathematical modeling of gene/protein interactions to the interactions observed in vivo during development, and thus restrict to processes that are biomimetic by default. In general, these concepts do not apply to one-step concerted manufacturing schemes, which therefore will not be examined further.

In section 5, we deal with the problem of having several stable states, instead of one, in which a gene network can settle as it is usual in the process of cell differentiation. The robustness of the gene/protein network, that is, its ability to give the correct gene expression observed in differentiated cells despite environmental disturbances, is reflected in the properties of the mathematical model, which gives the same correct solution despite changes in the initial conditions, as mentioned in section 2. This solution represents a stable state, called an attractor because it “attracts” initial conditions that become finalized through network dynamics. Such a representation is useful, especially when several stable states coexist. Which one will be realized at the end depends on the initial conditions that could direct cells to any of the several different stable states. Several examples are given from the literature referring to cell differentiation, where this representation has proven useful in organizing and explaining experimental observations. It is again the stability of the states and the determination of the initial conditions leading to each one that can solve the problem of process design in an accurate computational way.

In section 6, we introduce several concepts from the science of networks and show how they are important in process design. Network science, as systems biology, deals with networks. However, although systems biology focuses on component, for example, gene or protein, interactions in networks from the point of view of regulatory mechanisms as described in section 2, network science is a new scientific field that examines the topology of the networks in a more abstract form, trying to decipher common organizational principles among diverse networks. Work done to date suggests strongly that the mathematical behavior of networks will be instrumental in providing tissue engineering with a solid theoretical background comparable to that of other engineering fields. In any case, design questions could be answered with both systems biology and network science. For example, the mathematical model presented in section 2 from the point of view of systems biology, which includes several regulatory mechanisms in the details of protein-to-protein/gene interactions, has a corresponding abstract model consisting only of gene-to-gene interactions described exclusively by the network topology (see section 2).

In section 7, we present as an example an in vitro biomimetic process of pancreatic induction in endodermal cells by mesoderm, where we try to answer the process design questions from the point of view of network science.

Systems biology: From component interactions to systems behavior

The genocentric paradigm of biology, which placed the gene and its function as primary in biological studies, has provided an enormous amount of data concerning the individual cell components and their functions. However, only limited information about functions can be extracted directly from the genome. Biological functions cannot be attributed to individual molecules but rather to complex interactions between the numerous components of a cell for cell functions or between cells for tissue functions and so on, spanning the levels of hierarchical organization of organisms. For example, most diseases are not caused by a single gene defect but rather by a malfunctioning network of interacting genes and their coded proteins. More than 100 genes have been identified as contributing to the coronary artery disease. ⁵ It is therefore apparent that instead of dealing with single genes, in this case we instead have to generate information concerning the behavior of these genes in an ensemble of 100 interacting genes and proteins in a functional network. Even collections consisting of a small handful of components may display a behavior markedly different from those of the individual components. A very simple example is the markedly cooperative binding of oxygen to a system of four hemoglobin subunits, starkly different qualitatively from the noncooperative binding displayed by myoglobin. Of more direct relevance to cell fate decisions taking place during development, a cascade of three MAP kinases displays stimulus–response characteristics profoundly different from those expected for a single component. ⁶ A model composed of a set of ordinary differential equations based on accepted principles of enzyme kinetics, ⁷ with parameters that may be determined by experiment or estimated, predicts that in contrast to a graded stimulus–response relationship for a single component, the cascade exhibits switch-like, “all off/all on” behavior.

Currently, there is a gradual emergence of a new paradigm, which treats biological phenomena from the systems point of view with a bottom-to-up approach, trying to decipher the system properties from the properties and interactions of the components, ⁸ instead of analyzing the system to its components. The path for this change has already led to the development of systems biology. Systems biology seeks to understand how all the individual components of a biological system interact in time and space to determine the function of the system, be it cell, tissue, organ, or organism. It makes use of the large amount of data from molecular biology and genomics to develop mathematical models of the complex function of such systems. Systems biology will change research practice and lead to the integration of information from the molecular up to the organism level. For example, in experiments designed to elucidate the underlying pathophysiology of a disease, data collection and interpretation are of equally important. For complete interpretation in terms of physiology to be achieved, the use of mathematical models to integrate huge amounts of data describing gene expression, protein function, cell function, and whole-body physiology are needed. ⁹ Clinical efficacy of drugs can be also predicted using physiological models of disease and disease processes. ¹⁰ Several applications in health have been already published. Gadkar et al. have developed a mathematical model of the pathogenesis of type 1 diabetes, and they used it to study the effects of anti-CD40L therapy, determining optimal treatment protocols. ¹¹ Rullmann et al. have developed a mathematical model to describe the inflammatory and erosive processes in afflicted joints of people suffering from rheumatoid arthritis, including in this several processes such as the life cycle of inflammatory cells, endothelium, synovial fibroblasts, and chondrocytes, as well as their products and interactions, since it is actually the interplay between these processes that determines the clinically relevant measures for inflammation and erosion. ¹² The authors used the model to predict the therapeutic effect of modulating several molecular targets.

From the above examples, it becomes evident that the methods of systems biology are relevant when we have to integrate information at one level (e.g., gene expression), to find answers to questions referring to a higher level (e.g., pathophysiology at the level of organism or drug effect to the patient). There are many such questions in biology and medicine that are now being approached with systems biology methods. ¹³

Systems biology in tissue (developmental) engineering

In the case of tissue engineering, such questions are also critical for the process design, since the essence of the task is to find in vitro conditions in which the integration of gene and cell interactions will lead to differentiated cells that are in a stable state, or bioartificial tissues that are functionally integrated and robust entities. A random distribution of cells in scaffolds, even if cell viability is retained by the sophisticated methods/tools of tissue engineering, is much less likely to closely resemble an authentic living tissue that exhibits properties that arise from sequential cell interactions that occur during natural development and are qualitatively different. ¹⁴ Internal integration of developmental modules through cellular interactions that makes this cell collection behave as a distinct unit signifies that these entities are autonomous (see Part I). Such a living entity has distinct properties, which are qualitatively different from the properties of its component cells. For example, the control of growth plate elongation is not a chondrocyte property but a property of the growth plate module arising from the interaction between chondrocytes participating in the negative feedback loop of Ihh/PTHrP. ¹⁵ Similarly to the intermediate modular tissue forms, tissues in their final form constitute integral entities with properties arising from interactions between their cells. For example, glucose homeostasis in the liver is a function of the liver as a whole, not of isolated hepatocytes, which emerges from the metabolic cooperation of glycolytic (periportal) hepatocytes, which take up glucose during the absorptive phase, and gluconeogenic (perivenous) hepatocytes, which release glucose during the postabsorptive phase. ¹⁶ Another example is the controlled release of insulin, which is not a function of beta cells but a function of the organized beta cells within the islet structure.^17,18 Perhaps the clearest example is the neural tissue; its signal processing functions are based on the topology of the synaptic network instead of on single neurons. ¹⁹

Systems biology in development

In development, systems biology aims to extract the general design rules of the network of interactions of genes, proteins, and/or cells that are responsible for integrating the behavior of components—genes, proteins, or cells into the system, the most important property being the modularity/robustness of cell states or intermediate tissue forms. ²⁰ This effort has already provided important information allowing detection of general architectural characteristics of the networks in development, such as the feedback loops that ensure the progression of development and the repressing interactions that participate in spatial control. ²¹ As numerous recent studies show (e.g., Refs.^22–28 ), it is clear that systems biology has already gained wide acceptance by developmental biologists.

For many signaling pathways controlling cell specification, sufficient information now exists to construct and test models built from several individual pathways. These models open the door to a rigorously scientific, quantitative description of modular behavior observed in empirical experiments by developmental biologists for decades. This will allow direct experimental evaluation of modular properties in a tissue engineered in vitro, instead of approaches relying only on macroscopic phenomena, which might require extensive experimental work. This is no longer just a theoretical possibility: robust models have predicted behavior of systems as complex as the developing fly and frog embryo reliably^23,24; these will be discussed further. This suggests strongly that existing technology can model individual developmental modules determined to be needed in a given manufacturing scheme to guide the nature and extent of externally imposed controls and also set limits for measurable parameters consistent with process design objectives.

Below, we will describe how systems biology is applied to development in vivo, what kinds of questions can be answered and what methods are used (section 2), and then how to transfer these concepts and methods to the design of in vitro biomimetic processes (section 3), thereby transforming the questions of section 2 to process design questions.

The example presented in section 3, selected because of its simplicity, makes clear how a multicellular system can be robust because of the cell-to-cell interactions that maintain intracellular gene/protein interaction networks leading to spatially differentiated gene expression. This is the aspect of development most relevant to developmental engineering. The model presented is based on mathematics no more complicated than ordinary differential equations. This makes the incorporation of process design parameters, such as initial cell concentrations and cell-to-cell interaction through diffusing signaling proteins, a fairly straightforward exercise.

The segment polarity pattern of Drosophila is formed by the operation of a modular, robust gene/protein network

A particularly instructive example of a developmental module where such a formulation of modularity/robustness was achieved is the segmentation appearing during development of the fruit fly, Drosophila, one of the best understood developmental mechanisms. It will be apparent through the analysis of this example that macroscopic modularity of intermediate developmental multicellular forms observed empirically cannot be attributed to single genes or signaling pathways acting as separate entities. Modularity, as a global property of the spatially extended biological system, arises instead from the way members of a particular set of genes influence expression of other genes in the set, forming a complex network of mutual interactions. The gene expression network of one cell extends its action to neighboring cells, influencing their gene expression and activating various signaling pathways through secreted proteins. In turn, these cells respond to the first, activating its signaling pathways in a specific way, so that finally the gene expression is stabilized. The gene interaction network therefore spans the system extended and coordinated throughout the whole macroscopic developmental module, and determines which genes will be expressed and in what locations in the developing organism. Just as the developmental module is robust macroscopically, the same robustness is exhibited by the gene/protein interaction network in its operation autonomously, keeping active the interactions of its components and stabilizing the spatially differentiated gene expression pattern. Transferring the microscopic robustness of the gene and signal pathway network to macroscopic robustness observed experimentally in this way connects measurable/calculable variables referring to the gene/protein network to phenomena that can be observed directly. Thus, one can finally connect process design parameters with the robustness of the tissue form using a mathematical model and predict the necessary modifications in the in vitro process.

The various parts of the body of insects develop on particular segments, layers of cells that appear during embryogenesis. The genes expressed in each segment specify the correct number of body parts and the correct polarity of each one. In Drosophila, a complex network of gene interactions converts a single-celled Drosophila egg to a multicellular embryo with 14 segments, forming a spatial pattern of parallel stripes with each segment bounded by a stripe of cells expressing the engrailed gene, en (Fig. 1a). Various sets of genes are expressed in space differentially, in consecutive stages of development, before the final pattern of 14 segments appears. At each successive developmental stage, the pattern of differential gene expression becomes more precise, with the expression of genes at any given stage controlled by toses expressed in the previous stage This sequential pattern of events, characteristic of developing systems and also of mathematical models, that describe them (see below) is known as path dependence. Most of the genes that are expressed successively are transient. The segment polarity genes, the last ones in the cascade of spatially differential gene expression, are expressed permanently (Fig. 1b). These genes are activated in a spatially differentiated pattern that leads to 14 segments by the five pair-rule genes (even-skipped, hairy, odd-skipped, paired, and runt), which in the previous stage have been also expressed forming a prepattern of seven stripes, each one corresponding to a pair of the final segments.

OPEN IN VIEWER

FIG. 1.

The segment polarity pattern of Drosophila, a macroscopic developmental robust pattern, composed of 14 segments forming a spatial pattern of parallel stripes with each segment delimited by a stripe of cells expressing the engrailed gene, en (a). The segment polarity genes (en, wg, ptc, cid, and hh), the last in the cascade of spatially differentiated gene expression during Drosophila development, are expressed permanently and differentially in each cell of the segment, with sharp boundaries between segments implying that the segmental pattern is a macroscopic multicellular developmental module with robust regulation (b). The assembly of several gene interactions in adjacent cells that communicate through secreted wingless and hedgehog proteins form a spatially extended gene network that gives the pattern robustness (c). CID, cubitus interruptus; CN, repressor fragment of cubitus interruptus; en, engrailed; hh, hedgehog; PH, patched–hedgehog complex; ptc, patched; wg, wingless; small letters correspond to genes and capitals to proteins.

The sharp boundaries between segments, corresponding to cell layers expressing or not expressing particular genes, imply that the segmental pattern is a macroscopic multicellular developmental module with robust regulation. It is known that once it is triggered, it is self-maintained, or autonomous, ²⁵ which is not the case for prepatterned pair-rule genes that do not persist.

Mathematical model

The question from the point of view of systems biology is what particular genes are expressed in each location in the pattern or why these genes are expressed there. Thanks to the studies of Drosophila developmental biology, all these genes, their spatial expression patterns, and many of their interactions are known. What systems biology seeks to identify is the regulatory mechanism(s) responsible for the robustness of the segment polarity pattern. In other words, the question refers to a system property, robustness, and not to a component property, gene expression or gene regulation. This implies that the gene/protein network, that is, which genes or proteins and how they interact, is known to some extent. If we introduce a further criterion of robustness to the information already gathered from experimentation, we can identify still unknown interactions in the mathematical model of the network, as we will see below. This mechanism cannot be attributed to one gene, any of the dual interactions between genes, or any single signaling pathway, since the system is spatially distributed through cells that express different genes (which, however, communicate), determining each others gene expression. It is rather the assembly of several gene interactions in a spatially extended gene network that confers robustness to the pattern; in other words this network, by itself, is robust. As such, systems biology has to synthesize the various observations of molecular genetic analysis and see how these observations fit together to establish the system's properties. Von Dassow et al. have made that step by modeling the activated signaling pathways and relevant genes of the interaction between cells that secrete Wingless (wg) and adjusted cells that secrete Hedgehog (hh, Fig. 1c). ²⁶ The authors developed a mathematical model of 136 ordinary differential equations for the core network comprised of five genes (en, wg, ptc [patched], cid [cubitus interruptus], and hh) and their proteins. The equations describe the time evolution of the concentration of mRNA, protein, or protein complexes, and have terms for synthesis, decay, transformation, and transport. The authors incorporated the known gene interactions, assembled into several small-scale regulatory mechanisms such as positive and negative feedback loops, and are subsequently assembled in one integral network. ²⁷

The initial model could not predict the macroscopic characteristics of the segmental pattern for any value of the parameters—such as half-lives of mRNA and proteins, binding rates, cooperativity coefficients, and diffusion constants—or initial conditions, such as a prepattern specified by the pair-rule genes (expressed during the developmental stage immediately before segmentation). However, the incorporation of two additional gene interactions (wg auto-activation and the inhibition of en by CN, the repressor fragment of cubicus interruptus protein) gave a surprisingly robust model able to predict the correct pattern over a large range of parameters and initial conditions. This does not mean that the range of parameter values allowing the model to give the correct solution pattern is infinite. For randomly chosen parameter values, the authors observed the correct pattern in 1 in every 200. This is not small, since there is a 90% chance for a randomly chosen parameter, among the 48 involved, to be compatible with the existence of a solution (0.9⁴⁸ is approximately 1/200). In vivo, the variations of parameters can be thought of as analogous to mutations of small effect; variations of initial conditions correspond to developmental noise. It is interesting that there are no optimal values of parameters, since almost every parameter in the model can range over orders of magnitude and still give the correct solution. Thus, it is primarily the organization of the gene network that provides the stability of the model and not the details of molecular interactions. Von Dassow and Odell characterized this network as a module, “a device unto itself” because it can accomplish the task of maintaining the spatial gene expression pattern “without any persistent, extrinsic spatial, or temporal biases on any of its components.” ²⁷

Boolean networks represent clearly the behavior of gene networks

To clarify the concepts above, we will use a simple model as an example. The simplest modeling approach to simulate regulatory systems employs Boolean networks, first proposed for gene interactions by Kauffman.^28,29 In Boolean networks, the gene expression levels can be only in one of two states: either 1 (ON) or 0 (OFF). The inhibitory or stimulatory effect of each gene to the other is represented by arrows connecting the nodes/genes. The network behavior over time is modeled as a sequence of discrete synchronous steps. The value, 0 or 1, of a node–gene at the next time step depends on values of other genes that influence it. In Figure 2, we see a simple example of three nodes/genes, A, B, and C, which interact to influence each other's expression, thus forming a small network. Gene A is expressed when it is induced by gene B; gene B is expressed when it is induced by both A and C simultaneously; and gene C is expressed unless A is expressed, in which case expression of C is inhibited. If we assume that at this time the state of the network is presented by the vector (0, 0, 0), which means that none of the genes has been activated, at the next time step the vector will take the value (0, 0, 1) (Fig. 2b). This is because since A is not expressed there is no inhibition of C, which therefore could be expressed (taking the value 1). However, B still cannot be expressed and its value remains 0 because its expression requires induction by both A and C and A is still missing (i.e., not expressed). After reaching the state (0, 0, 1), no further changes can be made; therefore, this state can be considered as a solution. Since there is no A, C keeps the value 1 and B remains with the value 0 because there is no A. In the case of starting from another initial condition (1, 1, 1), the network again adopts (0, 0, 1) as final state, but in this case, it passes through more intermediate states. As we see, this is a small “device unto itself,” because its state at each moment, that is, which genes are expressed or not, depends only on these genes and their interactions. Different initial conditions provide different solutions, defining different final states. The final states are such that the next step gives the same state, and therefore we could say that the network operates after this point to reproduce the same states (or oscillate between two states, Fig. 2b, between (1, 0, 1) and (0, 1, 0)). If we assume that the state (0, 0, 1) corresponds to the correct solution, corresponding for example to an experimentally observed pattern or cell differentiation state, the developmental process has to assure that in the previous stages the network will not end in the states (1, 0, 1) or (0, 1, 0), which oscillate between themselves. Therefore, the role of prepatterning is to assure that correct initial conditions have been set up to ensure that the next stage the network reaches is the correct state. We could also express this by saying that the network is triggered to reach the correct state solution by certain sets of initial conditions but not by others. We could also say that the solution (0, 0, 1) is robust, while the (1, 0, 1) that oscillates with the (0, 1, 0) unstable because in the latter it is only two of the possibly eight initial conditions that realize it and any disturbance will lead to the solution (0, 0, 1).

OPEN IN VIEWER

FIG. 2.

A simple example of a Boolean network of three genes, A, B, and C, that interact, influencing each other's expression and forming a small network. Gene expression levels can be in only one of two states: either 1 (ON) or 0 (OFF). Gene A is expressed when it is induced by gene B; gene B is expressed when it is induced by both A and C simultaneously; and gene C is expressed unless A is expressed which inhibits expression of C (a). Different initial conditions such as (0, 0, 0), (1, 1, 1), or (1, 0, 1) may lead to the same or different final states such as (0, 0, 1) and the oscillating condition [(1, 0, 1), (0, 1, 0)], passing from different states that form a trajectory (b).

Though the Boolean networks are simplifications of the real situation, they clearly show the importance of network topology, that is, which genes influence the expression of which others and how. It is with these types of models that Kauffman has shown that constraints arising from the network topology, such as the gene inhibition, make large networks settle down to a finite number of solutions.^29,30 For example, if gene A unconditionally inhibits gene B, then there would be no possibility to have as a final state one in which A and B are both present. In a network that includes all the 25,000 human genes, theoretically 2²⁵⁰⁰⁰, that is, around 10⁷⁵²⁶, states are possible, which is an astronomical number. If we consider each final state as representing the gene expression of a cell state, then there would be almost a continuum of states instead of the approximately 200 that are observed. ³¹ This shows the importance of the topology of the network in restricting the possible solutions to the number consistent with observation.

Network topology, an important factor of modularity

Albert and Othmer have used a Boolean model for the Drosophila segment polarity gene/protein network based on a binary ON/OFF representation of mRNA and proteins, that is, the concentrations were not continuous variables. ³² The authors showed that the dynamics of the gene interactions are determined mostly by the network topology that is, which gene interacts with which other(s) and their type of regulation (e.g., inhibiting or activating), without details of the rate laws. Since this model does not include any kind of functional details and has only two states for the variables, but is nevertheless still able to reproduce the segment polarity network, it is clear that the topology of the network and the presence or absence of its components confer robustness, and not the functional details or the absolute values of component concentrations. As the authors mentioned, the fact that the segment polarity gene network could be successfully modeled by a simple Boolean model does not mean that other networks could be also modeled in such degree of simplification. They might require more detailed models, as for example incorporating asynchronous updating or multilevel instead of binary variables, or in some cases as in metabolic networks the Boolean approach might not be appropriate. In addition, they mentioned that the choice of Boolean model for the segment polarity network was motivated by the ON/OFF character of the experimentally observed gene expression patterns.

The conclusions drawn from the above presented examples relevant for the biomimetic process design are as follows:

The above studies show clearly that the underlying cause of the macroscopic developmental modularity can be attributed to the corresponding modularity of a network of interacting genes/proteins.

Systems biology can accurately answer difficult process design questions

With the attribution of the macroscopic modularity of intermediate tissue forms to the operation of a gene/protein network that is self-sustained, and thus which network is a relatively isolated subset of the whole network operating in the cells, that is, a module, several obscure issues of biomimetic process design can be illuminated. The degree of match between the in vitro and the corresponding in vivo process, questioned in literature, ² that can assure that the process is really biomimetic comes out directly. If the in vitro process succeeds to activate the gene/protein module, then it has achieved the degree required since the goal of the process is the formation of intermediate modular tissue forms, like the segmental pattern of Drosophila, resulting from the operation of the gene/protein module. Von Dassow and Odell have examined several simpler versions of the model presented in the previous section, not including all the gene/protein interactions, and they found that although some of these models may give the correct solution pattern, the solution becomes less frequent in the parameter range, that is, less parameter values were able to provide the solution. ²⁷ This means that the subset of activated interactions, as well as the resulting pattern, is less robust from those of the fully activated network. If experimental conditions in the in vitro system activate only a subset of the network, this will result in an intermediate tissue form that is less robust. Consequently, the match between in vivo and in vitro in this case reflects the missing interactions and it is therefore very well determined. If, however, the in vitro conditions fail to activate any of the subsets of gene/protein interactions of the module that can exhibit some robustness, then the degree of match is low and certainly less than the one required to call this process biomimetic, irrespectively of any macroscopic similarity of the in vitro and in vivo conditions that, importantly so, do not provide any absolute criterion for comparison. In any case, the lower limit of the degree of match from the one required to produce a modular/robust tissue intermediate can be determined from the relative robustness of the mathematical model that represents the partly activated operating subset of the network in comparison with the model of fully activated network. Von Dassow and Odell ²⁷ have also observed that if the small-scale regulatory mechanisms, such as positive and negative feedback loops, that had been assembled to construct the modular network with the robust behavior are not strongly connected to each other with interaction among their member genes/proteins, the loosely connected network exhibited less robustness than the one in which part of the interactions is missing. This is expected according to the concept of modularity we have primarily used for the tissue forms whose structural stability depends on the internal interactions and not on interactions with surrounding codeveloping tissues that are not stronger than the internal ones. In the case of the network too, the strengthening of the internal interactions is what makes this network a module with robust behavior, that is, in other words all these interactions should be active for the network to produce at each moment the same mRNA or protein concentration.

As for the possibility to take advantage of a relevant subset of developmental mechanisms instead of attempting the much more difficult task of complete recapitulation of tissue development, ³ it becomes obvious that this depends on the modularity/robustness of the gene/protein network that implements this mechanism. If this network is not robust, as for example in the case that is only a subset of another mechanism that corresponds to a robust larger network, any attempt to activate this subset is less probable to succeed. Inventing in vitro conditions to do so instead of activating the full network will probably result in the loss of the robustness of the process.

It seems from the above that a systems biology approach to in vitro developmental process is the most suitable to quantify the process design and make it biologically meaningful at the tissue level and computational based on measurable/calculable variables (gene expression, protein concentration, their spatial distribution, etc.). Systems biology not only treats information about developmental phenomena computationally but also primarily determines which information is relevant in each case. So the problem of the process designer becomes more specified, instead of looking to any information, such as gene expression or in vivo conditions, for a particular developmental stage that he/she has to implement in vitro, he/she has only to consider this information that is related with the activation of a gene/protein module, either known in developmental biology or suspected according to the existing evidence.

Robustness can be used as a design criterion for construction of the correct model

To make practical use of systems biology in process design, the first critical question is the relative ease of constructing useful mathematical models from available information. This information may be available for some model systems well studied by developmental biologists, but analogous data for tissues are less complete. Nevertheless, the information needed to develop adequate gene network models will not be impractically large, because model robustness can be used as a simplifying criterion for the model construction. Von Dassow et al. did exactly this, modifying the initial model, which was not robust, by adding two interactions to achieve a model that was thus corresponding to biological reality. ²⁶ This method has also been followed by Eldar et al. ³³ to identify the interactions that give rise to another developmental module of Drosophila, the robust bone morphogenetic protein (BMP) morphogen gradient. The criterion of robustness limits the possible design solutions for networks that could define the pattern. Thus, different patterning mechanisms could be distinguished according to their robustness. To understand the patterning mechanism, the authors developed a general model based on available molecular knowledge. They next used computational methods to screen for robust networks giving the correct pattern, varying the parameters as needed to fit these criteria. Although most of the networks evaluated gave the correct pattern, only a small fraction, less than 1%, was robust to twofold changes in gene dosage (parameter variation). It is noteworthy that the robust networks displayed several unique properties. For example, for all the robust cases, the ligand was diffusible after binding to an inhibitor, and only a bound protein, but not a free ligand, was cleaved by a protease, in accordance with one of the previously proposed mechanisms.

In addition to helping define qualitative characteristics, the robustness criterion imposes limits on the possible range of parameters, many of which have not been determined experimentally. In the segment polarity model, nearly 50 free parameters were not known. ²⁶ Appropriate parameters for the model were found by random sampling over a plausible range of parameter values and imposing the robustness condition. A similar approach has been used to explain the ability of Xenopus embryos to maintain appropriate proportions regardless of size. ²⁴ Thus, this approach is likely to have general utility. The knowledge of the complete pathway was not needed. Pathways can be represented in a simplified form based on an input–output relation. Neither the Drosophila segment polarity model ²⁶ nor the Xenopus scaling model ²⁴ used the complete pathway between input/signal and target, but instead used a dose–response relationship. Thus, the complete pathway between wg and en can be represented simply in the segment polarity model as a term in the equations that describe the induction of the gene en by the protein of wg, ³³ and scale-free patterning of the Xenopus gastrula can be described in terms of the concentrations of just four proteins. ²⁴

The difference of in vivo and in vitro modeling

There is, however, an important difference between the in vivo and in vitro situations. Contrary to the situation in vivo, where the pattern forms by default, establishing the desired pattern in vitro requires a process design that includes explicitly described stages of the process and their implementation in bioreactor/biomaterial systems. In the example below, we will show how a biomimetic model based on information already available regarding the developmental biology of the growth plate can be enhanced with the help of computational systems biology to define practical in vitro conditions for the activation of the appropriate gene/protein interaction module.

A hypothetical, but testable biomimetic process for growth plate

In Part I, using information from developmental biology, we gave as an example a biomimetic process that has reasonable chances to lead to a structure close to that of the growth plate. This process is assembled from four subprocesses that could achieve a gradual progression of cell differentiation and construct size through the subprocesses that, as we explained in Part I, is important for the maintenance of cell-to-cell communication. This communication could be disrupted easily if the construct size becomes disproportionate to the cell signaling capabilities, since these are determined by the cell differentiation stage. Such a gradual progression of cell differentiation and construct size also leads to a correspondingly orderly organization of the cells into the growth plate pattern. Chondrocytes are distributed almost randomly in early long bone primordia. ³⁴ They are arranged later, as the tissue develops into groups or clones in the presumptive growth plates, eventually becoming stacked in columns. ³⁵ This gradual organization or patterning, with chondrocytes at different positions within a column expressing different genes corresponding to different differentiation stages, is reminiscent of Drosophila segmentation, described earlier, which becomes more precise at each successive developmental stage, ultimately giving rise to the 14 segments (Fig. 2).

The four subprocesses are as follows (Fig. 3):

OPEN IN VIEWER

FIG. 3.

A biomimetic process for in vitro formation of a bioartificial growth plate assembled by four subprocesses that could achieve a gradual progression of cell differentiation and construct size along these subprocesses. This is important for maintenance of cell-to-cell communication leading to correct cell organization. Subprocess 1: cell expansion on a microcarrier surface without induction of cell differentiation. Subprocess 2: small size clusters of microcarriers bearing cells undergo chondrogenic differentiation induced by the high cell density. Subprocess 3: this subprocess cultures the cells for longer periods than the subprocess 2, giving as output larger clusters on the microcarriers, in which the cells have differentiated further, but not to the final stage of mineralization. Subprocess 4: two hydrogel layers, each one filled with the constructs of the subprocesses 2 and 3, forming a primitive zonal (not columnar) structure statistically similar to that of the growth plate (prepatterning) which, under appropriate conditions, could lead to a robust pattern for growth plate.

Subprocess 1. Mesenchymal stem cells attached on microcarriers grow inside a rotation bioreactor covering the microcarrier surface. The first subprocess achieves cell growth on the microcarrier surface and corresponds to cell expansion, which increases the number of cells without inducing their differentiation.

Subprocess 2. Mesenchymal stem cells attached to microcarriers from the output of subprocess 1 are the input for subprocess 2. Small size clusters of microcarriers bearing cells start to be formed as the cells proliferate and cover the microcarrier surface in multiple layers. At the same time, chondrogenic differentiation is induced in areas of high cell density. Subsequently, cell aggregation in areas between microcarriers takes place because after the coverage of the surface of microcarriers by cells, the formation of larger multimicrocarrier clusters will occur. This happens through microcarrier-to-microcarrier connection from cells protruding from the microcarrier surface, which form bridges between the microcarriers that hold them together, a phenomenon observed previously. ³⁶

Subprocess 3. The output from subprocess 2 is the input for subprocess 3. This subprocess cultures the cells for longer time periods than subprocess 2, giving larger size cell-microcarrier clusters as output. In such clusters, the cells have differentiated farther, but still not to the final stage of mineralization.

From the subprocesses 2 and 3, we have two constructs of different size and in different differentiation states, neither of which corresponds to the final state of the bioartificial growth plate. At this point, no cell organization has started, except that the differentiating cells are located in the center of the clusters.

Subprocess 4. In this subprocess we use two hydrogel layers, filling each one with the constructs of the subprocesses 2 and 3, and place them on top of each other, with the hydrogel layer having the differentiating cells in an earlier state on the top. This construct has a primitive zonal (not columnar) structure statistically similar to the one of the growth plate. The first layer contains less chondrocytes in advanced differentiation state, for example, more proliferating than prehypertrophic or hypertrophic, whereas the other contains more cells in more advanced differentiation state, for example, more hypertrophic cells. We could say, therefore, that the bilayered input of this subprocess is a construct of large size (gradual increase of the tissue size as the subprocesses continue), with a statistically “organized” structure based on the gradient of cell differentiation state from the first to the second hydrogel layer. Multilayered hydrogel systems allow cell signaling along interfaces between cells in different layers (Elisseeff in Mikos et al. ³⁷ ), and therefore cell-to-cell signaling through the negative feedback loop of Ihh/PTHrP, which is suspected to be responsible for the robust structure of growth plate, can take place.

The bioartificial construct, which is the input of the subprocess 4, has not yet established intrinsic control over the rate of chondrocyte differentiation. At a minimum, the negative feedback loop of Ihh/PTHrP, which retards this rate and leads to a columnar structure, will be required not only to facilitate cell self-organization (see Part I) but also to start applying systems biology methods. The construct with the two hydrogel layers provides a primitive/statistical zonal structure, or prepatterning, which allows the directed exchange of signals mediated by the Ihh/PTHrP proteins. Ihh will be released mostly by the second layer, which contains more mature chondrocytes; PTHrP will come from the less mature chondrocytes in first layer. Although of course Ihh and PTHrP will also diffuse within the layers where they originate, their target cells reside predominantly in the other layer. Under this statistically monodirectional signal exchange, the new chondrocytes leaving the proliferating state to enter their differentiation program will be under the influence of this loop and will be aligned along the signal gradient from the first to the second hydrogel layer, thus producing a columnar organization for the newly differentiating cells. This also is reminiscent of the prepatterning of Drosophila segmental pattern by the pair-rule genes (see above, section 3).

The transformation from statistical prepatterning to the columnar pattern is in agreement with developmental biology. Schipani et al. have generated transgenic mice in which constitutive expression of PTH/PTHrP receptors was targeted to the growth plate through a collagen II promoter. ³⁸ As expected, these mice showed a delay of chondrocyte differentiation because PTHrP, secreted by the proliferating chondrocytes, acts on prehypertrophic cells to retard their progression to the hypertrophic stage. Moreover, the growth plate was disorganized. The tibial hypertrophic chondrocytes appeared at the periphery of the diaphysis and not in the center as normally seen. However, when these mice were mated with mice that did not express PTHrP, “rescue” was seen; chondrocyte differentiation was accelerated, and the animals with both cell types had correctly patterned growth plates. The important point for our discussion is that a mixture of cells with “opposite” abnormalities, one expressing the receptor constitutively and the other not expressing the ligand, randomly distributed in the developing growth plate, cooperated globally in the whole tissue space to balance their effects and restore normal patterning. This is an example that shows that the partial organization of growth plate that we have at the beginning of subprocess 4 can be made complete by the effect of newly formed cells, which will behave differently in respect to their orientation by virtue of the fact that they end up in a different environment that facilitates their organization.

In conclusion, in subprocess 4, a prepatterning has been achieved, leaving open the possibility of further development to the final columnar pattern of the growth plate. However, for this to take place, the gene/protein module that gives this pattern has to be activated and therefore the conditions, that is, the parameters and initial conditions of the subprocess (e.g., hydrogel porosity or thickness, the number of layers, and expression levels of Ihh/PTHrP) have to be found for its activation. Some of these can be estimated easily. For example, tangential flow is the likely appropriate mode for the bioreactor fluidics, because a flow through the construct in subprocess 4 will washout any secreted protein, in agreement with the in vivo situation, where vascular regression is required. ³⁹

There are obvious similarities between the segment polarity pattern and the pattern of growth plate. As the segmental pattern results from the interaction between two adjacent cells exchanging wg and hh proteins to activate the signaling pathways in each cell that form an integrated gene/protein expression network, the proteins Ihh and PTHrP are exchanged between chondrocytes at different differentiation stages to different zones or column heights (c), forming a negative feedback loop (Fig. 4).

OPEN IN VIEWER

FIG. 4.

As the segmental pattern of Drosophila results from the interaction between two adjacent cells that exchange the protein signals coded by the genes wg and hh, activating the signaling pathways in each cell that together integrate to form a gene/protein expression network (a) similarly in the growth plate (c), the proteins Indian Hedgehog (Ihh) and parathyroid hormone-related protein (PTHrP) are exchanged between chondrocytes at different differentiation states and consequently located at different zones (b) or column heights (c), forming a negative feedback loop that is responsible for the balance of proliferation and differentiation that leads to the robust columnar pattern.

In vivo studies similar to those that identified the key genes/proteins involved in determining segment polarity have also been performed for the growth plate. Research on gene expression patterns identified genes/proteins involved in the organization of differentiating chondrocytes in the form of parallel columns (e.g., Ref. ⁴⁰ ). It was found that the proteins controlling the balance between chondrocyte proliferation and differentiation set limits on the column elongation rate. These limits are needed for the size increase of the growth plate (cell proliferation, differentiation, and hypertrophy) to be properly synchronized with column elongation and other processes taking place in parallel, such as secretion and organization of the extracellular matrix, both of which are needed for columnar structural integrity. For example, a complete disappearance of the columnar architecture was observed in transgenic mice having a mutation in the type II collagen gene. ⁴¹

These molecular signals are exchanged between chondrocytes at different differentiation states and consequently located at different zones or column heights. To a first approximation, a negative feedback loop consisting of Ihh and PTHrP is responsible for the balance of proliferation/differentiation and columnar organization. The signals are exchanged between proliferating and prehypertrophic chondrocytes located in distinct zones.^15,42 As a consequence, cell distances depend on the number of chondrocytes in the column in one zone, but in the other, it is these distances that determine the chondrocyte differentiation rate, since this rate depends on the transport of protein signals along the column (Fig. 5). In response to Ihh signaling, PTHrP is secreted by the proliferating chondrocytes, diffuses to the prehypertrophic zone, and retards differentiation of the cells in this zone to the hypertrophic state, thereby delaying column elongation. This in turn leads to a time-dependent decrease of Ihh secretion because of the decrease in the number of cells that produce it. Consequently, a decrease of PTHrP results from the decrease of Ihh. This negative feedback loop thus sets the limits in rate of chondrocyte differentiation and consequently column elongation.

OPEN IN VIEWER

FIG. 5.

Interactions between proliferating and prehypertrophic chondrocytes in the growth plate. In response to Ihh signaling, PTHrP is secreted by the proliferating chondrocytes, reaches the prehypertrophic chondrocytes, and retards their progression to the hypertrophic state, thereby delaying column elongation and synchronizing it with cell differentiation and other processes taking place in parallel, such as secretion and organization of the extracellular matrix, which is needed for columnar structural integrity. In addition, Ihh regulates the expression of BMP genes, which also upregulate chondrocyte proliferation. BMP, bone morphogenetic protein; P, proliferation; D, differentiation.

Applying the robustness criterion, ³³ using computational methods as described earlier to determine the minimal set of interactions required for the gradient, we could include additional signals known to influence the balance of proliferation/differentiation. For example, Ihh increases the proliferation rate of chondrocytes ⁴³ (Fig. 5). In addition, Ihh regulates the expression of BMP genes, which in turn can also upregulate chondrocyte proliferation ⁴⁴ (Fig. 5). In both of these cases the proliferating chondrocytes start differentiating, moving along the column until they escape from the range of PTHrP signaling. Leaving the proliferation zone, these chondrocytes can produce Ihh. van der Eerden et al. present a more detailed scheme for the major signaling pathways operating among chondrocytes ⁴⁵ that could provide additional input for the model. The optimal set to be included in the model can be identified by optimizing robustness using computational methods as described earlier.^26,33

Additional data related to known physical and kinetic properties of the growth plate in vivo and additional aspects of chondrocyte biology can simplify the process by restricting the computational analysis to a reasonable physiological range for the parameter values of the model. Examples of such data include diffusion coefficients, ⁴⁶ solute transport rates, ⁴⁷ volume increase of hypertrophic chondrocytes (it could be 10–15%^48,49), cell proliferation rates, ⁴⁹ cell growth by cellular division, matrix synthesis throughout the growth plate, and chondrocytic enlargement during hypertrophy,^50,51 and other relevant information from the literature.

In addition, expression of many genes along the zones of the growth plate has already been determined accurately.^52,53 Thus, we have sufficient data from in vivo developmental biology and an experimentally validated approach to provide further constraints using the robustness criterion to allow construction of a provisional computational model. We can even use the software developed by Meir et al., based on the program used by von Dassow, which is publicly available. ⁵⁴

To construct and refine the model, we should be checking model parameters and initial values for robustness of the model solution: a stable columnar pattern (or in the simplified version just one column). Therefore, different versions of the model will needed to be screened, adding or subtracting input genes/proteins, or their interactions and scanning plausible values of initial conditions and parameters until the simplicity and robustness of the model are optimized. In this case, we can start with a gene/protein interaction module that is known to be responsible for the columnar pattern. Since our interest is not the module itself but how we could use its model to design the in vitro subprocess 4 of Figure 3 to make it correspondingly robust, we proceed to the step of modifying the model parameters (e.g., physical parameters of hydrogels—instead of using diffusion coefficients of the in vivo growth plate we use those of hydrogel, proliferation rates) and initial values that will account for the hydrogel layer structure as provided by subprocesses 2 and 3 (initial number of cells, different cell differentiation states in the same zone, etc.). Again, we apply the criterion of robustness for parameters and initial conditions. When these are found, we return to subprocesses 2 and 3 and make the necessary modifications to implement them in the two hydrogel layers. In this way, we arrive at an organized, robust mathematical model for a biomimetic growth plate under conditions that assure corresponding robustness of its physical structure.

Although the above approach may seem overly complex, we should consider that the modular intermediate tissue forms are of particular importance in tissue engineering because they are the building blocks of complex tissues. In the case of the growth plate, after establishment of the columnar pattern we could combine subprocess 4 with other subprocesses to generate osteoblasts, and depending on the details of our manipulations, get either bioartificial osteochondral tissue or bioartificial bone, recapitulating endochondral ossification (see Part I for initial attempts at in vitro endochondral ossification) or we could just use the bioartificial cartilage alone. By including osteoclasts generated by another subprocess, we could restore the function of remodeling, to test for bone resorption before the bone is implanted. ⁵⁵ In this case, we would have to deal with interactions between osteoblasts that control the degree of osteoclastic activity and osteoclasts, which control osteoblasts differentially depending on their stage of differentiation. Sufficient information exists for the molecular details of these interactions that they have already been integrated in a mathematical model, ⁵⁶ though not from the point of view of systems biology. The determination of such a module and appropriate in vitro conditions may be of importance for bone engineering. This will need to be verified in vivo as this balance may be different in vivo. However, having a mathematical model of the module allows its robustness tested in terms of in vivo conditions.

In conclusion, there are several tissues for which sufficient information is available to allow design of biomimetic process using computational methods and identify optimal conditions using the criterion of model robustness. That will become apparent as long as tissue engineers familiarize themselves with the concepts of developmental biology, mentioned in Part I, and principles of systems biology relevant to process design.

The most important concepts are those of modularity and robustness. In the next sections we will encounter the same concepts again in cases where several modules can be activated from different initial conditions, as in the case of stem-cell differentiation. Although the methodology of systems biology cannot yet be applied in this case because of the lack of data, we will show how a different formalism, that is, the concept of attractors, makes it possible to design experiments and interpret the results in a way that allows the evaluation of modularity and robustness (stability).

Experimental evidence of trajectories converging to attractor states

While the adoption of stable network topologies that provide stable gene expression as gene interactions take place was shown in computer simulations several years ago,^29,30 recent experimental work has confirmed that the concept of attractors represents essential aspects of cell differentiation.

Huang et al. followed two different trajectories during the differentiation of human promyelocytic HL60 cells to neutrophils, triggered by two different factors, dimethylsulfoxide (DMSO) and all-trans-retinoic acid (atRA). The authors observed that the initially divergent trajectories finally started to converge, as estimated by gene expression patterns. ⁶² This could be represented in Figure 7 (representing graphically the changes of only two of the protein or mRNA concentrations along the trajectory), as a disturbance of the initial conditions, protein/mRNA concentrations, under the influence of the factors I (DMSO) or II (atRA), that induce gene expression changes. Either of these perturbations changes the initial state (state “a,” Fig. 7) to a state falling within the basin of attraction (“b” or “c,” Fig. 7) of the attractor corresponding to the mature neutrophils. The authors followed the expression of around 12,600 genes. After excluding the genes whose expression was too low or did not exhibit any significant change during the experiment, they selected 3841. Initially (12–18 h after treatment), the mRNA expression profile indicated that the two trajectories diverged. This is not surprising, because the different factors used targeted different genes. After the initial divergence, however, the trajectories converged to similar patterns by day 6, with 72% of the initial set of 3841 genes having identical expression levels. The authors considered that the final gene expression was sufficiently similar to consider that the initially diverging trajectories converged to a common state. They attributed the differences in expression of the rest of the genes to the possibility that some of the genes induced by DMSO and atRA may not be relevant to the macroscopic definition of neutrophils. There are probably other reasons why a “state” cannot be defined accurately in biological systems simply by a list of expressed genes. One property of biological systems is degeneracy. Contrary to redundancy, which requires an identical function to be performed by different elements, in degenerate systems structurally different elements give the same or different functions depending on the context of their expression. ⁶³ Such degenerate systems are the transcription factors, where different factors can generate similar patterns of gene expression. For example, Fambrough et al. studied the relationship between receptor tyrosine kinase–activated signaling pathways and the transcriptional induction of immediate early genes. Even mutant receptors lacking binding sites for activation of the PLCgamma, PI3K, SHP2, and RasGAP pathways still retain partial ability to induce 64 of the 66 fibroblast immediate early genes examined. ⁶⁴ It may happen therefore that either the macroscopic definition of neutrophils is not accurate and different states can be present, as in the example of Chang et al., ⁶⁵ or there is only one state but it can be realized with a different set of genes, as in the examples of Wieghaus et al. ⁶⁶ and Zhong et al., ⁶⁷ who found signaling pathways activated without change in the expression of major genes known to be related to these pathways, in accordance with the degeneracy concept (mentioned below in section 6).

OPEN IN VIEWER

FIG. 7.

An attractor state represents a stable topology of the network of gene/protein interactions (here a geometric simplification in two dimensions, e.g., two interacting genes). Surrounding the attractor state is its basin of attraction. If changes in the expression of genes such as those induced by factors I or II bring the initial protein/mRNA concentrations inside the basin of attraction, an autonomous sequence of further changes in the concentrations occurs, leading to the final concentrations represented by the attractor point. The sequence of the changes is represented by a trajectory that approaches (is “attracted” by) the attractor state in which the gene interactions have been stabilized, so that no further changes in gene expression and protein concentration take place. The larger the basin of attraction, the more robust the model and its solution (i.e., a spatial pattern).

In either case, however, the example mentioned above ⁶² clearly shows that the cell states are stable, since they attract different initial conditions, even though these states cannot be defined yet, as can be also seen from the lack of a formal definition of the “cell type.” Consequently, characterization of cell therapeutic products, which is required for regulatory purposes, cannot be achieved so easily. The examples described above indicate clearly that a list of expressed genes or activated pathways is not sufficient to characterize the cell type. Instead, it will be necessary to refer to more global information describing the structure of the gene network to determine whether basic cell differentiation functions have been restored, irrespective of the differences in the expression of particular genes. One level higher than the individual pathways, we encounter the gene/protein modular interaction networks, whose stability may provide better criteria for characterization of the cell state. Wieghaus et al. ⁶⁶ and Zhong et al. ⁶⁷ referred to the modular organization of the gene interactions to characterize the cell states, though the module definition used was not functional, but structural, based on the concept of modularity used in network science: internal connections denser than connections with genes outside the module, which implies robustness indirectly. However, even without the knowledge of the network structure, cell state characterization can be approached by determining the attractor states experimentally and estimating their stability to perturbations computationally, as illustrated at the end of this section.

Experimental evidence that trajectories move through sequential attractors

“Fate determination” is a term from developmental biology and refers to a differentiation process along a particular pathway. However, the transition of the cell state toward the final differentiation state is not continuous along the pathway, since of cell states are not continuous, but discrete. This was recognized over 50 years ago by Waddington,^58,59 who indicated that cell states “switch between distinct, well recognizable types” with intermediates rare and unstable. Each time the network of gene/protein interactions leads the cell state to that of an attractor, the cell adopts a state, with gene expressions or protein concentrations that bring it close to the basin of attraction of next attractor and makes it receptive to the influence of subsequent attractors. The successive states in vivo are reached with additional external interactions from neighboring cells or systemic factors that place the initial point of the trajectory inside the basin of attraction of the next attractor state. In vitro, externally added growth/differentiation factors or other manipulations could be used to bring the initial conditions inside the basin of attraction of each successive attractor. In Figure 8 the initial condition is “pushed” from attractor “a” (arrow a) into the basin of attraction of attractor “b” and converges to the corresponding attractor state. Next, it is pushed (arrow b) out of the basin of attractor “b” and into that of attractor “c.” These transitions are equivalent to the motion of a sphere through hills and valleys, with the bottom of the valleys representing the attractor states and the hills the borders of the basin of attraction (Fig. 8b, we have not included the “hills” unstable states in Fig. 8a). For the sphere to escape from a valley, an external force is needed; it then falls autonomously into the next valley. It is striking that the vocabulary of attractor basins, and the physical analogy of “hills and valleys,” was anticipated as early as 1957 by Waddington's concept of the “epigenetic landscape,” ⁶⁸ in which he presented the metaphor of a ball traveling down a landscape of branching valleys.

OPEN IN VIEWER

FIG. 8.

The path dependence of development can be perceived as the sequential convergence to successive attractors. The successive states are reached through changes (through manipulations during manufacture or natural developmental programs) in the initial conditions (e.g., protein or mRNA concentrations) that place them inside the basin of attraction of the next attractor state. Starting from attractor “a,” the initial conditions (e.g., protein concentrations) are “pushed” by factor A into the basin of attractor “b” and converge toward attractor state b. Subsequently, factor B initiates changes that place the system state within attractor “c” (a). These transitions are equivalent to the motion of a sphere through hills and valleys, with the bottom of the valleys representing the attractor states and the hills the boundaries of the basins of attraction (b).

Such a path-dependent transition between attractors has been shown experimentally by Chang et al., who concluded that cell differentiation is a “discontinuous switching between cellular states.” ⁶⁹ The authors used human promyelocytic HL60 cells differentiated to neutrophils with DMSO, monitoring the expression of CD11b, a well-established marker for mature neutrophils. During cell differentiation, CD11b expression appeared to increase gradually when analyzed at the population level. In contrast, CD11b expression at the individual cell level was either “all on” or “all off.” Cell-by-cell analysis revealed two distinct populations, one expressing low and the other expressing high level of CD11b. An additional round of DMSO stimulation converted low CD11b-expressing cells to CD11b high-expressing cells, showing that in the first round they had only partially proceeded to the final differentiation state. The authors attributed this behavior to the existence of two different cell states in the initial population that correspond to two different attractors, such as the attractors “a” and “b” in Figure 8b. The CD11b low-expressing subpopulation, under the influence of DMSO (arrow “a” in Fig. 8b), moves inside the basin of attraction of the next attractor (attractor “b” in Fig. 8b), the state corresponding to the CD11b high-expressing subpopulation. Further treatment with DMSO (arrow “b”) leads to the next attractor (attractor “c” in Fig. 8b), which is the final differentiation state. The CD11b high-expressing subpopulation has already progressed to attractor “b”; therefore a single treatment with DMSO can move it to attractor “c.”

One interpretation of these observations is that cell differentiation is indeed a multistep process realized through a sequence of discrete intermediate attractor states. These studies open the door to practical use of the attractor concept in tissue engineering. Transitions similar to those described above for myeloid cell differentiation take place in various developing tissues.^70,71

While these studies refer to state spaces describing the behavior of individual cells, they could be elaborated further by including attractors composed of heterogeneous cells distributed in space. In such a model, global, spatially distributed gene networks are composed of several different interacting cellular gene subnetworks. One example is the developmental module of the growth plate described earlier, whose cellular interactions are determined primarily by the Ihh/PTHrP feedback loop, which consequently determines the network of genes expressed in each stage of chondrocyte differentiation. Recently, in addition to the Drosophila embryo,^23,26 the behavior of the amphibian gastrula has been described accurately with a quantitative model based on a remarkably small number of parameters, and predictions made by the model were confirmed experimentally. ²⁴ Experimental evidence that such “parsimonious” models can predict the behavior of a complex vertebrate system accurately indicates that using similar approaches to characterize the behavior of engineered tissues in vivo will be feasible. Other efforts to connect the gene with the tissue level are in progress.^72,73

A manufacturing process can thus be seen as a sequential process in which the system is moved through a mathematically defined “state space” from one attractor to the next using defined stimuli, based on models that accurately predict the movement of the cell states through such a space.

The attractor map as a template for process design

Based on the concepts we have presented, we propose that differentiated cells would be manufactured via a series of subprocesses, each corresponding to a sequential attractor state. Since many sets of initial conditions will fall within a given attractor basin, multiple design solutions are possible (as in the case of the robustness of the segmental pattern to initial concentrations of proteins or mRNAs, i.e., to different prepatterning), and process controls can be conceived with the primary objective of placing the system state within the basin of attraction in a statistically robust manner. Moreover, parsimonious models may suffice to describe most attractors. From a bioengineering perspective, this would result in a substantial understanding and simplification of the overall design problem. Attempting to achieve a similar result with piecemeal, trial-and-error addition of individual growth factors implicated in differentiation of a tissue, organ, or specific pathway is likely to be difficult at best. In contrast, if several modules are operating sequentially and/or in parallel during differentiation, the primary focus would be to identify conditions sufficient to activate each module reliably. Modularity will result in autonomous coordination of the various signaling pathways along a differentiation trajectory, without external intervention. Moreover, process validation could rely to a substantial degree on confirmation that each in vitro stage, or subprocess, places the cell state inside the boundaries of the corresponding attractor to a degree determined by well-established quality systems methods and in complete accord with the concept of unit operations described previously.

The concept of attractors was used recently to explain erythroid or myeloid differentiation of a clonal population of mouse hematopoietic progenitor cells, showing the usefulness of an attractor differentiation map for process design. ⁶⁵ The authors showed that the clonal population fluctuates between metastable attractors, one with a low level of the stem-cell surface marker Sca-1, one with a medium level, and one with high level. Metastable attractors are those with weak attractive power; therefore, cells can escape their attraction easily. In other words, their basin of attraction is small, indicating that the cell state needs only a slight disturbance in gene expression or protein concentration to escape it and transit to another more stable state. Cells with low Sca-1 expressed higher levels of GATA1, which favors erythroid differentiation; those with high Sca-1 expressed high levels of PU.1, which favors the myeloid pathway. The fact that in a clonal population there are cells with different gene expression profiles is in accordance with the existence of different attractor states. Persistence and specificity could easily be explained by stochastic noise on the order of the potential difference between adjacent attractor basins corresponding to bistable states—this is exactly how rod/cone fate is decided in insect compound eyes, a very well understood system. ⁷⁴ When low or high Sca-1–expressing cells were isolated and cultured, they reconstituted the parental heterogeneity of the three different states. This means that the cells can transit between metastable states, since these states reappeared from any of the isolated populations. This experiment shows that useful information for the characterization of cell state can be provided by experiments designed under the conceptual framework of attractors.

The practical significance of having metastable states of the system in a differentiation map is to guide strategies for design of processes to move cells to a desired differentiation state that is highly stable (i.e., corresponding to a wide basin of attraction). In the above example, if erythroid differentiation is desired, cells expressing low Sca-1 should be selected. Otherwise, all three populations will be present in the product, added growth factors will not have the same effect, and undesired heterogeneity will be created.

Another interesting observation is that when the cells are committed to either phenotype, it is not only the upregulation of lineage-specific genes that takes place but also the extinction of lineage markers characteristic of the starting cell population. In other words, when a new module is formed, another one is disassembled.^75,76 As we will see in section 6, network theory predicts that the disintegration of a topology characteristic of biological networks should take place in a specific way dictated by that topology. It would be interesting for this to be confirmed experimentally in the above and other cases. Such a finding would lead to strategies that use network science to define external interferences attacking the components of specific modules, leading to their disintegration. Thus, not only will modules be switched on during differentiation toward specific lineages but also those of other lineages will be switched off in a methodical way directed to very specific genes of the module. In the example of the study of Xia et al., two modules, proliferation and differentiation, were identified and characterized well enough to assure switching one off and the other on as desired. This may also be of relevance to diseases, as in these conditions there may be activation of a module, which under physiological conditions is inactive, and thus provide innovative treatment targets. In the case of p53 protein, a well-known tumor suppressor protein, the three scientists who discovered its action published an article, after thousands of papers, devoted to this molecule, in which they claim that instead of focusing specifically on the molecule, we should instead focus on the network of interconnections with other molecules, comparing this network with the Internet. ⁷⁷

In this section, we examined the problem of having several stable states, instead of a single one, into which the gene/protein network can settle during cell differentiation. Current limitations in data concerning the role of genes/proteins in cell differentiation impose corresponding limits on the development of models for small, modular subsets of the entire network operating during differentiation that might lead to modular and robust cell states defining any desired cell type. However, estimation of the cell state robustness can be done using the concept of attractors, as we pointed out earlier. We have seen that a representation of the transitions of the cell state during differentiation under the framework of attractors can provide information about the robustness of the differentiation process and cell differentiation state that could even allow experimental comparison between the robustness of different states. In sections 3 and 4, we used this as a basic criterion for process design.

In the next section, we will see how network science can provide additional information about the topology of gene/protein networks operating in cells. These studies may also help provide a clear definition of cell state. As has been mentioned, certain cell states can be described by degenerate sets of activated genes and/or proteins. It is therefore reasonable to examine global features of the network that may provide structural information sufficient to define a specific cell state uniquely. This in turn might facilitate the detection of modules, so that further, more elaborate studies with systems biology can be undertaken, and provide more general (in several distinct networks) and more global (genome/proteome/transcriptome scale) information about the network topology. This more comprehensive approach might then open the door to the discovery of general principles that could become a solid theoretical basis for tissue engineering.

Process design through modularity/robustness

In the previous sections, we encountered two types of situations. The first scenario corresponds to processes for which sufficient data exist from in vivo developmental biology, such as the case of segment polarity pattern of Drosophila, for the construction of a mathematical model that includes a limited set of gene/protein interactions found experimentally to play a major role. In the other cases, however, as for many cell differentiation pathways, sufficient information is not available. Even in the latter case, viable approaches exist. We have already seen (section 5) that the robustness criterion can be used to fill in many of the gaps for process design, with limited experimentation as we will now explore how we can use descriptions of the gene interaction network topology based on the concepts/methods of network science to overcome limitations in the available experimental data.

In the case of growth plate, we designed the process in a series of subprocesses, trying to recapitulate the different developmental stages. Here experimental data provide evidence that the growth plate is a robust developmental structure or macroscopic modular tissue form; accordingly, the design of a biomimetic process should aim to form appropriate developmental modules precisely because of their robustness. Therefore, we addressed the design problem using methods of systems biology already shown by experimental work to explain the behavior of well-understood macroscopic developmental modules occurring during development in vivo (section 3). In this case, the concepts and methods of systems biology were employed to connect the macroscopic robustness or modularity of a developmental pattern (spatially differential gene expressions) to the gene/protein interactions inside the cells, as well as between cells, forming a network extended in the pattern. The prepatterning—the initial spatial distribution of concentrations of proteins—activates this network. As the prepatterned initial proteins switch various genes on or off, protein concentrations or the gene expression levels change with time. This in turn causes further changes until the network converges to its attractor, and the module therefore stabilizes. The change of the prepattern (initial conditions) toward the final state can be represented as a trajectory, a sequence that approaches (or is attracted by) the final, robust attractor. With the use of a mathematical model we can try to find the gene/protein interaction network making use of the concept of robustness. As the macroscopic pattern is robust to environmental noise, its underlying network should also be robust to perturbations in parameters, as reflected in the mathematical model. In other words, the model should give the same solution pattern for a wide range of parameter values and initial conditions, corresponding to robustness observed experimentally in vivo. In section 4, we applied the same concepts/methods of systems biology for the growth plate pattern. To do so, we had to design the process according to the rules of developmental engineering provided in Table 1 in a series of subprocesses that recapitulate the stages of in vivo development, leading to a process in which a prepattern is imposed artificially (subprocess 4 in section 4). A mathematical model can then be developed in the same way as for the in vivo pattern to describe the macroscopic developmental pattern in terms of gene/protein interaction networks. This in turn allows identification of in vitro conditions under which the artificial prepattern will converge to the final robust pattern of the growth plate. The criterion of robustness was used to find the in vitro conditions (range of parameters and initial condition values), as well as the structure of the model itself (i.e., which interactions should be included). Implementing these conditions in the process, we can achieve organization of chondrocytes into the robust pattern of the growth plate. Because of the information available from developmental biology, according to which a macroscopic developmental module can be attributed to a small number of gene and protein interactions, the molecular network responsible for the pattern can be modeled either with ordinary differential equations ²⁶ or with abstract Boolean networks ³² (section 3 for the models of segment polarity pattern of Drosophila).

In the other cases, with little information about the gene/proteins interacting and limited evidence of modularity at the level of macroscopic cell state, as is the case for stem-cell differentiation, we could not apply the same rigorous mathematical analysis. Mathematical models cannot be constructed to address only a small subnetwork of the whole network operating in the cell responsible for the robustness of the cell states during differentiation, because this subnetwork is not known. Instead the whole network or a large subnetwork should be considered. This makes it impractical to use computational approaches for determining the module from the dynamics of the system (continuous update of the gene/protein concentration through interactions until a robust solution is reached).

However, the same concept of robustness, expressed differently in terms of stable states (attractors) and initial conditions leading to them (basin of attraction), can be used to design a process for cell differentiation in vitro even in the absence of a mathematical model. As discussed in section 5, we can generate information from experiments designed to evaluate the existence of attractors (different initial conditions lead to the same attractor cell state; Fig. 7, DMSO or atRA treatment of human promyelocytic HL60 cells ⁶² ), their position in terms of the disturbance the cell state needs to adopt a new state (switching between attractor cell states; Fig. 8, differentiation of human promyelocytic HL60 to neutrophils in sequential steps), ⁶⁹ or their relative stability/robustness (metastable states of low stability between the erythroid and myeloid high stability states ⁶⁵ ).

The other alternative is to make use of the interconnections of genes in a large set around the genes of interest, looking at structural features only, not the dynamics that describe how the network behaves in response to perturbations in parameter or initial condition values. This approach to the extraction of useful information ranges from descriptive to computational. Starting from a set of genes of interest suspected to be involved in the process under consideration, identification of important interactions with others in the network can be facilitated by information accumulated in the literature for the interaction among genes, which is directly accessible in databases. Recent work in biology as well as biological phenomena relevant to tissue engineering make use of such methods to examine the concerted regulation of genes in an extended network of gene interactions that can be provided from databases. For example, network analysis tools were used recently to study genetic pathway regulation by phthalimide neovascular factor 1 (PNF1), a small molecule that induces angiogenesis, with the objective of applying this factor to promote neovascularization in bioartificial tissues. ⁶⁶ The authors used the Ingenuity Pathway Analysis software coupled to the Ingenuity Pathway Knowledge Base (IPKB; Ingenuity Systems, Redwood City, CA), which covers more than 23,900 mammalian genes and includes millions of pathway interactions from the literature. The dataset generated from PNF1 stimulation of human microvascular endothelial cells was used to select genes eligible to generate networks. These networks were extended to the full genome, making use of the information about the gene connections included in IPKB software. The authors argued that although examination of single-gene regulation was useful, the network approach gave more valuable information for two reasons: first, the differential expression of some genes is very small, and second, coordinated perturbation of several signaling pathways could be elucidated. Indeed, they observed several interconnected networks centered around the transforming growth factor-beta1 (TGF-β1) signaling pathway, which has many known effects on angiogenesis, and identified angiogenesis-related cellular processes that were activated. Strikingly, while PNF1 activated the network around TGF-β1, TGF-β1 differential expression was negligible.

To identify genes that play a critical role in pluripotency, Zhong et al. applied network analysis to gene expression data to compare the profile of pluripotent human embryonic stem cells (hESCs) and hESC-derived astrocytes. ⁶⁷ The authors found that the expression of multiple members of the p53 pathway was lower in hESC-derived astrocytes than in hESCs. However, p53 itself was not expressed differentially in their experiments. This finding illustrates a compelling advantage of looking at a more global scale for the effects of any intervention in the cell state instead of restricting the analysis to single genes or pathways.

From systems biology to networks: Looking for robustness

In making full use of the fact that the genes and proteins in a cell are interconnected in a huge network whose topology may reveal information suitable for process design, we should refer to the network science. Network science is an emerging scientific discipline that examines the topological characteristics of diverse types of networks, such as physical, informational, biological, or social networks, seeking to discover common principles of their architecture or topology. Network science provided evidence that some topological features are shared universally by networks of very different origin, suggesting that universal network design rules exist. ⁷⁸ The most important discovery of network science is that diverse networks, such as biological (gene/protein interactions, metabolism), technological (e.g., linked web pages in the internet), and social (e.g., scientists linked by coauthorship, or actors linked by playing in the same film), have an architecture that can be described with a few simple common design principles. The most remarkable characteristic is that these networks and the nodes (points of connection between network constituents) do not follow a random distribution in terms of links they have. In random networks there is a peak in distribution, meaning that the majority of the nodes have the same number of links, with very few nodes having much lower or higher links. Such networks have many nodes with the average number of links, that is, they have a scale. In contrast, most of the real networks examined to date do not follow a random distribution. In real networks, most nodes have very few links; there are a few nodes with a very high number of links; these are called “hubs” ⁷⁹ (Fig. 9). This means that there is no representative of the network node with an average number of links. Therefore, the network has no scale (“scale-free” networks; note that the network describing how Xenopus embryos maintain correct anatomical proportion regardless of size ²⁴ provides a direct physical correlate of this concept). In Figure 10a, we see a network with two hubs (black nodes), but no hubs exist in the network of Figure 10b. The nodes with few links, though many, cannot by themselves ensure that the network is fully connected, that is, that someone can “pass” through the network from node to node following the links as from node A to node B (Fig. 10); at some nodes there will be no link to proceed to the other part of the network. It is the hubs that can hold together the numerous nodes with few links and assure that the network is fully connected. Random networks can be easily broken into unconnected subnetworks by removal of just a few nodes (Fig. 9, transition from c to d). However, scale-free networks, such as those describing biological systems, are robust thanks to the hubs, and therefore do not disintegrate into unconnected small networks upon random removal of nodes, since most of the nodes are connected to very few others (e.g., transmission from a to b in Fig. 10). This explains why many mutations have little or no effect on phenotype. ⁸⁰ However, networks are extremely vulnerable to the removal of hubs^81,82 (transmission from a to c in Fig. 10). This is raised to a quantitative criterion for process design to retain the robustness of a developing intermediate tissue form. The elucidation of the topological features of the modular networks and the determination of their “error and attack tolerance” ⁸² (or simply robustness) can determine in vitro conditions that protect important nodes such as the hubs, which are critical to the integrity of the module. Alternatively, we may have to disintegrate an undesired module to restore or activate the desired alternative.

OPEN IN VIEWER

FIG. 9.

In scale-free networks (a), as the biological ones, most of the nodes have few links and there are few nodes with a high number of links (“hubs,” black nodes). In random networks there is a peak in distribution, meaning that the majority of the nodes have the same number of links and there are no hubs (b). It is the hubs that can hold together the nodes with few links in scale-free networks and assure that the network is fully connected. Random networks can be broken easily into unconnected subnetworks by removal of only a few links [(c) connected and (d) broken].

OPEN IN VIEWER

FIG. 10.

Hubs make scale-free networks robust. These networks do not disintegrate to unconnected small networks when random nodes are removed since most of the nodes are connected with very few others (transition from a to b). However, they are extremely vulnerable to the removal of hubs (transition from a to c) A and B are nodes.

Thus, we have again encountered the concept of robustness. In this case, it is not dynamic, as it was for the segment polarity pattern and its corresponding gene/protein network (section 3), but static, referring only to the topology of the network. Although the topological definition of these two concepts has limitations, it has been shown that processes like metabolic fluxes taking place in scale-free networks such as the metabolic reactions also exhibit scale-free topology. In the case of metabolism, the flux distribution is scale free with reactions having fluxes that span orders of magnitude, that is, most reactions have very small fluxes and coexist with a small number of reactions with very high fluxes. ⁸³

Modularity is another topological feature of diverse networks. The existence of highly interconnected groups of nodes that form modules has been observed not only in biological networks (e.g., gene, protein, or metabolic networks) but also in social (circles of friends) or technological (websites or discussion groups related to similar topics) networks. ⁸¹ The following example shows the utility of the modularity concept as defined in network science for in vitro process design. Xia et al. analyzed protein–protein interaction networks by transcription profiling and found two major modules, corresponding to proliferation or differentiation, that were negatively correlated. ⁶⁰ They have also detected proteins at the interface of the two modules, such as histone deacetylase and serum response factor, that according to literature data control the switch between the modules, that is, which will be activated and inhibited. Wang et al. have examined the network around Nanog protein, which acts in concert with other proteins such as Oct 4 and Sox2 to maintain the pluripotency of ESCs. ⁶¹ A protein interaction network was constructed composed of many proteins whose genes are putative direct transcriptional targets of the other proteins in the network. Because of the tight connection between the proteins of the network, the authors considered this to be a module responsible for cell pluripotency. Further, they raised the question of the vulnerability of the module to downregulation or inactivation of any one of its many components because of their high connectivity. Removal of a critical hub in this network could lead to disintegration of the modular topology and rapid loss of pluripotency. In a recent study, the same group analyzed the transcriptional network of pluripotency and detected the hubs—genes whose promoters are targets for a higher number of transcription factors and therefore have connections to many other genes. ⁸⁴ They have also found that there is a striking correlation between the number of bound factors and the possibility that a target gene is expressed in undifferentiated ESCs and then repressed during differentiation.

Another topological feature of networks is that modules are combined to form modules of higher order and these in turn are combined again for even higher order modules, leading finally to a nested organization of the network.^85,86 Returning to the question of the possibility of exploiting a relevant subset of developmental mechanisms instead of attempting the much more difficult task of complete recapitulation of tissue development, ³ we are in a position through network science to give a more accurate answer. The lack of detailed knowledge about the nested hierarchical structure of the cell gene network modular organization^27,87 and limited knowledge of the correspondence between modules and cell functions make highly specific activation of selected modules a difficult task at present. Compounding this, a theoretical danger exists that any deviation from the in vivo process could result in failure to activate a required developmental module, causing loss of the robustness that has been incorporated by evolution into the existing modular organization. Therefore, the approach of selective module activation should certainly be pursued, accompanied by further network studies for elucidation of modular organization, so that in the near future the external interventions in the process can be designed rationally from these concepts. Such efforts are already underway. In the future, artificial regulatory circuits and even modules with preselected inputs and outputs could be synthesized, constructing artificial developmental pathways that retain only the minimal essential characteristics for a particular process according to the spectrum of properties of bioartificial tissues needed in defined applications. ⁸⁸ Early efforts have demonstrated the ability to rewire signaling pathways that normally activate a module known to induce cell growth into a module that causes cell death. ⁸⁹ However, the long series of developmental events leading to the formation of a fully mature tissue and the corresponding activation/deactivation of modules during the process necessitates a detailed understanding of all the modules involved and their order of activation, information that at present is incomplete.

We have shown that the concepts of modularity and robustness can be applied to various types of processes with different demands and restrictions in terms of design. As in the case of the segment polarity network, where the robustness can be determined in a dynamic mathematical model, in stem-cell differentiation one could determine robustness from the topological analysis of the gene network without developing a systems biology mathematical model but instead by statistical analysis of the hubs and their distribution and specific location in the overall gene network. The importance of the concepts/methods and measures of network science for in vitro process design relies primarily on the universal character of the topological features of diverse networks, biological or not. These general features of networks might provide powerful principles for biological organization in general and therefore a sound scientific basis for tissue engineering.

Below we will refer to an example of biomimetic process design using the concepts/methods of network science to see how they could be integrated into process design computationally.