Minimal Model for Stem-Cell Differentiation
aa r X i v : . [ q - b i o . CB ] M a r Minimal Model for Stem-Cell Differentiation
Yusuke Goto ∗ and Kunihiko Kaneko † Research Center for Complex Systems Biology,Graduate School of Arts and Sciences,The University of Tokyo, 3-8-1 Komaba,Meguro-ku, Tokyo 153-8902, Japan (Dated:)
Abstract
To explain the differentiation of stem cells in terms of dynamical systems theory, models of inter-acting cells with intracellular protein expression dynamics are analyzed and simulated. Simulationswere carried out for all possible protein expression networks consisting of two genes under cell–cellinteractions mediated by the diffusion of a protein. Networks that show cell differentiation are ex-tracted and two forms of symmetric differentiation based on Turing’s mechanism and asymmetricdifferentiation are identified. In the latter network, the intracellular protein levels show oscillatorydynamics at a single-cell level, while cell-to-cell synchronicity of the oscillation is lost with an in-crease in the number of cells. Differentiation to a fixed-point type behavior follows with a furtherincrease in the number of cells. The cell type with oscillatory dynamics corresponds to a stem cellthat can both proliferate and differentiate, while the latter fixed-point type only proliferates. Thisdifferentiation is analyzed as a saddle-node bifurcation on an invariant circle, while the numberratio of each cell type is shown to be robust against perturbations due to self-consistent deter-mination of the effective bifurcation parameter as a result of the cell–cell interaction. Complexcell differentiation is designed by combing these simple two-gene networks. The generality of thepresent differentiation mechanism, as well as its biological relevance, is discussed. ∗ [email protected] † [email protected] . INTRODUCTION The differentiation of multipotent stem cells into lineage-specific cells is an importantprocess in developmental biology, for which an understanding in terms of dynamical sys-tems theory is desired. Stem cells can both proliferate (i.e., reproduce themselves) anddifferentiate into other cell types [1–4]. While the former indicates that the cellular stateis stable in retaining its composition, the latter implies that the original cellular state isalso unstable in that it moves toward a state with different compositions. How a stem cellsimultaneously supports these two conflicting features is an important question.Dynamical system approaches to cell differentiation have been developed by consideringthat abundances of cellular components are changed through intracellular reactions [5–11].Protein expression dynamics consist of mutual activation and inhibition, and the concentra-tion of each protein changes over time, while its composition determines the cellular state.A dynamical system in the state space that consists of each protein expression level canthen be considered. Based on this picture, it is natural to assign an attractor to each celltype, which allows the robustness of each cell type to be explained as the stability of anattractor against noise. If the system has multiple attractors, each differentiated cell typecorresponds to each attractor.Although this attractor picture is important for understanding the robustness of eachcellular state and how distinct cell types are formed, answers to the following questionsremain elusive: (i) Which attractor describes the two conflicting functions in stem cells, i.e.,proliferation and differentiation? (ii) How are initial conditions for different attracting statesselected through the course of development? (iii) How is the stability of the developmentalcourse, i.e., the robustness in the timing of cell differentiation, and the number distributionof each cell explained, which possibly includes regulation of the cell–cell interactions?To address these questions, coupled dynamical systems that include intracellular dynam-ics of protein abundances, cell–cell interaction, and an increase in the cellular number bycell division were developed [12–16]. Cells with oscillatory intracellular dynamics of proteinconcentrations were shown to differentiate to a novel type with an increase in the cell num-ber, and a cell model whose protein expressions are regulated by a gene regulatory network(GRN) with 5 genes has recently been studied [17]. From extensive simulations, GRNs thatinclude “stem cells” that can undergo both proliferation and differentiation were selected.2he cells always exhibited oscillatory expression dynamics as a single-cell level, and afterdivisions, synchronization of the oscillations among cells was lost. With further cell division,differentiation to a cell type that loses the oscillatory expression followed. Furthermore, it isinteresting to note that recent measurements of protein expression dynamics in embryonicstem cells support this oscillation scenario [18]. The importance of intracellular oscillationto cell differentiation is now being recognized in the developmental biology community [19].To analyze this differentiation mechanism in terms of dynamical systems theory, however,it would be useful to adopt a simpler system consisting of fewer genes. Here, we study sucha minimal system, i.e., dynamical systems consisting of expression levels of only two genes(proteins). After introducing the model in § II, we describe the results in § III. By simulat-ing all possible regulation networks consisting of two proteins with various parameters, weextract a minimal system that allows for stem-cell differentiation. In such a system, theprotein expression levels oscillate in time as a single-cell level. With an increase in the cellnumber, synchronization in the oscillation among cells is lost due to cell–cell interactionsand some cells then differentiate to fall into fixed-point dynamics losing the oscillation. Thisprocess is analyzed using bifurcation theory and explained in terms of a saddle-node on aninvariant circle (SNIC) bifurcation. We also show that the number of each cell type afterdevelopment is robust against noise. The extracted two-gene network leading to the SNICis shown to provide a universal motif for asymmetric differentiation from stem cells, whilecomplex hierarchical differentiation is designed by simply combining the extracted two-genenetwork in parallel or in sequence. The relevance of the present results to dynamical systemsand biological cell differentiation is discussed briefly in § IV.
II. MODELA. Intracellular protein expression dynamics
Here we describe our cell differentiation model under cell-to-cell interactions. Cell statesare represented by protein expression levels of two genes x and y . These x and y genescan regulate the protein expression level of both itself and the other gene. We considerthe dynamics of the expression levels (protein concentrations) of the two proteins, denotedas x i ( t ) and y i ( t ), of the i -th cell at time t . Following earlier studies, we eliminate the3RNA concentration synthesized from each gene and obtain the dynamical system onlyof the protein expression levels. The dynamics of the expression levels at a single cell aredescribed as dx i ( t ) dt = 11 + exp {− β ( J xx x i ( t ) + J xy y i ( t ) − g x ) } − x i ( t ) dy i ( t ) dt = 11 + exp {− β ( J yx x i ( t ) + J yy y i ( t ) − g y ) } − y i ( t ) . (1)The 2 × J mℓ ( m, ℓ = x, y ) gives the transcriptional regulation from protein ℓ to m ,where J mℓ = 1 if the gene ℓ activates the expression of m , − exp ( − β ( z − g )) is adopted [20]to represent the on–off type expression with β as the sensitivity parameter, which roughlycorresponds to the Hill coefficient in terms of cell biology, where for β → ∞ the functionapproaches the step function. The parameters g x and g y give the threshold values of theexpressions of the x and y genes, respectively. B. Cell–cell interaction
The cells interact with each other and this interaction is mediated by some signal. Weassume here that the interaction is mediated directly or indirectly by one of the proteins,which we take to be y . In its simplest form, we consider only the interaction by the diffusionof the y protein. By assuming that the diffusion is fast and by discarding spatial inhomo-geneity over the cells, we employ global coupling, i.e., an all-to-all mean-field interaction.Thus, the overall gene expression dynamics obey the following equation: dx i ( t ) dt = 11 + exp {− β ( J xx x i ( t ) + J xy y i ( t ) − g x ) } − x i ( t ) dy i ( t ) dt = 11 + exp {− β ( J yx x i ( t ) + J yy y i ( t ) − g y ) } − y i ( t )+ DN ( t ) ( N X k =1 y k ( t ) − y i ( t )) , (2)where D is the diffusion constant, and N ( t ) is the total number of cells at time t .4 . Cell division To examine the cell differentiation process through development, the number of cells N ( t )is increased over time by cell division. Starting with a single cell, this cell is divided into twocells that have almost the same protein concentrations at a certain time. The concentrations x i ( t ) and y i ( t ) are slightly perturbed by cell division, leading to cell differences that obey aGaussian distribution with the variance σ , which is set sufficiently small. Here, the behaviorto be discussed (i.e., whether cell differentiation appears) is independent of the noise level.Noise is included mainly to remove synchronization over cells, which is unstable but preserveddue to numerical computation with a finite bit. The cells are simply divided into two afterevery time span t d . D. Model parameters and criterion for differentiation
To numerically investigate this model, we set the sensitivity parameter as β = 40 andthe division noise level as σ = 10 − . The results do not depend on these specific values,as long as the former is sufficiently large (e.g., β >
5) and the latter is not too large (e.g., σ < . t d is chosen to be 25, but this specific value also does not affectthe results.We scanned the other parameters for all possible choices of J mℓ = ± , m, ℓ = x or y )to examine the possibility of cell differentiation. We also checked all possible configurationsof J mℓ as long as the matrix includes at least one off-diagonal component (i.e., as a mini-mum, the interaction between x and y exists). This gives a total of 36 configuration. Theparameters g x,y were varied from − D was chosen as either0.2 or 1.0. Thus, the model in Eq.(2) was simulated for a total of 36 × × × N = 32. Foreach cell, we computed the temporal average of the gene expression level for a sufficientlylong time after discarding the transient time. If the difference between the maximum andminimum averaged gene expression levels ( x i ( t ), y i ( t )) over the cells was larger than 0.1,then we concluded that differentiation had occurred.5 II. RESULTSA. Differentiation classification
Type Behavior Network ParameterTuring fixed point → fixed point 2/36 6/121032Oscillation death oscillation → fixed point 2/36 3/121032Asymmetric differentiation oscillation → with oscillation oscillation 2/36 9/121032The simulation results revealed cell differentiation cases that can be classified into thefollowing three types: 1) Turing, 2) Oscillation death, and 3) Asymmetric differentiationwith remnant oscillations. In all cases, the single-cell dynamics had only a single attractor.As the cell number increases, the cells take two distinct states. We first describe these threetypes and show that the mechanisms of the first two types are already well known. We thusfocus on the third type, which is the most relevant to differentiation of stem cells.
1. Turing type (fixed point → fixed point) Fig. 1: (a) Upper: Network structure of a model that shows Turing-type differentiation; fixed point → fixedpoint. Lower: The conjugate network form. The parameter values that show this differentiation are ( D, g x , g y ) =(0 . , , , (1 , , , (1 , . , . , (1 , − . , − .
05) for the upper network and (
D, g x , g y ) = (0 . , , − , (1 , , − , (1 , . , − . x = y in themodel here and is not shown). (b) The single-cell nullcline of the model. (c) Time series of the protein expression x i ( t ) as thecell number increases from i = 1 to 32 through cell division every 25 time units; (b) and (c) are drawn for the upper networkstructure with the parameter values ( D, g x , g y ) = (0 . , , The single-cell dynamics in this case has a unique stable fixed point. As the cell numberincreases, the cell population splits into two groups, each taking different fixed-point values6f x and y and a higher value of x or y . The numbers in each population are equal, i.e., 16cells each for N = 32. One gene expression network type is shown in Fig. 1, where one geneactivates the expression of itself and the other, while the other gene, which is diffusible,inhibits the expression of itself and the other. This mechanism is explained well by theclassic Turing pattern [21] by replacing the spatially local diffusive interaction therein withglobal coupling. Indeed, Turing’s seminal paper included the present case, as he discussedthe case with N = 2 (see also [22]). The conjugate network is also shown in Fig. 1 (inparentheses), where the differentiation mechanism is understood in the same way.
2. Oscillation death (oscillation → fixed point) The single-cell dynamics for this type of differentiation has a limit-cycle attractor. Withan increase in cell number, the cell population again splits into two groups of equal number,each of which shows distinct fixed points through the cell-to-cell interaction, in the samemanner as the Turing-pattern case. This loss of oscillation is known as oscillation death[21, 23–25]. Two forms of the networks that show this behavior are shown in Fig. 2. Themechanism of this type of differentiation will be explained later, as it is common to that ofasymmetric differentiation in the next section. x yx y ( )
A) B) C)
Fig. 2: (a) Upper: Network structure of a model that shows differentiation through oscillation death; oscillation → fixed point. Lower: The conjugate network form. The parameter values that show this differentiation are ( D, g x , g y ) =(1 , − . , . , (1 , . , D, g x , g y ) = (1 , . , − .
15) for the lower network. (The flow in the statespace and the nullcline is simply a mirror symmetry with regard to x = y in the model here and is not shown). (b) Thesingle-cell nullcline of the model. (c) Time series of the protein expression x i ( t ) as the cell number increases from i = 1 to 32through cell division every 25 time units; (b) and (c) are drawn for the upper network structure with the parameter values( D, g x , g y ) = (1 , − . , . . Asymmetric differentiation with remnant oscillation (oscillation → oscillation transition) In this class, the attractor of a single cell is again a limit cycle, but as the number of cellsis increased, synchronized oscillation over cells is unstable under the cell–cell interaction,and protein-expression oscillations are desynchronized between some cells [28]. After thisdesynchronization, some cells leave the original limit-cycle trajectory and enter an oscillatorystate with a tiny amplitude, while the other cells remain in the original limit cycle (with aslight modification) [30]. The states mutually stabilize their existence, e.g., if all the cellsin one state were removed, some of the remaining cells will make a transition to the otherstate. Thus, two distinct oscillatory states coexist independently of the initial condition ofthe cells.This class of behavior is observed for the two networks, shown in Fig. 3, over 9 setsof parameter values. The parameter interval in which the present behavior is observed isrelatively small compared with the previous two cases.
Fig. 3: (a) Upper: Network structure of a model that shows asymmetric differentiation with remnant oscillation; os-cillation → oscillation transition. Lower: The conjugate network form. The parameter values that show this dif-ferentiation are ( D, g x , g y ) = (0 . , − . , . , (0 . , − . , . , (0 . , − . , − .
25) for the upper network and (
D, g x , g y ) =(0 . , . , . , (0 . , . , . , (0 . , . , .
85) for the lower network. (The flow in the state space and the nullcline is simplya mirror symmetry with regard to x = y in the model here and is not shown). (b) The single-cell nullcline of the model. (c)Time series of the protein expression x i ( t ) as the cell number increases from i = 1 to 32 through cell division every 25 timeunits; (b) and (c) are drawn for the upper network structure with the parameter values ( D, g x , g y ) = (0 . , − . , . B. Mechanism of oscillatory differentiation
We note that none of the cells in the first two classes have the potential to both pro-liferate (reproduce the same type by division) and differentiate (switch to a cell type withdistinct behavior). Furthermore, the first two mechanisms are already well understood based8n Turing’s study, while the latter mechanism, which corresponds to the cell-differentiationmechanism in our earlier studies, is not fully understood in terms of dynamical systems the-ory. Hence, we focus here on the third class and explain the mechanism of the differentiationinto two states.The differentiation is based on two stages: desynchronization of the oscillations and SNICbifurcation induced by cell-to-cell interactions [26, 27]. Examining the nullcline of the single-cell dynamics, shown in Fig. 4, we see that while the two nullclines come close, they do notintersect. Thus, there are no fixed points, and the limit-cycle attractor exists as a single-celllevel, as already mentioned. (For the following numerical simulations, we have used theparameter values g x = − . g y = 0 .
2, and D = 0 . Fig. 4: The nullcline and limit-cycle attractor of the single-cell dynamics for the GRN in Fig. 3. The network is redrawn forreference, where the parameter values are J xx = 1 , J xy = − , J yx = 1 , J yy = 0 , g x = − . , g y = 0 .
2. The two nullclines for dx/dt = 0 and dy/dt = 0 are drawn in red and blue, respectively. The limit-cycle attractor is drawn as the dashed line. The twonullclines come close but do not intersect, which is an important property for differentiation to occur by cell–cell interaction.Indeed, the nullclines and the limit cycle are arranged in the same manner for the conjugate network (shown in parentheses ofFig. 3a). y x Fig. 5: Dynamics of the model in Figs. 3 and 4 with D = 0 .
14. (a) Locus of the orbit of the two-cell dynamics of ( x ( t ) , y ( t ))and ( x ( t ) , y ( t )). Blue: cell type whose state is close to the original limit cycle. Red: differentiated type with loss ofautonomous oscillation. (b) The time series of x i ( t ) for the two-cell dynamics ( i = 1 , t ∼
80, they split into two distinct steady states.
1. Desynchronization
Consider a coupled system with cell-to-cell interaction, where each cell shows a limit-cycleoscillation as shown above. As is known, the oscillations of the cells lose synchronization fora certain range of the interaction and parameter values. This desynchronization occurs inthe present case, as was numerically confirmed by computing the stability exponent for thesynchronization. Indeed, the magnitude of the tangential vector representing the deviationbetween the two cells δ X =( δx, δy ) increases exponentially over time.
2. Differentiation
After desynchronization, the cell-to-cell interaction leads to differentiation with an in-crease in the cell number. Here, the cell–cell interaction term that occurs through diffusionin our model is represented by α ( t ) ≡ D ( N P Nk =1 y k − y i ) = D ( y − y i ), so that the dynamicsare written as dx ( t ) dt = 11 + exp {− β ( x ( t ) − y ( t ) − g x ) } − x ( t ) dy ( t ) dt = 11 + exp {− β ( x ( t ) − g y ) } − y ( t ) + α. (3)If the oscillations are synchronized over cells, then α ≡
0, and the equation is reduced tosingle-cell dynamics. Due to the desynchronization, however, α is nonzero and functions asa time-dependent bifurcation parameter. We first examine the bifurcation of the single-cell10ynamics against changes in the constant parameter α and then discuss the time dependenceof α .The nullclines of Eq. (2) are shown in Fig. 6. As α is changed, only the nullcline dy/dt = 0moves vertically downward with a decrease in α without a change in form; there is no effecton the nullcline dx/dt = 0. The two nullclines intersect with a slight decrease in α , so thatthe SNIC [26, 27] appears, and the limit-cycle attractor is replaced by a stable fixed point.(In the example in Fig. 6, this bifurcation occurs at α c ≈ − . y ,i.e., the average of y over α , is represented as a constant value that can deviate from y i ( t ),the nullcline for dy/dt = 0 is represented by y = y c = 1 / { (1+ exp {− β ( x − g y ) } )(1+ D ) }− Dy and it decreases if D increases. Recalling that the x nullcline y = x + β log (1 − /x ) − g x ≈ x = 1(see figure), x and y nullcline intersect and a fixed point appearsvia the SNIC, if D is large enough. In the differentiation of two cells that originates indesynchronization, as shown in Fig. 5, one cell has a larger y i value, and the α ( <
0) valueof that cell is smaller and can be smaller than α c . The nullcline dy/dt = 0 for such a cellthen intersects the nullcline dx/dt = 0. In contrast, the other cell remains close to theoriginal limit cycle. The effective α , estimated from the diffusion term, shown in Fig. 7,demonstrates that α remains in the region of the fixed point( α < α c ), which supports thisargument. Thus, for a certain parameter region, the two types of behavior are differentiated:one is close to the original limit cycle and the other is close to the fixed point generated bythe SNIC. y y dx/dt=0dy/dt=0 ( Α = 0)dy/dt=0 ( Α = - 0.017) Fig. 6: (a) Change in the nullcline of dy/dt = 0 as α in Eq. (3) changes. As α is decreased, the nullcline dy/dt = 0 movesdownwards (from the red to green line) and touches the nullcline of dx/dt = 0 (blue), when α = − .
20 40 60 80 100 - - Α H t L Α c = - 0.017 Α Α Fig. 7: Time series of the diffusion terms D ( y ( t ) − y ( t )) (blue) and D ( y ( t ) − y ( t )) (red) of the two cells. The cell–cellinteraction induces desynchronization, and the difference between the cells is amplified. One cell then stays in the region ofthe fixed point. For reference, the value α c = − .
017 is plotted. Note that D ( y ( t ) − y ( t )) remains below α c for most of thetime, while D ( y ( t ) − y ( t )) remains above it. Since the variables x and y oscillate in time, this argument assuming a constant α isinsufficient. Still, on average, the trajectory of the differentiated cells that take a larger y value remains above the nullcline dy/dt = 0 in the state space ( x, y ). Although the trajectorystays around the intersection of the two nullclines, the state of this cell is not completely afixed point due to the interaction term with the other cell: since the other cell type keeps theoscillatory dynamics close to the original limit cycle, the differentiated cell type is driven bythis oscillatory dynamics so that it shows an oscillation with a tiny amplitude. The originaland differentiated cell types are distinguished by the ability for autonomous oscillation. Notethat oscillation death(case2) is understood similary as the asymmetric differentiation. Inthats case, SNIC bifurcation occur at two points around ( x, y ) = (0 ,
0) and (1 ,
1) on x - y plane, and the limit cycle is replaced by the two fixed points as a result.
3. Parameter dependence
Since the desynchronization and differentiation are due to the interaction, the dynamicsof the two cells depends strongly on the diffusion coefficient D . By taking a two-cell systemwith a fixed g x,y as above, we can examine the dependence of the two-cell dynamics on thevalue of D . In Fig. 8, the local maxima x i ( t ) of the two cells i = 1 , D , i.e.,12or weak cell–cell interaction, the oscillations of the two cells are synchronized. With anincrease in D , the synchronization loses stability, and the oscillations desynchronize (Fig.8); coexistence of the two oscillations occurs for larger D .The corresponding changes in the time series and trajectories in ( x, y ) space are shown foreach region (A1, A2, ..., A5) in Fig. 9. From A1 to A2, a period-doubling bifurcation causesthe oscillations of the two cells to be out of phase. From A2 to A3 a pitchfork bifurcationoccurs. Then, the period-doubling cascade to chaos appears from A3 to A4, where the twooscillations are chaotic with desynchronization. At the transition from A4 to A5, chaoticorbit touches with saddle, and crisis occurs. The chaotic orbit in the four dimensional spacewith two cells becomes unstable and is replaced by two stable periodic orbits. Each of theoscillations is now periodic; one with a large amplitude close to the original limit cycle andthe other with a tiny amplitude near the fixed point. Form the view point of single celldynamics, it can be regarded as occurrence of SNIC bifurcation, when we regard diffusionterm as time dependent parameter. In region A5, the asymmetric differentiation from thestem-type cell, which was discussed earlier, occurs. We also note that the bifurcation in thecase of oscillation death follows almost the same sequence as the present case. A1 A2 A3 A4 A5 l o ca l m a x i m u m o f x A2 A3 A4 l o ca l m a x i m u m o f x Fig. 8: Bifurcation diagram of the two-cell dynamics. (a) Local maxima of x i ( t ) for two cells ( i = 1 ,
2; red and blue)plotted at the steady state (attractor). (b) Details around the A3-A4 transition. As D is increased, the single-cell orbit isdestabilized ( A → A D , distinct cells with different oscillations coexist (A5). y
100 120 140 160 180 20000.51 time x y
100 120 140 160 180 20000.51 time x y
100 120 140 160 180 20000.51 time x y
100 120 140 160 180 20000.51 time x y x Fig. 9: Orbit and time series of the two cells. Left: Orbits of ( x i ( t ) , y i ( t )) ( i = 1 ,
2) after the cells reach a steady state(attractor). Right: The time series of x i ( t ). The colors correspond to the two different cells i = 1 ,
2. (i) D = 0 . D = 0 .
12 corresponding to A2. The orbits are desynchronized.(iii) D = 0 . D = 0 . D = 0 .
132 corresponding to A5. After desynchronization, the dynamics of the two cells are split into two distinctbehaviors. The differentiation discussed in section III.B.2 corresponds to this region, while D adopted in that discussion was0.14. . Robustness in cell number regulation N A (cid:144) N P H N A (cid:144) N L Fig. 10: Histogram of the ratio of the number of original type- A cells with a large-amplitude limit cycle to the total numberof cells. Here, N = 100 cells are used with 200 sets of randomly chosen initial values for ( x i , y i ) distributed homogeneouslybetween 0 and 1. After the expressions reached the steady state, the ratio is computed to construct the histogram. N A H t = L N A H t = L Fig. 11: The range of differentiated cell types obtained from an initial condition of two types of cells. The abscissa is theinitial number of type A cells (the original cell type with a large-amplitude limit cycle) among 100 cells, and the ordinate is thenumber of such cells after reaching the steady state at t = 100. For each initial number of type A cells, 200 initial conditionsare chosen for the ( x i , y i ) cells. The dark blue line shows the average of the ratio of type A cells over the 200 samples, whilethe blue bars show the range of the standard deviation from the distribution. We have identified three classes of cell differentiation in terms of dynamical systemstheory. The first class is that in which multiple attractors exist, where the switch amongthe attractors by noise leads to cell differentiation. The second is a Turing-type class in15hich the initial cell state is unstable and leads to two different states. The third class isasymmetric differentiation from an oscillatory state. Only the third class has a type of cellthat can both proliferate and differentiate, which reflects the nature of stem cells.In the multiple-attractor case, the differentiation is due to noise so that the time courseof the development is stochastic; the number ratio of each cell type is unregulated. In theTuring-type case, the original cell type just disappears, so that there is no stem-cell-like cell.One merit of asymmetric differentiation lies in both the existence of stem-cell-type cellsand the robustness of the ratio of each cell-type number against noise [16, 17] or a changein the initial conditions. This robustness is expected, as the differentiated cell type appearsas a result of the instability of the state consisting of the first cell type only, and the twocell types stabilize each other through cell–cell interactions. In this section, we call theoriginal cells with a large-amplitude limit cycle as type-A cells, and differentiated cells witha small-amplitude limit cycle as type-B cells and study the robustness in the number ratioof type A and B cells. To examine this robustness, we first simulated our model startingwith a single cell for a given initial condition and progressively added noise. We confirmedthat the ratio of one type of cell to the total number of cells has a narrow distribution evenin the presence of noise, once the total number of cells reached a particular number (e.g.,32).We then started with N initial cells with random initial values of ( x i , y i ) distributedhomogeneously between 0 and 1. As shown in Fig. 10, the ratio has a narrow distributioncentered around 0.43. To check the possible range of numbers of the two cell types, wesimulated the model with initial conditions of the expression levels so that there were N A type- A cells and N − N A differentiated cells. If N A is initially too large, then the stateis unstable, and some of the N A cells show a SNIC bifurcation to lose the oscillation anddifferentiate to type- B cells with fixed-point behavior. If N A is initially too small, (i.e., theproportion of differentiated cells of type N B is too large), then the type- B state is unstableand some of these type- B cells regain the oscillation to de-differentiate to type- A cells. Thus,there are upper and lower bounds on the ratio of the number of type- A cells. For example,in the parameter values used in Fig. 11, this range is 0.2 to 0.55. Therefore, the robustnessof the two types of cells is a consequence of the present asymmetric differentiation based onthe SNIC. 16 . Complex cell differentiation by combining two-gene motifs Fig. 12: Basic cell lineage diagram. Complex differentiation is composed of a combination of these two processes. (a) Parallelcase: two types of cells are created from one type of cell. (b) Sequential case: cells differentiate in series.Fig. 13: The parallel case. (a) The GRN. (b) The dynamics of the three cell types are plotted as orbits in x, y, z space. Whenasymmetric differentiation occurs in the x – y network, the protein z is expressed, which leads to a Turing instability. As aresult, cells with small-amplitude oscillations in the x – y plane split into two states with regard to the value of z (type- B and- C cells).The parameter values that show this differentiation are( D y , D w , g x , g y , g z , g w ) = (0 . , . , − . , . , . , x z Fig. 14: The sequential case. (a) The GRN. (b) The dynamics of the three cell types are plotted as orbits in x, y, u space.Asymmetric differentiation occurs in the x – y network and some cells change from type- A to type- B cells. The type- B cellsshow asymmetric differentiation in the u – v gene network, which results in some type- B cells changing into type- C cells. Theparameter values that show this differentiation are ( D y , D w , D v )= (0 . , . , . g x , g y , g z , g w , g u , g v ) =(-0.1,0.2,0.99,0.0,-0.07,-0.2). J mℓ = 1 if the gene ℓ activates the expression of m , − J uz = − . Most of the present cells consist of a large number of cell types that are generated throughsuccessive differentiations. The entire cell lineage can be constructed by combining thefollowing two differentiation processes.(i) The parallel case: bifurcation into two types of cells from one type of cell, as shownin Fig. 12(a).(ii) The sequential case: cells differentiate in series to form a hierarchical differentiation,as shown in Fig. 12(b).Complex cell differentiations as observed in a hematopoietic system, for example, areshaped by combining these parallel and sequential differentiations. We show here that thesetwo basic forms are designed straightforwardly by combining the Turing-type module andthe asymmetric differentiation module.The parallel case can be achieved by a gene network as shown in Fig. 13. The Turing-type module is regulated by the asymmetric differentiation module. There are four proteins x, y, z, and w in total whose concentrations are represented by the corresponding variables.The protein levels of the cells initially show oscillations, forming a large-amplitude limit cyclein the x – y plane (type- A cells). As the number of cells increases, asymmetric differentiationoccurs in the x – y plane to form a state with constant, higher expressions of x and y viathe mechanism mentioned earlier. With this activation of x , a Turing-type bifurcation in18he z and w protein expressions is induced. Both z and w are expressed in the beginning,but following the suppression of z , the state differentiates into two types of cells with higherand lower expressions of both z and w based on the Turing-type mechanism. Hence, thedifferentiation from stem-cell-type A to two types of cells B and C occurs, as shown in Fig.13.The sequential case, however, involves the gene network shown in Fig. 14 in whichtwo asymmetric differentiation motifs (( x, y ) and ( u, v )) are connected in sequence withthe intervening Turing motif ( z, w ). Here again, all the cells initially show oscillation andform a large-amplitude limit cycle in the ( x, y ) plane (type- A cells). As the number ofcells increases, asymmetric differentiations occur and some cells differentiate into type- B cells, which fall on a small-amplitude limit cycle in the ( x, y ) plane. After this bifurcation,the expression of x is constantly activated in the type- B cells, which then suppresses theexpression of z and activate the expression of u , triggering the asymmetric differentiation ofthe expressions of ( u, v ) and leading to the differentiation to type- C cells with a constantactivation of u . Thus, hierarchical, sequential differentiation from type- A to B and then to C is generated, through which the oscillation amplitude decreases accordingly. (Here, theintermediate Turing module is used to suppress the expression level of z when the expressionof x is constantly activated and to activate the expression otherwise. This is not an essentialcomponent, as we believe other network forms can be adopted.) IV. DISCUSSION
We have studied in this paper an interacting cell model consisting of two genes (proteinexpressions) and extracted, from extensive simulation, minimal gene networks that showdifferentiations. These differentiations were classified as Turing, oscillation death, or asym-metric differentiation with remnant oscillations (see also [33, 34] for other possible types fordifferentiation including an inhomogeneous limit cycle in which two limit cycles coexist).Only asymmetric differentiation was shown to allow for cells with stemness, i.e., compatibil-ity with both proliferation and differentiation. The differentiation is understood as a SNICbifurcation which is triggered by cell–cell interactions with the cell number of the originalcell type being the bifurcation parameter. It was shown that each cellular state is stabilizedaccording to the cell–cell interaction, which depends on the number distribution of each19ell type. The number ratio is thus regulated autonomously. In this sense, the effectivebifurcation parameter due to cell–cell interactions is “self-consistently” determined [31]. Atheoretical analysis for such a self-consistent bifurcation should be developed in future.With the present study, we can now answer the questions addressed in the Introduction.(i) Which attractor describes the two conflicting functions in stem cells, i.e., proliferation anddifferentiation? – A single-cell attractor providing stemness is a limit-cycle that is close tothe point of SNIC bifurcation. The limit-cycle attractor provides stability, while the cellularstate is easily switched to a different state with the aid of SNIC bifurcation, due to cell-cellinteraction. (ii) How are initial conditions for different attracting states selected throughthe course of development? – With the cell-cell interaction, desynchronization of oscillationsfollows, which diversifies the cellular states. Then, with the cell-cell interaction, states ofsome cells are kicked out from the original basin of attraction, triggered by SNIC bifurcation.(iii) How is the stability of the developmental course explained, which possibly includesregulation of the cell–cell interactions? – Since the desnchronization and bifurcation occuras a result of the increase in cell number, the differentiation timing is almost deterministicthrough the developmental course. The cell types thus generated are stabilized with eachother through cell-cell interaction, which depends on the number ratio of each cell type.Hence the ratio of each cell-type is regulated so that it stays within a certain proportion,and it is robust to noise.We were able to extract the minimal gene expression network for the stemness and demon-strated that it consists of two genes; one activates and suppresses the other, while the otherhas an activation path to itself. Due to the simplicity, the network can work as a motif[32] for complex cell differentiation in general. In fact, several networks previously studiedfor three or more genes [17] include the present minimal network motif as their core com-ponent. By including feedback or feed-forward path(s) for gene expression networks to theextracted minimal network, stem-cell-type behavior as well as the differentiation processfurther enhances the robustness against change in the parameter values.Although we have adopted a simple form of the threshold expression dynamics, the bi-furcation analysis developed here shows that the present mechanism for differentiation ispossible in other forms of the expression dynamics, as long as the SNIC occurs due to thecell–cell interaction or signal molecule. For example, by using the Hill form for the expres-20ion, ǫ dx i ( t ) dt = ( x i ( t ) K x ) n x i ( t ) K x ) n
11 + ( y i ( t ) K y ) n − x i ( t ) + I dy i ( t ) dt = ( x i ( t ) K x ) n x i ( t ) K x ) n − y i ( t ) + I + D (cid:16) N X k =1 y k ( t ) − y i ( t ) (cid:17) , (4)the present differentiation progresses for appropriate values of the parameters (with Hillcoefficients n i ( i = 1 , ,
3) of sufficiently large values).Real GRNs, however, include more than a thousand genes, and actual networks are quitecomplicated. In spite of this, the present network motif, as it is so small, can be easilyincluded in the networks of the present cell. A combination of the present two-gene networkmotifs can lead to differentiations of a complex cell lineage, as was demonstrated here. Itwill be important to extract such network motif combinations in relation to the observedcomplex differentiation [35].Acknowledgement:The authors would like to thank Shuji Ishihara, Chikara Furusawa, Benjamin Pfeuty,and Narito Suzuki for useful discussions. This work was partially supported by a Grant-in-Aid for Scientific Research (No. 21120004) on Innovative Areas “Neural creativity forcommunication” (No. 4103) and the Platform for Dynamic Approaches to Living Systemfrom MEXT, Japan. [1] Lanza R, Gearhart J, Hogan B, Melton D, Pedersen R, Thomas ED, Thomson J, West M(2009)
Essentials of Stem Cell Biology . 2nd edition, Academic Press, San Diego, USA[2] Potten CS, Loeffer M (1990) Development 110: 1001-1020[3] Ramalho-Santos M, Yoon S, Matsuzaki Y, Mulligan RC, Melton DA (2002) Science 298: 597-600[4] Slack JM (2002) Nat. Rev. Genet. 3: 889-895[5] C. H. Waddington (1957)
The Strategy of the Genes . George Allen & Unwin, London[6] Forgacs G, Newman SA (2006)
Biological Physics of The Developing Embryo . CambridgeUniv. Press, Cambridge, UK
7] Goodwin BC (1963)
Temporal Organizations in Cells . Academic Press, San Diego, USA[8] Kauffman SA (1993)
The Origins of Order: Self-Organization and Selection in Evolution .Oxford University Press, Oxford, UK[9] Kauffman SA (1969) J. Theor. Biol. 22: 437-467[10] Glass L, Kauffman SA (1973) J. Theor Biol. 39: 103-129[11] Chang HH, Hemberg M, Barahona M, Ingber DE, Huang S (2008) Nature 453: 544-548[12] Kaneko K, Yomo T (1994) Physica D 75: 89-102[13] Kaneko K, Yomo T (1997) Bull. Math. Biol. 59: 139-196[14] Kaneko K, Yomo T (1999) J. Theor. Biol. 199: 243-256[15] Furusawa C, Kaneko K (1998) Bull. Math. Biol. 60: 659-687[16] Furusawa C, Kaneko K (2001) J. Theor. Biol. 209: 395-416[17] N. Suzuki, C. Furusawa, K. Kaneko (2011), PLoS One 6, e27232.[18] Kobayashi T, Mizuno H, Imayoshi I, Furusawa C, Shirahige K, Kageyama R (2009) GenesDev. 23(16): 1870-1875[19] Furusawa C, Kaneko K (2012) Science 338: 215-217[20] Mjolsness E, Sharp DH, Reinitz J (1991) J. Theor. Biol. 152: 429-453[21] Turing A M (1952) Phil. Trans. Roy. Soc. B 237: 37-72[22] Mizuguchi T, Sano M (1995) Phys. Rev. Lett. 75: 966-969[23] Prigogine I, Lefever R (1968) J. Chem. Phys. 48: 1695-1700[24] Bar-Eli K (1985) Physica D 14: 242-252; Aronson D, Ermentrout GB, Kopell N (1990) PhysicaD 41: 403-449[25] Koseska A, Volkov E, Kurths J (2010) Chaos 20: 023132[26] Izhikevich EM (2006)
Dynamical systems in neuroscience: the geometry of excitability andbursting . MIT Press, Cambridge, MA[27] Ermentrout GB, Kopell N (1986) SIAM J. Appl. Math. 46: 233-253[28] Kaneko K (1990) Physica D 41: 137-172[29] Okuda K (1993) Physica D 63: 424-436[30] Pattern formation based on this type of oscillation dynamics is discussed in Cooke J, ZeemanEC (1976) J. Theor. Biol. 58: 455-476; Gierer A, Meinhardt H (1972) Kybernetik 12: 30-39[31] Nakajima A, Kaneko K (2008) J. Theor. Biol. 253: 779-787[32] Milo R, Shen-Orr S, Itzkovitz S, Kashtan N, Chklovskii D, Alon U (2002) Science 298: 824-827
33] Ullner E, Zaikin A, Volkov EI, Garc´ıa-Ojalvo J (2007) Phys. Rev. Lett. 99: 148103[34] Koseska A, Ullner E, Volkov EI, Kurths J, Garc´ıa-Ojalvo J, (2010), J. Theor. Biol. 263: 189-202[35] Morrison SJ, Shah NM, Anderson DJ (1997) Cell 88: 287-298; Loh YH, et al. (2006) Nat.Genet. 38: 431-44033] Ullner E, Zaikin A, Volkov EI, Garc´ıa-Ojalvo J (2007) Phys. Rev. Lett. 99: 148103[34] Koseska A, Ullner E, Volkov EI, Kurths J, Garc´ıa-Ojalvo J, (2010), J. Theor. Biol. 263: 189-202[35] Morrison SJ, Shah NM, Anderson DJ (1997) Cell 88: 287-298; Loh YH, et al. (2006) Nat.Genet. 38: 431-440