Symbiotic Cell Differentiation and Cooperative Growth for the Emergence of Multicellularity

The origin of multicellularity is a fundamental open question in biology. For multicellular organisms to evolve from an aggregate of unicellular organisms, cells with an identical genotype must first differentiate into several types. Second, this aggregate of distinct cell types should show better growth than that of an isolated cell in the environment. Third, this cell aggregate should show robustness in the number distribution of differentiated cell types. To reveal how an ensemble of primitive cells achieves these conditions, we developed a dynamical-systems model of cells consisting of chemical components with intracellular catalytic reaction dynamics. The reactions convert external nutrients to internal components for cellular growth, and the divided cells interact through chemical diffusion. We found that cells sharing an identical catalytic network spontaneously differentiate induced by cell-cell interactions, and then achieve cooperative division of labor, the mutual use of products among differentiated cell types, enabling a higher growth rate than that in the unicellular case. This symbiotic differentiation emerged for a class of reaction networks under the condition of nutrient limitation and strong cell-cell interactions. Then, robustness in the cell type distribution was achieved, while instability of collective growth sometimes emerged even among the cooperative cells when the internal reserves of chemical products is dominant. The simplicity and generality of the present mechanism suggests that evolution to multicellularity is a natural consequence of interacting cells with limited resources, being consistent with the behaviors and forms of several extant primitive forms of multicellularity, such as certain bacteria.


INTRODUCTION
The origin of multi-cellular organisms (MCOs) is considered one of major transitions in evolution [1], and thus the mechanisms for their emergence from simple aggregates of reproducing cells have gathered much attention from both theorists and experimentalists over the last few decades [2][3][4][5][6][7]. Indeed, multicellularity is known to have evolved at least 25 times independently [8], and primitive forms of MCOs are found in bacterial biofilms and slime molds, in which the fundamental properties of multicellularity, i.e., cell differentiation, division of labor, or apoptosis, are consistently observed [9]. Such primitive forms of MCOs have also been found to emerge in the course of experimental evolutions of unicellular yeast [10] and algae [11]. These observations suggest that the transition to an MCO is a universal route from unicellular organisms, as the number of cells increases to form a crowded aggregate. There are at least three requisites for the emergence of a primitive form of robust MCO from an aggregate of cells: (i) Cell differentiation: starting from a single cell, cell states become differentiated as their number increases, and they coexist and grow. (ii) Cooperative growth: these different cell types must cooperate for their growth, possibly by division of labor, in order to achieve higher survival than the isolated unicellular case. (iii) Robustness in the distribution of the number of differentiated cells: balance in the population distribution of cell types must be maintained so that cells of different types can coexist stably and grow together. Indeed, the first property, differentiation from a single cell type, has been studied using a dynamical system of interacting cells, which led to proposal of the concept of isologous diversification [12]. However, such differentiation capacity does not necessarily guarantee an evolutionary advantage of an MCO. For an MCO to have greater fitness in a certain condition than a unicellular organism, growth rate as a cell ensemble should be higher than that of an isolated cell. A previous mathematical model [7] demonstrated that an ensemble of cells sharing a common genotype could achieve niche differentiation through cell differentiation, and thereby relax the strength of resource competition. However, in this case, the differentiated cells do not help each other to grow faster than isolated cells, and thus the requisite (ii) is not satisfied. For a benefit of multicellularity over isolated cells (unicellular organisms), the differentiated cells have to cooperate for their growth. Simply achieving cooperative growth, however, does not necessarily imply that this state is robust, since if one cell type reproduces faster than any other type, the fastest type would dominate the population and cell type diversity would be easily lost. Therefore, coexistence of diverse cell types has to be sustained for MCOs to emerge, leading to condition (iii). A related question has been discussed in multi-level evolution theory, by introducing a fitness parameter at the cell level and at the MCO level and investigating how these two fitness parameters are aligned [6,13,14]. This is an important question to address, since conflict between these two levels might be an obstacle for the emergence of an MCO. However, studies based on prescribed fitness values are not sufficient to elucidate the mechanism. In nature, such fitness is determined by the cellular growth rate, which depends on two main factors: the available resources and the strength of cell-cell interactions. Therefore, fitness is not predetermined, but rather changes according to the distribution of cells of different types. Thus, it is natural to consider cells with intracellular reaction dynamics, cell-cell interactions, and uptake of resources, by which the fitness is determined as the cellular growth rate under these conditions rather than being prescribed in advance. In the present paper, by using a simple model of cells that contain diverse components and interact with each other through the exchange of chemicals, we address the question of whether multicellularity is a necessary outcome of an ensemble of interacting cells. Specifically, we show that a cell ensemble under strong cell-cell interactions with limited resources fulfills cell differentiation, cooperative growth, and robustness in the distribution of cell types, and thus a prototype of robust MCOs is achieved.

MODEL
We consider a mathematical model proposed in [7,12,15], which describes a simple, primitive cell that consists of k chemical components {X 0 , . . . , X k−1 }. As illustrated in Fig. 1, we assume that n cells globally interact with each other in a well-mixed medium. The internal state of each cell is characterized by a set of variables (x is the concentration of the i-th chemical X i and v (m) is the volume of the mth cell (m = 1, . . . , n). As a simple model, we consider a situation with only catalysts and resources, where these k components are mutually catalyzed for their synthesis, thus forming a catalytic reaction network. A catalytic reaction from a substrate X i to a product X j by a catalyst X l , as X i + αX l →X j + αX l , occurs at a rate ǫx where α refers to the order of the catalytic reaction and is mostly set as α = 2. Here, ǫ is the rate constant for this reaction, and, for simplicity, all the rate constants are equally fixed at ǫ = 1.
Cell states change through intracellular biochemical reaction dynamics and the in-and outflow of chemi- where P (i, j, l) takes the value 1 if there is a reaction X i + αX l →X j + αX l , and is 0 otherwise. In Eq.
(1), the third term describes the influx of X i from the medium, and the fourth term gives the dilution owing to the volume growth of the cell, so that the constraint of Here, only a subset of chemical species is diffusible across the cell membranes with the rate of diffusion D. X i is transported from the medium into the m-th cell at a rate Dσ i (x and is 0 otherwise. Therefore, the cells interact with each other through the diffusible chemicals. In addition, the m-th cell grows in volume according to the rate ), and its volume dynamics are given by dv (m) /dt = µ (m) v (m) . Besides catalytic chemicals, there are also nutrient chemicals that cells require for growth. In our simulation, only X 0 is defined as the nutrient. The nutrient chemical X 0 is supplied into the medium from the external environment according to the rate D med (C − x concentration of X 0 . Therefore, the temporal change of x (med) i is given by where σ ′ 0 takes unity only if i = 0, i.e., if X i is the nutrient. For simplicityD med was set as D med = D, though the results reported here do not greatly depend on the value of D med . According to these processes, each cell grows until its volume doubles, and then the cell divides into two cells with almost the same chemical compositions. Here, the catalytic network in daughter cells is identical to that in their mother cell. As the initial condition, only a single cell exists with a randomly determined chemical composition. In addition, we set the carrying capacity of a medium N , which is an upper limit to the number of cells that can coexist in the medium. When the cell number exceeds its upper limit N due to cell division, the surplus cells are randomly eliminated. Hereafter, this model is referred to as the N -cell model.

Cell differentiation in the N-cell model: brief summary
We simulated the N -cell model over hundreds of randomly generated reaction networks (see Supplemental Material for details).
We were particularly interested in if and how the cells differentiate, and whether the growth rate would increase as a result of differentiation. For this purpose, cell differentiation is defined as the emergence of cells with different chemical compositions within the population that share an identical catalytic network. For the case where the concentrations asynchronously oscillate in time, we evaluated whether cells have different compositions even after taking the temporal average over a sufficiently longer time scale than the oscillation period. To evaluate the growth enhancement, we compared two different situations, "interacting" (N = 100) and "isolated" (N = 1) cases, and then we computed R µ , the ratio of the growth rate µ of interacting cells to that of isolated cells. The growth enhancement is defined as R µ > 1. The behavior of the N -cell model is classified into four categories. In category (a), interacting cells differentiate into two or more types and grow faster than isolated cells, i.e., R µ > 1 (Figs. 2 and S1). In category (b), interacting cells differentiate but their growth is slower than that of isolated cells (R µ < 1), which implies that some cell types exploit the other types (Fig. S2). In category (c), cells do not differentiate with respect to the average composition, but chemical concentrations asynchronously oscillate in time. In category (d), the behavior of each cell is identical, regardless of the presence or absence of cell-cell interactions, and therefore R µ = 1.
To characterize a prototype of MCO, we are mainly concerned with category (a), as this case enables both cell differentiation and cooperative growth. We found four common properties in this category. (1) A state with homogeneity among cells becomes unstable as the cell number increases, and is replaced by two (or more) distinct cellular states. tions are concentrated for only a few chemicals, whereas the concentrations of the other chemicals are nearly zero; i.e., each cell type uses only a sub-network of the total reaction network. (3) Different cell types share only a few common components, and the other components mostly exist in one cell type. (4) The components that predominate in one cell type diffuse to the other cell type, where they function as catalysts, and vice versa. Thus, the two cell types help each other to achieve higher cooperative growth.

Reduction of the N-cell model
After examining a number of networks in category (a), we extracted a common core structure in the reaction network topology, designated as networks 1-3 (Figs. 3A, 3B, and S3). In these networks, cells in the N -cell model differentiate into two types, type-1 and type-2, as exemplified in Fig. 3C. In type-1, x 1 is high while x 2 is close to zero, and in type-2, x 2 is high and x 1 is close to zero. Consequently, X 3 (X 4 ) can be produced only in the former (latter) type, and the two types of cells complement each other by exchanging X 3 and X 4 . Furthermore, the differentiated cells grow faster than the isolated cells (Fig.  3D).
To analyze the mechanism of this cooperative differentiation, we reduced the N -cell dynamics to two effective groups of cells represented by (x where v (i) denotes the total volume of each cell group (i = 1, 2). Considering that the total cell number is sustained at its maximum N , the total cellular volume is also bounded. Therefore, v = v (1) + v (2) is regarded as a constant in the reduced version of interacting cells, termed the reduced-2cell (r2cell) model. This model obeys Eqs. (1)- (2), and the effect of random cell elimination associated with cell division is implicitly incorporated into dilution due to volume growth. Besides, by considering the symmetry in networks 1-3, we can assume v (1) = v (2) = v/2 for symmetric differentiation with the same number of cells of the two types, while the case with v (1) =v (2) will be discussed later. Likewise, we also consider the reduced-1cell (r1cell) model corresponding to the "isolated" case of the N -cell model, by ignoring cell division and assuming that the cellular volume is constant at v (iso) = v.

(i) Cell differentiation
The behavior of the r1cell and r2cell models (i.e., isolated and interacting cells) can be classified into several phases, depending on parameters (C, V, D), where V ≡V med /v is the volume ratio between the cells and the medium. The phase diagram with network 1 for D = 1 is shown in Fig. 4A, and Fig. S5 shows phase diagrams of networks 1-3 for various D values. The blue area in Fig. 4A represents phase (I), in which the cells cannot differentiate, and always reach a single fixed point attractor in both the r1cell and r2cell models. In phase (II), differentiation into two fixed points occurs in the r2cell model from a stable fixed point in the r1cell model, as shown in Fig. 4B. In phase (III), the r1cell model exhibits oscillation, while two cells in the r2cell model reach two distinct fixed points (Fig. 4C). In terms of dynamical systems theory, this loss of oscillation is referred to as oscillation death [16,17]. In phase (IV), both "oscillationdeath" differentiation and synchronous oscillation (i.e., non-differentiation) can occur depending on the initial condition, whereas the r1cell model always exhibits oscillation. Thus, differentiation occurs in phases (II)-(IV) (i.e., at the left of the green line in Fig. 4A), while stable differentiation without falling into synchronized oscillation is achieved only in phases (II)-(III), i.e., with small C and V values, representing a limited resource and strong cellcell interaction condition. In Fig. 4A, phases (II)-(III) are divided by the red line, and the red and green lines are determined according to linear stability analysis (see Supplemental Material for details).
With respect to the network structure, the catalytic reactions X 1 + αX 2 →X 5 + αX 2 and X 2 + αX 1 →X 5 + αX 1 function as two mutually repressive reactions, i.e., forming a double-negative feedback loop. Further, the product X 5 consumes X 1 -X 2 and is maintained within the cell, which enhances the dilution of X 1 -X 2 . Thus, each of these reactions works as a composite negative feedback loop, leading to instability of the homogeneous cell state. Since nonlinearity is a necessary condition for multi-stability, a high order of catalytic reactions α tends to facilitate cell differentiation.
(ii) Cooperative growth by differentiation Figure 5A shows the dependence of R µ on parameters in network 1, exemplifying that differentiation increases the growth rate. Surprisingly, this differentiationinduced growth enhancement was always observed for any set of parameters in networks 1-3. We next sought to determine the mechanism contributing to the faster growth of differentiated cells. An intuitive explanation is as follows. On one hand, an isolated cell must contain all chemical components required for selfreproduction (e.g., X 0 -X 5 in the upper panel of Fig. 3A), leading to lower concentrations of each chemical on average. On the other hand, differentiated cells can achieve division of labor; each type of differentiated cell exclusively produces a portion of the required chemical species, and cells exchange these chemicals with each other. Since catalytic reactions occur only in a sub-network of the original network (e.g., a network in the lower panel of Fig. 3A), the chemicals are concentrated on fewer components, which increases the efficiency of chemical reac-tions and promotes cellular growth. This suggests that stronger cell-cell interactions support higher growth. Indeed, Fig. 5B shows that a smaller V , i.e., stronger cell-cell interaction, causes larger R µ . A smaller V also increases R p , the ratio of the total production of X 3 -X 4 in the r2cell model to that in the r1cell model (Fig. 5C); that is, the production of exchanged chemicals is enhanced. To conclude, stronger cell-cell interactions reinforce the division of labor, whereby differentiated cells can grow more efficiently. The rate of growth enhancement through cell differentiation R µ can be roughly estimated by recalling that the growth rate of a cell is given by the average influx of the nutrient chemical it receives. We compared the growth rate of an isolated cell µ (iso) to that of a differentiated cell µ (dif ) by assuming that the concentration of each chemical species is equally distributed, except for the nutrient chemical. Considering that an isolated cell has a catalytic network with k chemical components and q reaction paths from the nutrient X 0 , each concentration of X 1 -X k−1 is calculated as Hence, for the steady state, the growth rate is estimated by µ (iso) = qx (iso) 0 ). On the other hand, the sub-network in a differentiated cell is considered to have k ′ chemicals and q ′ reaction paths from the nutrient (k ′ < k, q ′ < q). Then, each chemical concentration and the growth rate are given by Here, we also assume that x , because these concentrations mostly depend on the supplied nutrient concentration C rather than on the internal dynamics of individual cells. From these assumptions, the growth ratio R µ ≡µ (dif ) /µ (iso) is calculated as R µ = (q ′ /q)[(k−1)/(k ′ −1)] α . For example, with network 1 or 2, k = 6, q = 4, k ′ = 4, q ′ = 2, and thus R µ = (1/2)(5/3) α , which is greater than unity, at least when α≥2. Although this estimate is not strictly accurate, it nevertheless demonstrates how cell differentiation can enhance cellular growth, which is facilitated by a greater α. Even when the chemical concentrations were non-uniform, division of labor could accelerate growth when α was sufficiently large.

(iii) Robustness in the number distribution of differentiated cells
The cells in our models achieved (i) cell differentiation and (ii) cooperative growth. However, this is not sufficient for the transition to an MCO, since if one cell type grows faster than the other type, the cooperation between the differentiated cells collapses. Thus, the third condition is necessary: the growth rate of each cell type needs to be in conformity, through mutual regulation by cell-cell interactions. Thus far, we have considered the case with equal populations of the two cell types by imposing the condition v (1) = v (2) . Here, we examine the case with v (1) =v (2) for fixed v (1) and v (2) , to evaluate whether the increases in cell volume (or number) are identical between the two types to meet this third requirement. Therefore, Figs. 6A and 6B show plots of the growth rate versus r (1) , where r (i) ≡v (i) /(v (1) + v (2) ) is the volume proportion between type-1 and type-2 cells. Now, let us denote the dependence of µ (1) on r (1) by a function F (r (1) ). Then, the growth rate of the type-2 cell µ (2) is given by G(r (2) ) = G(1 − r (1) ), which is equal to F (1 − r (1) ) due to symmetry in the catalytic network. Since differentiated cells help each other, balanced growth would be expected; if the volume or relative number of one cell type is larger than the other, a larger (smaller) amount of chemicals would be supplied from the majority (minority) type to the minority (majority) type, so that the growth rate of the minority type is enhanced compared to that of the other type. This is the case for network 2, where F (r (1) ) < F (1 − r (1) ) for r (1) > 1/2, and the difference in volume (or number) decreases over time, leading to a balanced cell distribution (Fig. 6B). Nonetheless, this is not always the case. Figure 6A shows that the growth rate of the majority type is larger than that of the minority type in network 1; that is, F (r (1) ) > F (1 − r (1) ) for r (1) > 1/2. Accordingly, the difference in volume increases over time, and thus the two different cell types cannot stably coexist. This instability of collective growth differs from a scenario of parasitic behavior, because µ (tot) ≡r (1) µ (1) + r (2) µ (2) is higher than µ (iso) for almost the entire range of r (1) . The condition for the stability is analytically expressed as follows. First, the temporal change of r (1) is repre- /dt = 0 is fulfilled in a cell with steady growth, and D med V ≫D. Then, from the definition of the growth rate µ, we get Since the first term is always negative, as described in the Supplemental Material, Eq. (3) shows that if the difference in exchanged chemicals between the majority and minority cells [ i )] increases in proportion to the increase in volume ratio of the majority type, then F ′ (1/2) is negative and thus the collective growth is balanced. The stability and instability of collective growth are also observed in the original N -cell model. With an "unstable" network 1, the N -cell model repeats the following dynamic behavior, as shown in Fig. 6C: the medium is dominated by cells of one type, and cells of the minority type become extinct. Then, their differentiated compositions cannot be maintained with a single cell type, leading to de-differentiation. Thus, the coexistence of differentiated cells is temporally regained.
In contrast, in a "stable" network 2, the two differentiated cell types stably coexist and their growth is balanced (Fig. 6D), in which a perturbation to increase the population of one cell type leads to a decrease in the growth rate of that type.

DISCUSSION
We have revealed how (i) cell differentiation, (ii) growth enhancement, and (iii) robustness in a cell population can emerge in a simple system consisting only of intracellular reaction dynamics and cell-cell interaction through chemical diffusion. First, cells sharing a common genotype (i.e., an identical reaction network and identical parameters for reaction and diffusion) differentiate into several types with different chemical compositions, as a result of instability of a homogeneous cell state, induced by cell-cell interaction. This differentiation is facilitated under a condition of limited resources and strong cell-cell interaction, given a high order of catalytic reactions. This dynamicalsystems mechanism has also been proposed in previous studies using models of metabolic networks [7] and gene regulation networks [18]. Second, the differentiated cells can achieve cooperative division of labor. The chemical abundances of each cell type are concentrated in only a few species, which increases the speed and efficiency of catalytic reactions, while the chemicals that are needed but lacking in a given cell are received from the cells of another type. Thereby, the differentiated cells enhance their growth cooperatively. This growth enhancement due to division of labor is clearly distinguishable from a scenario of relaxation of the competition for resources through niche differentiation reported previously [7], in which the growth rate is not increased relative to that of an isolated cell. The present mechanism of division of labor is similar to the economic theory of comparative advantage proposed by Ricardo [19], which claims that the mutual use of surplus from a different country is more advantageous than producing all necessary resources in a single country, unless the transport cost is too high. In fact, our mechanism works best in the case of a high density of cells, where chemicals are easily exchanged without much loss within the medium. Remarkably, this cooperative differentiation is not sufficient to satisfy condition (iii), robustness in the number distribution of differentiated cells. If one cell type begins to dominate the population, production of the chemicals needed by the minority type will increase. Thus, a feedback mechanism to reduce the majority population is expected. However, if the fraction of non-diffusible chemicals is increased for the majority cell type, this storage of chemicals within a cell would suppress the supply of chemicals to the other cell types. Consequently, the majority cell type would further increase its population. This suggests that to achieve a balanced population state, the mutual transport of necessary chemicals must work efficiently beyond any possible increase in internal reserves. Our results demonstrate that an aggregate of simple cells consisting only of catalytic reactions and diffusive transport of chemicals can fulfill the requisites for evolution and stability of an MCO according to a benefit in cell growth. This suggests that the emergence of MCOs does not require a fine-tuned mechanism, but is rather a natural course in evolution. Further, our results imply that the emergence of MCOs is facilitated by a strong cell-cell interaction, limited resources, and a high order of catalytic reactions. From this point of view, it is interesting to compare the present results to some characteristics of current, primitive MCOs. When deprived of nutrients, cellular slime molds [20,21] or myxobacteria [22] gather together and exhibit stable cell differentiation, similar to an MCO. Second, in Volvox, a mutant possessing only homogeneous cells grows slower than the wild type [23], which is consistent with our observation that differentiation leads to growth enhancement. Moreover, filaments of the cyanobacterium Anabaena are known to differentiate, with each filament metabolically specializing in photosynthesis or nitrogen fixation, enabling more efficient growth [24]. Third, the biofilm of Bacillus subtilis achieves division of labor and metabolic co-dependence between interior and peripheral cells by chemical oscillation [25], suggesting the relevance of nonlinear dynamics and cell-cell interactions for differentiation. In this context, our results are also related to the experimental emergence of MCO behaviors from unicellular organisms [10,11]. In the fitness alignment across different levels, a game theory approach is sometimes adopted to address the evolution and dynamics of conflict between individuals and society. From a game theory perspective, the cellular growth rate is regarded as a measure of fitness or score. Hence, when two cell types stably coexist in network 2, Nash equilibrium is achieved at r (1) = 1/2, while in network 1 no such equilibrium is achieved for 0 < r (1) < 1, and thus different cell types cannot stably coexist. In contrast, the Pareto optima, where cells grow at the highest rate possible, are achieved at r (1) = 1/2 with both networks 1 and 2. It follows that the unstable collective growth can be interpreted as an inconsistency between the Nash equilibrium and Pareto optimum. Interestingly, both cell types gain some benefit from their unstable coexistence, even though the benefit of one type is greater than that of the other. The instability of collective growth emerges even in a system containing only cooperators, and thus differs from observations in the well-studied cheater-cooperator systems [6,26]. This novel type of dilemma emerges from a dynamical system to alter the internal state depending on the cell type distribution. Such unstable behavior in a system of cooperating units will be observed in an artificial symbiosis experiment with Escherichia coli and diffusible amino acids [27,28]. Considering the difference between networks 1 and 2, the stability and instability of the system can be switched by even a slight change in the diffusibility of a single chemical species. This implies that slight epigenetic changes and transcriptional errors occurring during the lifetime of an organism can lead to instability in the cell distribution, which may relate not only to the evolution of MCOs, but also to the phenomena of metamorphosis during development and carcinogenesis.