Structural requirements and discrimination of cell differentiation networks
Abstract
Mathematical models of stem cell differentiation are commonly based upon the concept of subsequent cell fate decisions, each controlled by a gene regulatory network. These networks exhibit a multistable behavior and cause the system to switch between qualitatively distinct stable steady states. However, the network structure of such a switching module is often uncertain, and there is lack of knowledge about the exact reaction kinetics. In this paper, we therefore perform an elementary study of small networks consisting of three interacting transcriptional regulators responsible for cell differentiation: We investigate which network structures can reproduce a certain multistable behavior, and how robustly this behavior is realized by each network. In order to approach these questions, we use a modeling framework which only uses qualitative information about the network, yet allows model discrimination as well as to evaluate the robustness of the desired multistability properties. We reveal structural network properties which are necessary and sufficient to realize distinct steady state patterns required for cell differentiation. Our results also show that structural and robustness properties of the networks are related to each other.
1 Introduction
Stem cells and their potential to give rise to multiple cell types have become a focus in systems biology and mathematical modeling throughout the last years (Peltier and Schaffer, 2010; MacArthur et al., 2009). Starting from a multipotent stem cell state, they undergo the process of cell differentiation which is completed when the cell ends up in a mature cell state. The eventual goal in stem cell research is to guide the differentiation of stem cells toward a desired cell type safely, and to maintain this state robustly. Therefore, a promising yet challenging task is to reveal the processes driving cell differentiation by mathematical modeling.
Cell differentiation is commonly viewed and modeled as a sequence of cell fate decisions (Foster et al., 2009). Thereby, the individual decision steps can be represented by modules, each driven by a small gene regulatory network of “master” transcriptional regulators (TRs) in which feedback motifs play a key role (Xiong and Ferrell, 2003; Huang et al., 2007; Foster et al., 2009). The development of models for these modules however is often severely hampered by the lack of detailed knowledge about molecular processes, shifting the focus toward qualitative models in order to understand the fundamental mechanisms. A variety of publications has contributed to the modeling of (stem) cell differentiation, covering a broad range of diverse modeling approaches, e.g. Roeder and Glauche (2006); Huang et al. (2007); Narula et al. (2010); MacArthur et al. (2008). Usually, specific stem cell systems such as hematopoietic or mesenchymal stem cells are considered. Yet, similar motifs and system properties can be found among them. For example, all these approaches have in common that they represent distinct cell types by distinct stable steady states, declaring multistability as a universal system property.
One of several aspects within this field concerns the question of stable cell states, and minimalistic models that are capable of reproducing experimentally observed cell states. Although it has been established to view cell differentiation as subsequent decisions driven by multistable switching modules (Foster et al., 2009), the network structure of the single modules is often still unclear. We therefore perform an elementary study of small networks consisting of three TRs. These networks should be able to exhibit one progenitor state and two competing differentiated cell types. We also compare two different hypotheses about the steady state levels of the TRs in the individual cell types. In order to reveal properties of 3node networks that provide candidate models for the differentiation process, we address the following questions:

Given three cell typespecific TRs, which networks of interactions between the TRs can reproduce the required stable steady states?

Among the networks that fulfill the stable steady state requirements, are there differences in robustness of the required multistability properties with respect to perturbations in the interactions?
As these are very basic questions, this work builds upon a qualitative modeling framework which needs only very few assumptions about the reaction kinetics. In order to model the influence of a TR on its target, monotonous activation and inhibition functions are used, but no exact kinetic parameters have to be specified. This framework has been previously used for example in Chaves et al. (2008); Breindl et al. (2010).
The outline of this paper is as follows: In Section 2, we describe the modeling framework, and define several criteria used to classify the investigated networks. In Section 3, we present and discuss the results of the model analysis. Section 4 gives a short summary and concludes with an outlook on future investigations.
2 Methods
The modeling framework used for this work is briefly recalled first. We then specify the stable steady states to be reproduced by each candidate model, and introduce some characteristic measures and a robustness measure.
2.1 Modeling framework
The framework is based on ordinary differential equations and makes use of only very few modeling assumptions which shall be explained next. As we are considering networks of TRs, the influence of such a regulator on its target gene has to be described mathematically. It is assumed that this influence can be described by monotonous activation and inhibition functions. The definitions of these functions are given next.
Definition 1
An activation (inhibition) function is a function () with and:

() is continuously differentiable,

and as
( and as ), 
() is monotonously increasing (decreasing).
Note, that commonly used kinetics such as MichaelisMenten or Hill kinetics are covered by this definition. Furthermore, it is assumed that each of the TRs undergoes a linear degradation. With this, the dynamics of a network of TRs are given by
(1) 
In this, is the vector of concentrations of the TRs and models the combined influence of all TRs in a network on the TR . In this study we only allow multiplicative combinations of activation and inhibition functions. As generally the parameters and the exact shapes of the activation and inhibition functions are not known, (1) only captures the interaction structure.
As mentioned in the introduction, our goal is to study all possible networks consisting of three nodes , . The interaction structure of such a network is represented by an interaction matrix
(2) 
To give an example, the interaction matrix
encodes the network
in which the parameters and the shapes of the activation and inhibition functions are not specified.
2.2 Specifications of stable steady states
We are specifically interested in the stable steady states of these networks as they determine the cell type. However, in this context we will not consider stable steady states as individual points but as forwardinvariant sets in the state space to account for biological variability. By doing so, only high and low levels of a TR are distinguished. To be more precise, we assume concentrations , and with to be known, and a concentration in the interval is considered as low. Equivalently, a concentration is considered as high. With this, a stable steady state of (1) is considered as a forwardinvariant hyperrectangular set with . For ease of notation, the ntuple of labels will be used to refer to this hyperrectangular set.
Let us now turn to our system of three TRs that provide a decision module in the cell differentiation process. In this setup, three steady states are of special interest: The progenitor state , a differentiated cell type , and a competing differentiated cell type . Each of these cell types is characterized by a specific 3tuple of high and low TRlevels . One can however find several plausible ways to characterize the three different cell types via levels of the three TRs. We assume the differentiated cell types to be characterized by a high level of the typespecific TR or , respectively. Regarding the progenitor state , we aim to compare two hypotheses (“specifications” of stable steady states, in the following for short SSSspecifications) against each other, which represent conceptionally distinct mechanisms:

corresponds to a progenitor factor, maintaining the progenitor state. is high in the progenitor state , whereas it is low in the differentiated states . It has to be downregulated in order to achieve cell differentiation. This means the required stable steady states are:
(3)  

corresponds to a differentiation factor, enabling differentiation. is low in the progenitor state , whereas it is high in the differentiated states . It has to be upregulated in order to achieve cell differentiation. This corresponds to requiring the stable steady states:
(4) 
Besides the specified stable steady states, we impose a further requirement on the candidate models: In order to investigate “truly interacting” nodes rather than assemblages of isolated nodes or subgraphs, we require each 3node network to be weakly connected.
2.3 Characteristic measures
Several characteristic measures are introduced to classify networks with respect to their interactions, including selfloops, and the types of these interactions (activating / inhibiting). The number of nonzero entries in the interaction matrix , which equals the zeronorm , is used as a measure of the network’s connectivity. To measure the preponderance of activating versus inhibiting interactions, we use the number of positive entries in minus the number of negative entries in , which we call the positivenegative weight.
The number of signinconsistent (signconsistent) loops is used to measure the degree of inconsistencies (consistencies, respectively) arising from network interactions. A loop is either a feedback loop from to itself, or a feedforward loop from to . Signinconsistencies arise when a feedback loop has negative sign, or the two paths of a feedforward loop from to have opposite signs. Signconsistent loops are feedback loops with positive sign, as well as feedforward loops with the two paths having the same sign. Similar motifs have been outlined e.g. in Ma’ayan et al. (2008).
To measure how strongly the connectivity is distributed among the three nodes, the maximum indegree is computed for each network:
(5) 
which is the largest rowsum, thus the norm .
2.4 Robustness measure
Since the goal of this study is to compare networks that can reproduce the desired cell types, we computed for each network structure a robustness measure that quantifies if and how well this structure can generate a set of forwardinvariant sets as specified in (3,4). We applied the algorithm introduced in Breindl et al. (2010) which is summarized next.
In order to explain the underlying idea of this method, let us use the symbol to denote both, activation and inhibition functions. These functions are indexed such that denotes the regulated TR, i.e., , and enumerates the regulators of . The number of regulators of is denoted as .
Also, let us measure a perturbation of a monotonous function as the norm of the difference of the original monotonous function and the perturbed monotonous function , i.e.,
(6) 
Given a set of desired forwardinvariant sets , , the goal of the method is to assign a robustness value to the interaction structure itself, i.e., to the matrix as in (2). This robustness value is defined such that it has the following three properties: (i) If a value is assigned to a system , there exists a realization of with monotonous functions , such that the desired sets , , are forwardinvariant. (ii) If no monotonous function is perturbed by more than , i.e., , it can be guaranteed that forwardinvariance of the sets , , is maintained. (iii) For every , there exist perturbed functions such that forwardinvariance of the sets , , is lost, and for at least one index it holds that .
The procedure to compute for a given network structure and specification , , is outlined next. To this end, two definitions are given next.
Definition 2
The 3tuple of pairs of positive real numbers
such that
and is called tube for
an activation function.
Equivalently, the 3tuple of pairs of positive real numbers
such that
and is called tube for an
inhibition function.
Definition 3
An activation function (inhibition function ) is said to satisfy a tube (), denoted as () if the following inequalities hold.
(7)  
(8)  
(9) 
This means that a tube for a monotonous function restricts the space where the graph of this function can evolve (see also Figure 1). Note that the values of the tubes are given by the specification of the forwardinvariant sets. It is now possible to formulate conditions on the values of the tubes such that the specified sets are forwardinvariant if all functions lie inside their tubes, i.e. , where can mean or as required, and the tubes are indexed as the according monotonous functions.
Proposition 1 (Breindl et al. (2010))
Given a set with and given tubes that satisfy the conditions
(10) 
where , and, with denoting the argument of ,
and
(11) 
where , and, with denoting the argument of ,
If , then the set is forwardinvariant for system (1).
Starting from this result, the robustness measure can now be introduced. As forwardinvariance for the specified sets can be guaranteed as long as , the goal is to find, for each tube , the function , which is best centered to the tube . In this context, best centered means that has to be perturbed more than any other function in order to violate at least one of the constraints (79) (denoted as ). The minimal perturbation of this best centered function in order to achieve a violation of the tube constraints is given as the result of the optimization problem
(12) 
with
(13) 
Then, the robustness measure is obtained by maximizing the smallest value of all tubes in the network over all tubes that satisfy Proposition 1.
Definition 4
Given a system (1) with unspecified activation and inhibition functions. The robustness measure for the system is defined as
(14)  
s.t.: 
This measure has the interpretation given at the beginning of this section. Furthermore, if there exists no realization of such that the desired sets are forwardinvariant, the optimization problem (14) is infeasible. It could furthermore be shown that for networks where the functions are products of activation and inhibition functions, problem (14) results in a convex optimization problem (Breindl et al., 2010).
3 Results
3.1 General observations
Constructing all possible 3node networks with interactions as specified in (2), i.e. each interaction having exactly one value , yields possible networks. Requiring the networks to be weakly connected leaves networks. Out of these candidate networks, the algorithm from Section 2.4 identifies models that meet the specification (S1), i.e. out of all networks; and models that are able to reproduce the specification (S2), i.e. out of all networks.
The distributions of the characteristic measures as outlined in Section 2 are shown in Fig. 2. Models for (S1) (Fig. 2 upper row, red) and models for (S2) (Fig. 2 upper row, green) show similar distributions of nonzero entries. The distribution of the positivenegative weight for models of (S1) is shifted to the left, i.e. toward more negative entries, whereas for models of (S2) it is slightly shifted to the right, i.e. toward more positive entries. These differences are in accordance with the number of inconsistent loops which require negative entries: Models for (S1) show a distribution of inconsistent loops up to high numbers and only intermediate amount of consistent loops, whereas models for (S2) tend to low numbers of inconsistencies and more consistent loops. Also, models for (S1) (red) are found for all numbers of inconsistent loops, whereas models for (S2) (green) are restricted to up to 8 inconsistent loops. In contrast, the (S1)models are more restricted to certain numbers of consistent loops than the (S2)models are.
Summarizing, there are less networks supporting the hypothesis of a progenitor factor (S1), and they show a tendency toward negative (inhibiting) interactions and hence inconsistent loops. The number of networks providing models for the hypothesis of a differentiation factor (S2) is higher. On average, these networks have about equal number of positive and negative entries, thus reducing the amount of inconsistent loops in favor of consistent loops.
3.2 Structural properties
Regarding the network structure, we found that all models that can reproduce a SSSspecification can be classified leading to sufficient and necessary conditions on the interaction structure. These findings are summarized in the following observations.
Observation 1
Let be an interaction matrix (2). If the network represented by is such that each node , together with its incoming links is contained in Fig. 3, then there exists a realization of that can reproduce the SSSspecification (S1). Also the other direction holds: For each weakly connected 3node network that can be assembled from nodes together with its incoming links contained in Fig. 3 there exists a realization that can reproduce the SSSspecification (S1).
An equivalent observation can be made for the network that can reproduce the SSSspecification (S2).
Observation 2
Let be an interaction matrix (2). If the network represented by is such that the node and the nodes , together with their incoming links are contained in Fig. 3, then there exists a realization of that can reproduce the SSSspecification (S1). Also the other direction holds: For each weakly connected 3node network that can be assembled from nodes and , , together with their incoming links contained in Fig. 3 there exists a realization that can reproduce the SSSspecification (S2).
From Fig. 3 it can also be explained why the positivenegative weight of models for SSSspecification (S1) tends toward a larger number of negative entries: All interactions between two different nodes in the network need to be negative (inhibiting). Only selfloops are allowed to be activating. For SSSspecification (S2) on the other side also activating links from to the other nodes are possible.
Using these “building blocks” as listed in Fig. 3 it is possible to construct all possible 3node models that can, for an appropriate choice of activation and inhibition functions, reproduce the SSSspecification (S1) or (S2). To see possible applications of our approach, note that a switching module for mesenchymal stem cell differentiation based on literature data and analyzed in detail in Schittler et al. (2010) satisfies specification (S1) and can indeed be composed from the building blocks presented in Fig. 3. By applying the algorithm from Section 2.4, it is in principle possible to compute the building blocks for any SSSspecification. The algorithm is also directly applicable to larger networks.
3.3 Robustness properties
The mean robustness value for all characteristic measures is shown in Fig. 2 (lower row). Two important observations can be made: There is a clear negative correlation of with increasing nonzero entries, and a clear positive correlation of with increasing positivenegative weight. Both observations are true independent of the specific SSSspecifications (red, green). In contrast, no clear correlation of versus the number of inconsistent or consistent loops can be detected.
It is stated by Kwon et al. (2008) that perturbations on nodes involved in fewer (more) than the average number of feedback loops have a lower (higher, respectively) impact on the overall network robustness. Thus, we investigated whether there is a similar relationship between the and the robustness value . The results of our analysis are shown in Fig. 4. Among all (S1)reproducing models, only four different robustness values have been assigned while for the (S2)reproducing models seven different robustness values have been computed. Furthermore, the (S1)models are more abundant in higher robustness values, whereas the majority of (S2)models has a lower robustness value (cf. Fig. 4). That is, models for hypothesis (S2) require a finer “tuning” than models for (S1), biologically pointing toward more vulnerable regulatory networks.
Yet, despite these differences, the similarities between the two distributions are eminent which suggests a strong relation between the network’s structural properties and their robustness: Independent of the specification (S1) or (S2), networks with a higher tend to have lower values. This means that, the more incoming interactions a node has, the more vulnerable the network gets to perturbations of these interaction. This is in accordance with the findings of Kwon et al. (2008) about network structure and their notion of robustness.
4 Conclusion
The main focus of this work was to analyze which networks consisting of three transcriptional regulators are in principle able to reproduce a set of required stable steady states, to search for interaction patterns that are necessary in order to meet these requirements, and to compute and compare the robustness properties of these networks. Still, the investigations raise new questions and possible future work, which we point out briefly in this concluding outlook.
Signinconsistencies make a system nonmonotone, and the more inconsistent loops are contained in a network, the “farther” the system is from monotonicity. Ma’ayan et al. (2008) argue that intracellular systems will be “closetomonotone”. Kwon et al. (2008) report that a higher number of positive feedback loops and a smaller number of negative feedback loops result in a higher robustness. These arguments are in accordance with the result here that positive (negative) entries increase (decrease) the robustness of a network. But in contrast to the number of positive / negative entries, the number of inconsistent and consistent loops does not show a clear effect on the robustness as defined here. What is more, inconsistent loops are even required for some SSSspecifications to be met. These studies use different notions of robustness, and their interrelation should be worth further investigations.
Furthermore, one might argue that especially in cell differentiation, stimulus inputs and transitions between stable steady states rather than just the existence of stable steady states play an important role. However, although these aspects were not considered in this study, the approach in this paper provides a first selection step of networks with respect to one necessary condition (out of several), namely by checking whether the required stable steady states can be reproduced. These aspects provide plenty of opportunities for separate investigations, and we hope to address these issues in future work.
Achknowledgement
The authors would like to thank for financial support from the German Research Foundation (DFG) within the Cluster of Excellence in Simulation Technology (EXC 310/1) at the University of Stuttgart, and from the German Federal Ministry of Education and Research (BMBF) within the SysTec program (grant nr. 0315506A), and from the Helmholtz Society within the CoReNe project. D.S. acknowledges financial support by The MathWorks Foundation of Science and Engineering.
References
 (1)
 Breindl et al. (2010) C. Breindl, S. Waldherr and F. Allgöwer. A robustness measure for the stationary behavior of qualitative gene regulation networks. Proceedings of the 11th symposium on computer applications in biotechnology (CAB), pp. 36–41, 2010.
 Chaves et al. (2008) M. Chaves, T. Eissing and F. Allgöwer. Bistable Biological Systems: A Characterization Through Local Compact InputtoState Stability. IEEE Transactions on Automatic Control, Special Issue on Systems Biology, vol. 53, pp. 87–100, 2008.
 Foster et al. (2009) D.V. Foster, J.G. Foster, S. Huang and S.A. Kauffman. A model of sequential branching in hierarchical cell fate determination. Journal of Theoretical Biology, vol. 260, pp. 589–597, 2009.
 Huang et al. (2007) S. Huang, Y.P. Guo, G. May and T. Enver. Bifurcation dynamics in lineagecommitment in bipotent progenitor cells. Developmental Biology, vol. 305, pp. 695–713, 2007.
 Kwon et al. (2008) Y.K. Kwon and K.H. Cho. Quantitative analysis of robustness and fragility in biological networks based on feedback dynamics. Bioinformatics, vol. 24, pp. 987–994, 2008.
 Ma’ayan et al. (2008) A. Ma’ayan, A. Lipshtat, R. Iyengar and E.D. Sontag. Proximity of intracellular regulatory networks to monotone systems. IET Systems Biology, vol. 2, pp. 103–112, 2008.
 MacArthur et al. (2008) B.D. MacArthur, C.P. Please and R.O.C. Oreffo. Stochasticity and the molecular mechanisms of induced pluripotency. PLoS ONE, vol. 3, e3086, 2008.
 MacArthur et al. (2009) B.D. MacArthur, A. Ma’ayan and I.R. Lemischka. Systems biology of stem cell fate and cellular reprogramming. Nature Reviews Molecular Cell Biology, vol. 10, pp. 672–681, 2009.
 Narula et al. (2010) J. Narula, A.M. Smith, B. Gottgens and O.A. Igoshin. Modeling reveals bistability and lowpass filtering in the network module determining blood stem cell fate. PLoS Computational Biology, vol. 6, e1000771, 2010.
 Peltier and Schaffer (2010) J. Peltier and D.V. Schaffer. Systems biology approaches to understand stem cell fate choice. IET Systems Biology, vol. 4, pp. 1–11, 2010.
 Roeder and Glauche (2006) I. Roeder and I. Glauche. Towards an understanding of lineage specification in hematopoietic stem cells: A mathematical mode for the interaction of transcription factors GATA1 and PU.1. Journal of Theoretical Biology, vol. 241, pp. 852–865, 2006.
 Schittler et al. (2010) D.Schittler, J. Hasenauer, F. Allgöwer, and S.Waldherr. Cell differentiation modeled via a coupled twoswitch regulatory network. Chaos, vol. 20(4):045121, 2010.
 Xiong and Ferrell (2003) W. Xiong and J.E. Ferrell, Jr.. A positivefeedbackbased bistable memory module that governs a cell fate decision. Nature, vol. 426, pp. 460–465, 2003.