FEEDBACK REGULATION IN MULTISTAGE CELL LINEAGES

Studies of developing and self-renewing tissues have shown that differentiated cell types are typically specified through the actions of multistage cell lineages. Such lineages commonly include a stem cell and multiple progenitor (transit amplifying; TA) cell stages, which ultimately give rise to terminally differentiated (TD) cells. In several cases, self-renewal and differentiation of stem and progenitor cells within such lineages have been shown to be under feedback regulation. Together, the existence of multiple cell stages within a lineage and complex feedback regulation are thought to confer upon a tissue the ability to autoregulate development and regeneration, in terms of both cell number (total tissue volume) and cell identity (the proportions of different cell types, especially TD cells, within the tissue). In this paper, we model neurogenesis in the olfactory epithelium (OE) of the mouse, a system in which the lineage stages and mediators of feedback regulation that govern the generation of terminally differentiated olfactory receptor neurons (ORNs) have been the subject of much experimental work. Here we report on the existence and uniqueness of steady states in this system, as well as local and global stability of these steady states. In particular, we identify parameter conditions for the stability of the system when negative feedback loops are represented either as Hill functions, or in more general terms. Our results suggest that two factors – autoregulation of the proliferation of transit amplifying (TA) progenitor cells, and a low death rate of TD cells – enhance the stability of this system.


Introduction
Tissue development is charged with the difficult task of specifying correct types of terminally differentiated (TD) cells (whose function typically characterizes the tissue), at the same time it must specify correct numbers and proportions of the various TD and progenitor cell types that compose the tissue. Studies have shown that multistage cell lineages, which typically include a multi-potent stem cell that gives rise to successively more differentiated progenitor cell types (transit amplifying cells; TA cells), underlie the specification of TD cells in many tissues [22,28,37]. In some tissues, control and coordination of development is mediated by feedback signals, which regulate proliferation and differentiation at specific lineage stages [38]. Although often established during embryonic development, the continued presence of feedback controls in adult organisms is believed to enable tissues to respond robustly to injury, rapidly and accurately regenerating appropriate cell numbers and proportions [11,18].
The olfactory epithelium (OE) of the mouse is a useful model system for understanding how feedback regulation controls tissue development and regeneration [26]. The OE is a neuroepithelium similar in structure to the early embryonic central nervous system, but simpler in that it possesses only one major type of neuron, the olfactory receptor neuron (ORN). ORNs, the TD cells of the OE, extend apical (dendritic) cytoplasmic processes to the surface of the OE; these dendrites are in direct contact with the external environment via the nasal cavity. ORNs also extend axonal cytoplasmic processes basally; these axons traverse the basal lamina of the OE, grow through the connective tissue stroma that underlies the OE, and then make contact with the olfactory bulb of the forebrain; olfactory sensory information is conducted from the periphery to the central nervous system via the axons of the ORNs [4,10]. The OE is subject to constant environmental impact due to its exposed position in the body, and therefore must regenerate ORNs continually, throughout the lifetime of the organism [7].
Studies of mouse OE have supported the view that ORNs derive from a multistage lineage with three proliferating cell types (reviewed in [5]): (1) self-renewing stem cells, which give rise to (2) TA progenitors that can be identified by their expression of the proneural gene Mash1. Mash1-expressing progenitors give rise to (3) immediate neuronal precursors (INPs), which comprise a second TA population, whose members expresses Ngn1, a proneural gene distinct from Mash1. INPs divide to give rise to daughter cells that undergo terminal differentiation into (4) ORNs. The four stages in the ORN lineage are depicted in Figure 1.
It has long been known that cell proliferation in the OE is rapidly up-regulated when death of a large number of ORNs is induced (reviewed in [7]). More recently, in vitro observations showed that proliferation and generation of ORNs by cultured OE neuronal progenitors is inhibited by the presence of large numbers of ORNs [24]. These results suggest that ORNs produce signals that normally inhibit neurogenesis, such that the elimination of such signals upon ORN loss leads to an increase in ORN production. Subsequently, it has emerged that growth and differentiation factor 11 (GDF11), and ActivinβB, two members of transforming growth factor-β (TGFβ) superfamily, are strong candidates for such factors [13,18,31,32,38]. In particular, in vitro and in vivo experimental data indicate that GDF11, which is produced by INPs and ORNs, inhibits the proliferation of INPs, while ActivinβB, which appears to be produced by all cells in the lineage, inhibits proliferation of the stem cell [13]. There is evidence that both molecules slow the rates at which their target cells divide and increase the probability that the products of those divisions differentiate into cells at the next lineage stage [13,18]. Other types of molecules, such as Follistatin (an inhibitor of GDF11 and ActivinβB), and bone morphogenetic proteins (BMPs), may also regulate the lineage [13,18,31,32,38]. It has been proposed that a multiplicity of feedback mechanisms and target lineage stages is required for tight control of cell population sizes and regeneration rates [18].
Evidence for feedback strategies analogous to those used in the OE has been found in many tissues and organs. For example, muscle cells produce myostatin/GDF8, a molecule closely related to GDF11, which exerts an inhibitory effect on myogenesis [19,23]. In well-studied hematopoietic and epithelial lineages, other proliferative feedback mechanisms have been proposed or described [1,29].
Various models have been proposed for different cell lineages. One of the influential works by Tomlinson and Bodmer [36] was a study of tumorigenesis based on discrete cell models, later followed by Johnston et al. [16] for cell population dynamics in the colonic crypt and colorectal cancer. Continuum models were also proposed for the same biological systems [15,16]. Recently, using stochastic and discrete modeling, Mangel and Bonsall [21] studied variability and replacement in multistage cell lineages, and Glauche et al. [12] considered lineage specification of hematopoietic stem cells. Most of these works primarily used numerical simulations for model exploration. Other biological systems involving multiple cell types, such as hematopoiesis for myelogenous leukemia [8,9] and T-cells in autoimmune diabetes [20], were studied using dynamic systems approach. On mathematical analysis, a linear stability analysis was performed in [3,30] for a class of cyclic systems and a cellular control process with positive feedback; a global stability analysis was carried out for monotone cyclic systems [34] using the Poincaré-Bendixon Theorem in multi-dimension based on the discrete Lyapunov function method [33]. In particular, global stability may be studied using ω-limit set [2], Lyapunov function method [3] or standard Poincaré-Bendixon Theorem in two-dimension [14].
In this paper, we study the OE lineage (and analogous multistage cell lineages) using a continuum dynamic model. In the system, the feedbacks considered are negative feedbacks of the type supported by experiments [13,38]. First, we study a system with a constant stem cell population, and two feedbacks on proliferation of TA cells. For this simplified system, the existence and uniqueness of the steady state, and local and global stability of the steady state are analyzed. It is found that a unique positive steady state exists, and both self-regulation of the proliferation of TA cells, and a low death rate of TD cells, enhances global stability of the steady state. Next, we consider more general lineages involving potentially dynamic stem, TA and TD cells, with a general form of negative feedback. For existence and stability of three possible steady states, conditions on the feedback are derived, with detailed analysis on parameter ranges when the feedbacks take the form of Hill functions. We find that selfregulation on the proliferation of stem cells and TA cells, as well as a large ratio of stem cell cycle length to TA cell cycle length, can enhance local stability. While the existing models typically used computational tools to explore the mechanisms of the cell population control, our work utilizes mathematical analysis to explain the functions of different feedback loops. This paper is organized as follows. In Section 2, we introduce the mathematical model; in Section 3, we study the model with a lineage in which the size of the stem cell population is constant (not influenced by feedback); in Section 4, we consider the more general case in which the populations of all three cell types are dynamically regulated; in Section 5, we discuss and summarize.

A cell lineage model with endogenous negative feedback
The cell lineage system is assumed to have three types of cells: stem cells, transit amplifying (TA) cells and terminally differentiated (TD) cells, with the populations at time t denoted by χ 0 (t), χ 1 (t) and χ 2 (t) respectively. In the model, the stem cells and TA cells proliferate with a probability p i and differentiate with its complementary probability 1 − p i , where 0 < p i < 1, i = 0, 1. The TD cells undergo apoptosis (programmed death) with rate δ. ζ i (i = 0, 1, 2) represents the cell cycle length divided by log 2 for each type of cell.
In many systems, such as the OE lineage, the proliferation of cells may be inhibited by substances produced by the cells at the same or later stages in their lineage. By assuming the amounts of substances are proportional to the sizes of the cell populations that produce the substances, the proliferative probabilities become functions of the cell populations, p 0 = p 0 (χ 0 , χ 1 , χ 2 ) and p 1 = p 1 (χ 1 , χ 2 ). Figure 2 shows a schematic relationship and possible negative regulation among different types of cells in the cell lineage.

A model with a constant stem cell population
A special case of system (4)- (6) is that the stem cell population does not change in time and the proliferation of the stem cells is not regulated. This can be achieved if and only if, in equation (4), p 0 (χ 0 , χ 1 , χ 2 ) = 0.5 with a positive χ 0 (see Figure 3). In this case, half of the stem cells progeny remain stem cells and the other half differentiate; as a result the stem cell population remains unchanged. The three-equation system (4)-(6) therefore becomes a two-equation system: a given positive constant. The reduced system (7)-(8) can be regarded as a two-cell lineage with a constant supply. In this model, only two negative feedbacks occur at the proliferation probability of the TA cells, and this simplicity allows a comprehensive analysis for the existence and uniqueness of the steady state, and their local and global stability of the system (7)- (8). As a result, roles of the two negative feedbacks are better understood.

Existence and uniqueness of the steady state
Because χ 1 and χ 2 represent the size of cell populations, only non-negative steady state solutions are considered in this paper.
The steady state equations of system (7)-(8) can be obtained by setting their left-hand side zero. Adding the two steady state equations yields (9) and substituting equation (9) into the first steady state equation yields (10) The steady state can be found by first solving equation (10) for χ 1 and then using (9) to derive χ 2 . The number of the steady states of system (7)-(8) equals the number of roots in equation (10), and the function p 1 determines the number of steady states.

3.1.1.
The steady state is unique if the minimal proliferation probability is smaller than its corresponding differentiation probability-As shown below, a necessary and sufficient condition for the existence and uniqueness of a positive steady state of system (7)-(8) is: (11) With the assumption that p 1 is non-increasing and using (11), one-half is an upper bound for the minimum of p 1 if and only if the unique postive steady state exists. Biologically, this condition means that the inhibition on the proliferation probability, from either the TA cells or the TD cells, must be strong enough such that at a certain time, more than half of the TA cells differentiate to the next stage per cell cycle.
To prove (11), consider the function (12) With δ > 0 and ζ > 0, it is straightforward to show , and lim χ1→∞ K(χ 1 ) = −∞ if and only if (11) is satisfied. By continuity of K, there exists a positive such that , if and only if (11) is true; therefore, (11) is necessary and sufficient for the existence of a positive root for K. Next, we claim that is the only positive root of K = 0. Consider (13) the first term is strictly decreasing and the second term is monotonically decreasing in χ 1 so is strictly decreasing in χ 1 . Therefore, has at most one root for χ 1 > 0. Since and K(χ 1 ) share same positive roots, K(χ 1 ) has at most one root. It implies that (11) is a necessary and sufficient condition for the existence of the steady states, and if a steady state exists, it must be positive and unique.

A system with negative feedback represented by Hill functions has a unique steady state-The feedback in biological models is often represented in terms of
Hill functions [25]. Two possible Hill function forms of p 1 are (14) and (15) with the feedback coefficients g 1 , g 2 ≥ 0 and at least one of them being positive, p̄1 > 0 the maximal differentiation probability, and m, n denoting the Hill exponents. Hence, p 1 (χ 1 , χ 2 ) in (14) and (15) satisfies the conditions (11) and the system has a unique positive steady state. In the rest of this section, when the Hill function form of p 1 is referenced, only the first expression (14) is considered. However, the same results can be obtained for (15) through similar arguments.

Local stability analysis
Assume the system has a non-negative steady state , (i.e. p 1 satisfies (11)). Denoting , the Jacobian matrix evaluated at the steady state is (16) The characteristic equation of the Jacobian matrix (16) is (17) Then, the steady state χ* is locally asymptotically stable if and only if both the roots of equation (17) have negative real parts, namely, if and only if (18) and (19) Because is a steady state, by equation (10), is negative. Also, by the fact that and , inequality (19) is always satisfied. Therefore, equation (18) dictates the local stability. (18), here we consider p 1 in the Hill function form (14). The partial derivative of the function p 1 (χ 1 , χ 2 ) is (20) Because the first two terms are negative, the inequality (18) is proved if (21) By equation (9), we have , so equation (21) is true if which is equivalent to

The inhibition by TD cells determines stability-To find p 1 satisfying
When the Hill exponent n = 1, by we obtain When n = 2, by p̄1 < 1 we obtain Therefore, when n = 1, 2, condition (21) (and hence (18)) is satisfied, the steady state is unique and locally asymptotically stable.

Negative regulation from TA cells increases the range of parameters for
local stability-For n ≥ 3, the above analysis cannot be carried out directly. In fact, unstable steady states for n ≥ 3 are observed through direct numerical simulations. Here, we study the linear stability for two cases: 1) n = 3 and g 1 = 0, a case with feedback only from the TD cells; and 2) n = 3, m = 1, a case with both feedbacks. Figure 4 plots the stability regions for p̄1 and δ at four different g 2 , a parameter measuring the feedback coefficient of the TD cells and thus the gain of the feedback. The area of unstable region increases as g 2 decreases, as shown in Figure 4. A low death rate for the TD cells leads to stability of the steady state without any oscillations. A high death rate for the TD cells requires a low maximal proliferation probability p1 to guarantee a stable steady state without oscillations (because the perturbed solutions away from the steady state converge to the steady state in the Stable (2C) region with oscillations). When both the death rate of TD cells and the proliferation probability of TA cells are high, the steady state becomes unstable. This instability can be removed when a negative feedback from TA cells to its own proliferation probability is included (see Figure 5). This suggests both negative feedbacks help improve the stability of the steady state.
Finally, a limit cycle exists when the steady state is unstable for the systems in Figure 4 and Figure 5. This can be proved by Poincaré-Bendixson Theorem [27] using the fact that the real parts of the two eigenvalues of the Jacobian matrix at the unique steady state have the same sign. The detailed proof is not presented here.

Global stability analysis
Local stability analysis studies the behavior of a system when the solution is a small perturbation away from a steady state. For an arbitrary perturbation of the steady state solution, a global stability analysis [35] is needed. The following theorem gives a condition on p 1 (χ 1 , χ 2 ) for the global stability of a positive steady state of system (7)-(8).
The proof of the theorem is a direct application of Poincaré-Bendixson Theorem [27] and it is sketched as follows with details presented in Appendix A.1.
3. If the system satisfies (22), the steady state solution is locally stable.

4.
Using the Poincaré-Bendixson Theorem and the above three statements, one can prove the global stability of the steady state.
Based on the theorem, a simpler condition for the global stability can be obtained: Corollary 1-If p 1 satisfies (23) then it satisfies the condition (22). By Theorem 3.1, the positive steady state is globally stable.
Because the stability analysis is on the long-time asymptotic behavior of the solution, Theorem 3.1 and Corollary 1 remain true if the two variables χ 1 , χ 2 are replaced by lim inf t→∞ χ 1 (t) and lim inf t→∞ χ 2 (t).

Low death rate of TD cells leads to global stability-It
can be shown that the steady state is globally stable if (24) The proof is presented in the Appendix A.2.
One immediate result from (24) is that when n = 1 the steady state is always globally stable regardless of choice of other parameters.
When n > 1, one possibility of increasing the stability region for p̄1 is decreasing the death rate δ. Another possibility is to increase the inhibition strength g 2 , that is, to increase the level of the negative feedback from the TD cells to the proliferation probability of the TA cells. Both results are consistent with the linear stability analysis shown in Figure 4 and Figure 5. For the OE system, the δ is usually very small [26], and this fact is consistent with the mathematical strategy to enhance the stability region.
To better analyze (24), one may use a more strict yet simpler inequality to replace (24): (25) which says p̄1 is bounded by a decreasing function of the Hill exponent for the negative feedback from the TD cells on proliferation probability of TA cell. For example, n = 2 implies p̄1 < 2/(1.5) 2 = 8/9 for existence of a globally stable steady state while n = 3 implies p̄1 < 3/2 2 = 0.75, a more restricted condition than the case for n = 2. where (i = 0, 1, 2) is positive. We use S1, S2 and S3 to represent each type of steady state, respectively. By direct examination of the the full system (4)- (6), it is easy to show that the five other types of solutions, such as , cannot be the steady state of (4)- (6).

Full model with stem, TA and TD cells
In the rest of the section, we first discuss each type of steady state, without specifying any functional form for the negative feedback, and present the conditions for the existence and stability for that type of steady state. Then, by using Hill functions for the feedback with Hill exponent taken to be one, we are able to classify four different parameter regimes based on the maximal proliferation probabilities of the stem and TA cells to identify all steady states in each parameter regime, and provide stability properties for some of the steady states. Finally, we present a linear stability analysis for the steady state S3.

• S1 steady state
It is easy to see from the system (4)-(6) that S1, a case with zero population for each cell type, is always a steady state for any parameters. The interesting question is on its stability. Because the eigenvalues of the Jacobian at S1 are (26) S1 is locally stable if and only if p 0 (0, 0, 0) < 0.5 and p 1 (0, 0) < 0.5. This shows that given small populations of each cell type, the populations will decrease and vanish eventually if the proliferation probabilities of stem cells and TA cells are both lower than their differentiating rates.
Regarding global stability, S1 is a global attractor if p 0 (χ 0 , χ 1 , χ 2 ), p 1 (χ 1 , χ 2 ) < 0.5 for all (χ 0 , χ 1 , χ 2 ), because the right-hand side of equation (4)-(6) will eventually become negative under this condition. Since we assume that p i are decreasing functions, the condition of global stability is equivalent to that of local stability, namely, p 0 (0, 0, 0) < 0.5 and p 1 (0, 0) < 0.5. This suggests that in order to avoid the S1 steady state, the proliferation probability for either the stem cells or the TA cells must be larger than its corresponding differentiation probability at some time during the evolution.

• S2 steady state
In the S2 steady state, there are no stem cells and the cell populations are maintained by the proliferation and differentiation of TA cells. At S2, the proliferation probability of TA cells has to be the same as the differentiation probability as seen from the equation (4)- (6). Mathematically, the existence of S2 thus depends on whether p 1 can achieve 0.5 at some nonzero TA and TD populations. Because p 1 is assumed as a nonincreasing function, a necessary condition for existence of S2 is that the maximal proliferation probability (27) and the feedback is strong enough such that p 1 (χ 1 , χ 2 ) ≤ 0.5 at some (χ 1 , χ 2 ). If p 1 is assumed to be a strictly decreasing function in either χ 1 or χ 2 , then there exist a unique S2.
The eigenvalues of the Jacobian matrix at S2 are and the two roots of the quadratic equation So S2 is locally stable if (28) These conditions may be achieved by allowing a longer cell cycle for TA cells than stem cells along with a low stem cell proliferation probability.

• S3 steady state
By first adding the steady state equations of system (4)-(6), one obtains (29) Substituting (29) into the first two equations of (4)-(6), one obtains (30) The number of S3 steady states is determined by the number of roots of system (30)- (31). From the discussion of S1, S3 exists only if the proliferation probability of stem cells is larger than its differentiation probability at some time. One of the necessary condition for existence of S3 is (32)

Parameter regimes for existence of the three types of steady states
In order to derive more specific conditions for existence and stability of the three steady states, we use the following Hill functions for the proliferative feedback: (33) (34) where k i (i = 0, 1, 2) ≥ 0 and g i (i = 1, 2) ≥ 0 with at least one of them being positive, and 0 < p0, p̄1 < 1 are constants. Based on the values of p̄0 and p̄1, we classify the system into four regimes with results summarized in Table 1.
• p̄0 ≤ 0.5 and p̄1 ≤ 0.5: S1 is globally stable and there is no other non-negative steady state.
Because the sign of the right-hand side of equation (4)-(5) becomes negative as t → ∞ under this condition, χ 0 , χ 1 tend to zero. Consequently, so is χ 2 by (6). Because of (27) and (32), both S2 and S3 cannot be the steady states. This implies that S1 is a unique and globally stable steady state of the system.
• p̄0 ≤ 0.5 and p̄1 > 0.5: S1 is unstable and there exists a unique globally stable S2 but no other non-negative steady state. S1 is locally unstable by using (26). There exists a unique S2 steady state given by (35) The eigenvalues of the Jacobian at this S2 steady state can be easily shown to be negative using a similar discussion as in Section 3.2. So S2 is locally stable. This steady state can further be proved globally stable for any initial conditions (χ 0 , χ 1 , χ 2 ) when χ 0 > 0 or χ 1 > 0, using a proof similar to that in Section 3.3. In addition, it can be shown that for an initial condition (0, 0, χ 2 ), the solution converges to S1, meaning that if there are only TD cells which cannot proliferate and only undergo death, the whole population will vanish eventually.
• p̄0 > 0.5 and p̄1 ≤ 0.5: S1 is unstable and there exists a unique S3 but no other non-negative steady state. S1 is locally unstable by linear stability analysis, but it is stable if the initial condition is of the form (0, 0, χ 2 ), as explained in the first case. With 2p̄1 ≤ 1, there exists no non-negative S2. The existence of a unique steady state S3 is proved in Appendix A. 4. The stability of S3 is discussed in subsection 4.3.
• p̄0 > 0.5 and p̄1 > 0.5: S1 is unstable and either there exist a unique stable S2 and no other non-negative steady state, or a unique unstable S2 and a unique S3 if (36) Similar to the previous cases, S1 is locally unstable because p̄0 > 0.5, but it is stable with initial conditions (0, 0, χ 2 ). If (36) holds, there is a unique, unstable S2 and a unique S3. On the other hand, if (37) then S2 is locally stable, and steady state S3 does not exist. The proof is given in Appendix A.4. When the left-hand side of (36) equals the right-hand side, S3 does not exist, and in our numerical simulations, S2 is found to be stable. On the plane, a bifurcation takes place when the inequality in (36) becomes equality, as shown in Figure 6 with two sets of parameters.

Stability conditions for the S3 steady state
When a unique S3 steady state exists, we consider the cubic characteristic polynomial of the Jacobian at S3, (38) where . Let λ 1 , λ 2 , λ 3 be the three roots of C(λ), the S3 steady state χ* is locally asymptotically stable if and only if These conditions are equivalent to that λ 1 , λ 2 and λ 3 are all negative, and its proof is shown in the appendix A.5. From the property of the cubic polynomial C(λ), one can easily obtain: So (39)-(41) give a set of conditions on the parameters for the local stability of χ*. For more specific results, we consider several Hill feedback forms with the corresponding stability results described below and summarized in Table 2.
• p 0 (χ 2 ) is a Hill function and p 1 is a constant less than 0.5 The proliferation probability of stem cells is assumed in terms of Hill function: (42) with p̄0 > 0.5 and k 2 > 0.
This function and parameters satisfy the conditions (39) and (40), and the condition (41) is equivalent to (43) When n = 1, implies that the steady state is locally stable for any p̄0.
In general, the range of p0 which allows local stability increases as increases, that is, if the cell cycle of the TA cell is shorter than that of the stem cell, or the death rate of TD cells is larger, the S3 steady state is stable for a larger range of proliferation probabilities.
• p 0 (χ 1 , χ 2 ) is a Hill function with exponent one, and p 1 is a constant less than 0.5 The proliferation probability of stem cells is assumed as (44) with p̄0 > 0.5 and k 1 > 0.
The function and parameters satisfy the conditions (39) and (40), and by (41), we derive a sufficient condition for the local stability: (45) which is similar to the first example where p 0 only depends on χ 2 . Also, if the condition (45) is not satisfied, S3 can be stable if (46) Therefore, if ζ is small (the cell cycle of the TA cell is much shorter than that of the stem cell) or is small (the feedback strength from TA cells is much stronger than that from TD cells), S3 is locally stable.
Assume a unique S3 exists, using a similar proof in Section 3.2, we obtain (47) Consequently, the conditions (39)-(41) hold for any parameters, and the unique S3 is locally stable unconditionally.

Conclusion and discussion
In this paper, mathematical models for the Olfactory Epithelium (OE) lineage were developed and analyzed. Three different types of cells -stem cell, transit amplifying (TA) cell, and terminally differentiated (TD) cell, were considered, and the dynamic modulation of proliferation by secreted signaling molecules (such as GDF11 or ActivinβB) was modeled as negative feedback between lineage stages. For a simplified system with a constant stem cell population size, the conditions for the existence, uniqueness, local and global stability of the steady state were derived for negative feedback represented in terms of both general decreasing functions and Hill functions. It was found that a unique positive steady state exists, and both self-regulation of the proliferation of TA cells, and a low death rate of TD cells, can enhance global stability of the steady state. This might explain in the OE system why GDF11, which is secreted by TA and TD cells, inhibits the proliferation of TA cells. For the full system in which all three cell populations can vary in time, conditions on existence and linear stability of three possible steady states were obtained for feedback represented by Hill functions. It was found that self-regulation of the proliferation of stem cells and TA cells can enhance local stability. In addition, having the cell cycle lengths of stem cells be larger than those of TA cells also makes the local stability region larger.
The model developed in this study can be used as a generic model of cell lineages in which three sequential lineage stages cells coupled through negative feedback. Our analysis of the steady states and their stability in such a system could improve understanding of interactions and regulation among different cell types, and the dynamics and stability of different cell populations, and lead to new experimental approaches.
In OE, as in other epithelial tissues such as the epidermis and the colonic mucosa, cells at different lineage stages are usually stratified along basal-apical axis, leading to the expectation that the distributions of feedback signaling molecules (e.g. GDF11) will also be non-uniform [18]. In the future, it will also be important to incorporate spatial inhomogeneity and spatial dynamics into models such as the ones described here. Although numerical simulation will likely be the most appropriate tool for exploring such models, the analytical results derived in this paper will provide insights into the choice of parameters, and can serve as a guide for systematic computational exploration.
This work was supported by the NSF/NIH initiative on Mathematical Biology (R01 GM75309 and R01 GM67247), and NIH grants P50 GM76516 (to ADL, ALC, QN and FW) and R01 DC03583 (to ALC). KKG was supported by a training grant in neuroscience (T32 NS07444-08), and the UCI medical scientist training program (T32 GM08620).

d.
If χ 2 ≫ 0, then goes to some negative constant because in (c) we showed that χ 1 is bounded. So χ 2 will decrease when χ 2 ≫ 0. With (b), χ 2 is bounded in non-negative region.
3. If the system satisfies (22), then the steady state solution is locally stable. Since is the steady state, the right-hand side and we get If (22) holds, then By and (18), the steady state is locally stable.

A.2. Prove that (24) implies global stability of system (7)-(8)
Proof. Substitute p 1 as given by (14) into the left-hand side of (23) and get Inequality (23) will be satisfied if The aim here is to find conditions for p̄1 such that (51) holds for χ 2 > lim inf t→∞ χ 2 (t).
First, we can prove that the unique steady state is globally stable if (25) holds: The condition (51) can be simplified by taking y = (g 2 χ 2 ) n and considering the function (53) which assumes the minimal value at . Thus, if (52) is true, then h(y) ≥ h(y min ) > 0 for any y ≥ 0, equivalent to (51), then the steady state is globally stable.
The upper bound in (52) can be relaxed due to the fact that only the asymptotic behavior of the solution is relevant to global stability, so (51) has to be true only when t is very large. Hence, it is not necessary to require h(y) > 0 for all y ≥ 0 in (53); instead, we only need to find a lower bound ỹ, associated with the asymptotic behavior of χ 2 , with h(y) positive only for y ≥ ỹ. The derivation of the upper bound in (24) uses the following lemma to be proved in A.3: Suppose x(t), F(t) and G(x(t)) are functions satisfying and the following properties:

G(x) is continuous and strictly increasing in (−∞, ∞).
Then, (54) With this lemma and p 1 in Hill function form (14), we get by equation (7) and (8), and , it is easy to verify that (55) and the two properties in the hypotheses are satisfied, and hence by Lemma A.1 (57) Similarly we take x(t) := χ 2 (t), and G(χ 1 ) := δχ 2 , and apply Lemma A.1 can be applied with (57) to get (58) In (53), h(y) > 0 if y > n − 1, and by (58), the unique steady state is globally stable if (59) Combining (52) and (59), we proved that the unique steady state is globally stable if (24) is true.

A.3. Proof of Lemma A.1
Proof. Suppose (54) is not true, that is, we would like to prove Lemma A.1 by contradiction.
By property (2), G(x(t 1 )) < G(x 0 ) < G(x) and by the definition of x, we get (62) Therefore, there exists t 1 > t 0 such that (63) Here we discuss two possible cases: 1. If there exists t 2 > t 0 such that x(t 2 ) > x 0 , then by the same argument as above, there exists t 3 > t 2 such that , and by the continuity of x, there must exist t 4 such that t 3 > t 4 > t 2 such that x(t 4 ) = x 0 and x′(t 4 ) ≤ 0. Since There exists a constant C such that (66) so for large enough t 0 , we have If t is large enough, by (C − G(x 0 )) > 0, x(t) will be greater than x 0 . This is a contradiction to the assumption that x(t) < x 0 for all t > t 0 .

A.4. Condition for the existence of unique S3
1. When p0 > 0.5 and p̄1 ≤ 0.5, there always exists a unique S3 steady state.
Proof. Using p 0 and p 1 in Hill form (33)- (34) and substituting into (30) and (31) with a and b defined as It is easy to verify that a > 0 and 1 − bδζ > 0. To get positive steady states, we require that , which is equivalent to (71) Therefore, the number of S3 depends on the number of as roots of (70) satisfying (71). Let be the left-hand side of (70) and we consider the following three cases: • If g 1 + bg 2 > 0, then and has only one positive root. Since , (70) has only one positive root which satisfies the inequality (71).

•
If g 1 + bg 2 < 0, then , so has two roots with the same sign. Because , there must be a root larger than the positive number . Therefore, has only one positive root that satisfies the inequality (71).
• If g 1 + bg 2 = 0, then has one positive root. By has only one positive solution satisfying the inequality (71).
The above shows that in any case, (70) only has one positive root satisfying (71), so there is a unique S3 steady state when p̄0 > 0.5 and p̄1 < 0.5. When p1 > 0.5, (72) is equivalent to condition (36). Therefore, there is a unique S3 if and only if (36) holds.

S2 is locally unstable if (36) is true.
Proof. One can easily verify that if (36) is ture, , resulting in a positive eigenvalue of the Jacobian at the S2 steady state , so it is locally unstable.

Moreover, if
by checking (28), one can show that S2 is unstable.

Does not exist Exists
Stability -Section 4.3