A Mass Conserved Reaction–Diffusion System Captures Properties of Cell Polarity

Cell polarity is a general cellular process that can be seen in various cell types such as migrating neutrophils and Dictyostelium cells. The Rho small GTP(guanosine 5′-tri phosphate)ases have been shown to regulate cell polarity; however, its mechanism of emergence has yet to be clarified. We first developed a reaction–diffusion model of the Rho GTPases, which exhibits switch-like reversible response to a gradient of extracellular signals, exclusive accumulation of Cdc42 and Rac, or RhoA at the maximal or minimal intensity of the signal, respectively, and tracking of changes of a signal gradient by the polarized peak. The previous cell polarity models proposed by Subramanian and Narang show similar behaviors to our Rho GTPase model, despite the difference in molecular networks. This led us to compare these models, and we found that these models commonly share instability and a mass conservation of components. Based on these common properties, we developed conceptual models of a mass conserved reaction–diffusion system with diffusion–driven instability. These conceptual models retained similar behaviors of cell polarity in the Rho GTPase model. Using these models, we numerically and analytically found that multiple polarized peaks are unstable, resulting in a single stable peak (uniqueness of axis), and that sensitivity toward changes of a signal gradient is specifically restricted at the polarized peak (localized sensitivity). Although molecular networks may differ from one cell type to another, the behaviors of cell polarity in migrating cells seem similar, suggesting that there should be a fundamental principle. Thus, we propose that a mass conserved reaction–diffusion system with diffusion-driven instability is one of such principles of cell polarity.


Introduction
Eukaryotic cells such as neutrophils and Dictyostelium cells respond to temporal and spatial gradients of extracellular signals with directional movements [1][2][3][4][5][6]. This process, known as chemotaxis, is a fundamental cellular process [5,[7][8][9]. In a migrating cell, specific molecular events take place at the front and back edges [1,2,5,10]. The spatially distinctive molecular accumulation inside cells is known as cell polarity. The front-back polarity usually has one axis, and this uniqueness is an important property because a migrating cell with two fronts could not move effectively [11]. Another behavior of the front-back polarity is higher sensitivity of the front to a gradient of extracellular signals [10,12]. This would also be important because the direction of movement should be controlled at the front edge.
Many mathematical models that account for gradient sensing and signal amplification in cell polarity have been proposed [12]. The local excitation and global inhibition model has been proposed to explain spatial gradient sensing [6,18]. Some models involve positive feedback loops for amplified accumulation of signaling molecules [19][20][21][22][23]. A reaction-diffusion model that includes local self-enhancement and long-range antagonistic effects has been proposed for directional sensitivity [24]. Most of the reported models of cell polarity, which involve the detailed parameters such as concentrations or rate constants, have been constructed with many parameters and equations. Although these detailed models are at least partially successful in reproducing experimental observations in cell polarity, the theoretical essence underlying cell polarity has not been explicitly demonstrated; thus, a simple conceptual model that can be used for analytical study is needed to extract common principles in cell polarity. Although the reported models consist of distinct molecular species or networks, it should be especially emphasized that many of them are able to exhibit similar behaviors of cell polarity regardless of their different frameworks. This fact indicates that a common principle should underlie the models, and a conceptual model is suitable for extracting common principles in cell polarity.
Because the Rho small GTPases are key regulators for cell polarity [16,17], we first developed a reaction-diffusion model of the Rho GTPases on the basis of an earlier model [25] to examine the spatial properties of the Rho GTPases. We found that the interaction of the Rho GTPases per se can generate specific spatial accumulation of the Rho GTPases, and that our model shows important behaviors of cell polarity. We also found that our model exhibits behaviors similar to the model by Narang and Subramanian [22,23], which is based on the molecular networks that are different from ours. This suggests that common principles should underlie both models. We found that a mass conservation of components and diffusion-driven instability are commonly conserved in the Narang and Subramanian models and in our model. Based on these common properties, we established conceptual models of a mass conserved reaction-diffusion system, and found that such properties can account for the critical behaviors of cell polarity. These results strongly suggest that a mass conservation of components with diffusion-driven instability is one of the fundamental principles of cell polarity.

The Rho GTPases Model
We developed a reaction-diffusion model of the Rho GTPases (Rac, Cdc42, RhoA) on the basis of the earlier model of the Rho GTPase [25], which explains the temporal behaviors, to examine whether the interaction of the Rho GTPases can generate the spatial behaviors in the cell polarity of migrating cells. The Rho GTPases exhibit guanine nucleotide-binding activity and function as molecular switches, cycling between an inactive GDP(guanosine 59-bis phosphate)-bound form and an active GTP-bound form. The Rho GTPases in active forms are located in the plasma membrane, and those in inactive forms are in the cytosol ( Figure 1A) [26]. It is likely that molecules in the cytosol have larger diffusivity than those in the plasma membrane. According to some studies, Cdc42 activates Rac [27][28][29], and RhoA has mutual inhibitory interactions with Cdc42 and Rac [29][30][31][32][33]. In addition, Rac plays a dominant role in a positive feedback loop, which involves PI3K, PIP 3 , and F-actin [13,[34][35][36]. Based on these experimental findings, we developed a diagram of the Rho GTPases interaction ( Figure 1B). We assume that molecules of Rac, Cdc42, and RhoA are activated by guanine nucleotide exchange factors (GEFs; ka i ) and are inactivated by GTPase-activating proteins (GAPs; ki i ), and that interactions between molecules (k ij ) are additive to GEFs or GAPs. Some molecule-molecule interactions are stimulation dependent. Activations of molecules by the stimulation (ks i ) are also assumed to be additive to GEFs. As in many previous models [18][19][20][21][22][23], we describe the spatial kinetics of molecules by simple diffusion equations. A recent study in which the diffusion coefficients of the Rho GTPases in the plasma membrane are determined [37] may support this assumption. The model of the interaction of the Rho GTPases is as follows (see also Materials and Methods):

Author Summary
Eukaryotic cells such as neutrophils and Dictyostelium cells respond to temporal and spatial gradients of extracellular signals with directional movements. In a migrating cell, specific molecular events take place at the front and back edges. The spatially distinctive molecular accumulation inside cells is known as cell polarity. Despite numerous experimental and theoretical studies, its mechanism of emergence has yet to be clarified. We first developed a mathematical model of the Rho small GTP(guanosine 59-tri phosphate)ases that cooperatively regulate cell polarity, and showed that the model generates specific spatial accumulation of the molecules. Based on our Rho GTPases model and other models, we further established a conceptual model, a mass conserved reaction-diffusion system, and showed that diffusion-driven instability and a mass conservation of molecules that have active and inactive states are sufficient for polarity formation. We numerically and analytically found that molecular accumulations at multiple sites are unstable, resulting in a single stable front-back axis, and that sensitivity toward changes of a signal gradient is specifically restricted at the front of a polarized cell. We propose that a mass conserved reaction-diffusion system is one of the fundamental principles of cell polarity.
where Rac, Cdc, and Rho with suffixes m and c denote the concentrations of Rac, Cdc42, and RhoA in the active state (membrane) and inactive state (cytosol), respectively. The numerical suffixes represent the following: 1, Rac; 2, Cdc42; and 3, RhoA. Dm i and Dc i denote the diffusion coefficients of molecules in the active state and inactive state, respectively (Dm i , Dc i ). The position-dependent parameter, S, denotes the intensity of stimulation. Because the parameters have not been fully obtained experimentally, we set parameters to reproduce the behaviors of cell polarity (Figure 2), and further analyzed the generality of such behaviors in detail with conceptual models (see below). Reversible accumulation of the Rho GTPases. Cell polarity accompanies exclusive reversible accumulation of Cdc42 and Rac, or RhoA, at the fronts or backs of migrating cells, respectively [4,10,[13][14][15]. We examined the accumulation of the Rho GTPases in response to a shallow gradient of stimulation (Figure 2A-2C). The basal level of stimulation (S ¼ 0.1) did not induce exclusive accumulation of the Rho GTPases ( Figure 2A); however, the larger stimulation (averaged S ¼ 0.4, with a slight gradient) induced accumulation of Cdc42 and Rac toward the stimulation point with the maximal intensity (x/L ¼ À0.2), and exclusion of RhoA from the point ( Figure 2B). The responses of accumulation and exclusion were switch-like against the intensity of stimulation (unpublished data). The accumulation and exclusion were reversed when the intensity of stimulation reduced to the basal level (S ¼ 0.1) ( Figure 2C).
Multiple transient peaks and a single stable peak. We next gave two points simultaneous stimulation. Cdc42 and Rac accumulated toward the stimulation points with the maximal intensity (x/L ¼À0.25 and 0.25) at t ¼ 60 ( Figure 2E). However, one of the accumulation peaks eventually disappeared, and only a single peak remained ( Figure 2F). This result indicates that the model shows specific behaviors in cell polarity; ''uniqueness of axis'' (see below). The selection of the remaining peak differed between simulation tests (unpublished data).
Sensing of the stimulation gradient by the polarized peak. We next examined the behavior of the accumulation peak (polarized peak). The polarized peak moved toward the new stimulation point of the maximal intensity, rather than generating a new peak ( Figure 2G-2I). This result indicates that the polarized peak can sense a gradient of stimulation (see below).

Conceptual Model for Cell Polarity
The behaviors in the Rho GTPases model, such as switchlike reversible accumulation, uniqueness of axis, and sensing of the stimulation gradient by the polarized peak, were similar to those observed in the models of Subramanian and Narang [22,23]. Despite the differences of the molecular species and networks among the models, the similarity in behaviors among them raises the possibility that a common principle could underlie them. Therefore, we examined whether common properties can be seen.These models belong to reaction-diffusion systems with a periodic boundary condition and exhibit switch-like response, implying that instability is important for accumulation of the components. In addition, these models involve components whose masses are conserved. The total amount of phosphoinositides between the plasma membrane and the endoplasmic reticulum is conserved in Narang and Subramanian's model [22,23], and the total amounts of the Rho GTPase between the membrane (e.g., Rho m ) and cytosol (e.g., Rho c ) are conserved in our Rho GTPases model. Based on these common properties, we derived a new concept (i.e., mass conserved reaction-diffusion system with diffusion-driven instability) and hypothesized that this system is a fundamental principle of cell polarity. We therefore developed a simple conceptual model with two components (u and v), which belong to a mass conserved reaction-diffusion system with instability, and examined whether the model can cause the behaviors of cell polarity to emerge sufficiently.
Mass conserved reaction-diffusion system. We investigated the following class of reaction-diffusion system, which is composed of two components, u and v: where u and v denote the concentrations of u and v, respectively, at time t and at position x, and D u and D v denote the diffusion coefficients of u and v, respectively. The reaction term f is given by a function of u and v. Because the total quantity of u and v is conserved in this case, we refer to this system as the mass conserved reaction-diffusion system. We developed several models that belong to the mass conserved reaction-diffusion system and show similar results (unpublished data), and here focused on the following model because of its simplicity for computation.
where a 1 and a 2 are parameters of the model. The position-  (A-C) Reversible accumulation of the components. We set S ¼ 0.2 at t ¼À100, and ran the simulation until t ¼ 0 to reach the stationary state (A). Then, the stimulation S ¼ 1f1 þ 0.01 cos[2p(x/L þ 0.2)]g was given, and the simulation was run until t ¼ 1,000 to reach the stationary state (B). The stimulation was reduced to the basal level (S ¼ 0.2) again, and the simulation was run until t ¼ 1,200 to reach the stationary state (C). (D-F) Multiple transient and a single stable accumulation peaks. Starting with the conditions in (A), we gave random perturbation (60.02) to obtain the homogenous initial conditions (D). Then, at t ¼ 0, the stimulation S ¼ 1f1 þ 0.01 cos[2p(x/L þ 0.25) 3 2]g was given, and the simulation was run for the indicated time (E,F). The stimulation points with the maximal intensity were at x/L ¼ À0.25 and 0.25. (G-I) Sensing of the stimulation gradient by the polarized peak. We set S ¼ 1f1 þ 0.01 cos[2p(x/L þ 0.2)]g at t ¼À1,000, and the simulation was run to reach the stationary state where the polarized peak was seen at x/L ¼ À0.2 (G). The new stimulation S ¼ 1f1 þ 0.01 cos[2p(x/L)]g was given, and the maximal intensity point was shifted to x/L ¼ 0. The simulation was run for the indicated time (H,I). The arrow indicates the direction of movement of the polarized peak (H). doi:10.1371/journal.pcbi.0030108.g003 polarized peak ( Figure 3G-3I). This finding indicates that the conceptual model retains the essential behaviors in the Rho GTPases model. We further used this conceptual model to examine in detail two behaviors of cell polarity: uniqueness of axis and sensitivity of the polarized peak.
Numerical evaluation of uniqueness of axis. We analyzed unstable homogenous solutions to examine the characteristic behaviors of this system by setting large homogenous S (S ¼ 1) ( Figure 4A and 4B). The small perturbation was given to the homogenous stationary state. This perturbation transiently induced two peaks in Model I with L ¼ 10. However, one of the peaks eventually disappeared, and only a single peak remained ( Figure 4A). When we doubled the system size (L ¼ 20), four peaks appeared; however, as in the case with L ¼ 10, only a single peak remained and others disappeared ( Figure  4B). No peak arose after the first transient peaks appeared. This result suggests that the number of transient peaks depends on the system size, and that multiple transient peaks are unstable, which results in a single peak at the final steady state.
We analyzed the relationship between the number of transient peaks and the system size by linearization analyses around the homogenous states (see Materials and Methods). For Model I (Equation 2), the wave number that grows most rapidly from homogenous state k Ã h is obtained as k Ã h ¼ 1.32; the most growable wavelength is 2p/k Ã h ¼ 4.78. Therefore, the system with L ¼ 10 most likely generates two peaks and the system with L ¼ 20 four peaks from the homogenous states.
We further confirmed that multiple peaks are unstable. We set a two-peak profile as the initial condition, and examined its stability by giving a small perturbation. As shown in Figure  4C, the two peaks were not stable; one disappeared, and only a single peak remained. This result strongly suggests that multiple peaks are unstable, resulting in a single stable peak.
Numerical evaluation of localization of sensitivity. We examined whether the sensitivity is localized at the polarized The system size was taken as L ¼ 10 (A) or L ¼ 20 (B). We gave random perturbation (60.01), and the simulation was run for the indicated time (A,B). (C) Next, we obtained a stable single peak in Model I by setting L ¼ 5 and taking the initial state as u ¼ 1.65 and v ¼ 0.35. Because we applied the periodic boundary condition to this system, we could set the center of the concentration peak at x ¼ 0 by translation. By duplicating and coupling this profile (L ¼ 5), we obtained a new profile (L ¼ 10) with two peaks. We gave small perturbation (60.01) to this profile (L ¼ 10), and the simulation was run for the indicated time (C). (D,E) Numerical evaluation of localization of sensitivity. We obtained a stable single peak in Model I by setting L ¼ 10 and S ¼ 1, and taking the initial state as u ¼ 1.65 and v ¼ 0.35. We set the center of the concentration peak at x ¼ 1.0 (D) or x ¼À1.6 (E) by translation. Then we gave a new stimulation S ¼ 1 (x , 0) and S ¼ 1f1 þ 0.04cos[2p(x/L À 0.25)]g (x ! 0). The stimulation gradient is given only at the area of x . 0. The simulation was run for the indicated time (D,E). doi:10.1371/journal.pcbi.0030108.g004 peak. We first induced the stable polarized peak as indicated (t ¼ 0), and gave the indicated position-dependent stimulation ( Figure 4D and 4E). When the polarized peak (peak at x ¼ 1.0) overlapped with the gradient of the stimulation, the polarized peak moved along the gradient of S ( Figure 4D); a new peak did not appear. In contrast, when the polarized peak (peak at x ¼À1.6) did not overlap with the gradient of the stimulation, the polarized peak did not move ( Figure 4E); a new peak did not arise. These results indicate that the polarized peak shows sensitivity to a gradient of S, whereas the rest of the position does not show such sensitivity, and suggest that the sensitivity is restricted only at the polarized peak.

Analytical Examination of Unique Axis and Localized Sensitivity
To better understand the results of the numerical simulations, we investigated the following model (see Equations 1a-1c), by analytical approximations: where D u ¼ aD v . Here, a 1 , a 2 , and a are parameters of the models. This model belongs to the mass conserved reactiondiffusion system, and is more advantageous for analytical examination. The homogenous solution was unstable regardless of the values of parameters in this model (see Materials and Methods), so this model did not show reversible accumulation. However, the model still retained the important properties such as uniqueness of axis and localization of sensitivity, so we can use this model to analytically examine whether these behaviors can emerge from a mass conserved reaction-diffusion system with instability. In the following sections, we show that: (1) the model has one-peak stationary states, regardless of the system size (if not too small); (2) multiple-peak stationary states are unstable; and (3) the polarized peak moves depending on the gradient of the parameter value and the sensitivity is localized. Finally, (4), we verified our analyses by comparing analytical results with the values obtained by numerical simulations.
Analysis (1): Existence of a one-peak stationary solution. We define the following variables and function: Equations 1a and 1b are rewritten as the following set of partial differential equations for N and P: Under the periodic boundary conditions, the stationary solutions of Equations 5a and 5b, N e (x) and P e (x), satisfy the following equations: Consider the case of Model II given by Equation 3. For any P e (.0), Equation 6b with the substitution of Equation 3 has a family of periodic solutions with periods between k min , k , ' (see Materials and Methods). Here k min is given by (1/a 1 a 2 )g 1/2 . Thus, for an arbitrary large system with L . k min , one can choose a one-peak solution that satisfies the boundary conditions at x ¼ 6L/2 by taking k ¼ L.
Note that P e is related to k and the average mass of N e (x), For a sufficiently long period (k ! '), this is expressed as Thus, for a sufficiently large system, N e (x) and P e (x) are approximated as: where x p denotes the center of the peak and b, N 0 , and P e are constants given by: ð9bÞ Here we obtain the one-peak solution for Model II (Equation 3) by setting an arbitrary L and N.
Analysis (2): Stability of periodic solutions. We examined the stability of periodic solutions by linearization analysis for a single-peak solution (Analysis (2.1)) and multiple-peak solutions (Analysis (2.2)), whose periods are equal to the system size and the 1/n of the system size, respectively.
Analysis (2.1): Stability of a one-peak solution. First, we consider the stability of a one-peak solution given by Equations 8a and 8b. Without loss of generality, we set , and the stability is estimated by the linearized equations of Equations 5a and 5b around Equations 8a and 8b, given as follows: where h N (x) and h P (x) are partial derivatives of f*(N, P) by N and P, respectively, at the solution Equations 8a and 8b, that is, h N (x) ¼ @f*(N e (x), P e )/@N and h P (x) ¼ @f*(N e (x), P e )/@P. Let us represent (DN, DP) as (e lt n l (x), e lt p l (x)) and consider the case of Model II. Equations 10a and 10b lead to the following: We illustrate n l and p l in Figure 5A and 5B, respectively, for the case with l . 0 (see Materials and Methods). In the region where jxj is larger than 1/b, p l can be approximated as follows: where C 1 is a constant originated from the integral constant. This approximation is verified numerically, as shown by the dashed line with an arrow in Figure 5B. Considering periodic boundary conditions, p l (ÀL/2) and p l (L/2) should connect smoothly; thus, Equation 12 never satisfies Equations 11a and 11b except for l ¼ 0. Analysis (2.2): Instability of multiple-peak solutions. Here we consider the case of multiple-peak solutions. As long as L/n is larger than 1/b, an identical n-peak periodic solution is approximated by: on each jth domain defined as x j p À L/2n x x j p þ L/2n, where x j p is the position of each peak j (j ¼ 1, 2, ..., n) given by x j p ¼ÀL/2 þ (j À 1/2)L/n. We set DP ¼ P(x, t) À P j e (x) ¼ e lt p j l . As in the case of the one-peak solution, we can obtain for each jth domain. For multiple-peak solutions, n j 0 and C j 1 are determined by the boundary conditions such that p j l connect smoothly between adjacent domains. This piecewise linear function takes the value p j ¼ Àn j 0 n 2 P e /N 0 at x ¼ x j p , and at the point the slope changes by Àl(2N 0 /n 2 P e b)p j . A smooth connection between them means the slope is given as (p jþ1p j )/(L/n) for x j p , x , x jþ1 p . Thus, p j satisfies the following relationship: Considering periodic boundary conditions, we obtain the following solutions: where k is an integer (1 k n/2), corresponding to a mode, and h 0 is an arbitrary constant. All modes have positive l, and therefore can grow. Figure 5C-5F shows the most growable mode for the n-peak solution (n ¼ 2, 3, 4, 5), that gives the largest value to l. Based on this analysis, the one-peak solution is stable, whereas the multiple-peak solutions are unstable, regardless of the parameter values.
Analysis (3): Response of the polarized peak to a positiondependent parameter. Here we study the case in which the system has a position-dependent parameter by substituting a Ã 2 (x) ¼ a 2 þ ea e (x) for a 2 in Equation 3, where jea e (x)j is sufficiently smaller than a 2 . It is useful to represent the explicit dependence of f* on the parameter, as f*(N, P, a Ã 2 ). In this case, an existing peak moves, as is the case with Model I (Figure 3G-3I). To derive the velocity of the peak (v p ), we define the following variable and functions: Because we cannot determine the values of l and C 1 without boundary conditions, we arbitrarily set l ¼ 0.06 and C 1 ¼ 0.3 here. The dashed line with an arrow in (B) indicates the approximation to a piecewise linear function. Next, we show the perturbations for multiple-peak solutions. We can describe a perturbation by a set of p j , where p j is the value of p l at the center of the jth peak. As shown in the text, we obtain p j as p ðkÞ j ¼ cos[(2kp/n)j þ h 0 ], where k is an integer (1 k n/2), corresponding to a mode, and h 0 is an arbitrary constant. (C-F) Show the most growable perturbations for two-, three-, four-, and five-peak solutions, respectively. For each n, the mode, k, that gives the largest value to l, is determined by l (k) ¼ [(12n 3 D v a 2 )/(L 3 N 2 b)]sin 2 (kp/n). doi:10.1371/journal.pcbi.0030108.g005 Equations 5a and 5b are rewritten as the following set of partial differential equations for N* and P*: Under the periodic boundary conditions, the stationary solutions of Equations 18a and 18b, N Ã e (z) and P Ã e (z), satisfy the following equations: Considering that e ¼ 0 and v p ¼ 0 lead to N Ã e (z) ¼ N 0 sech 2 (bz) and P Ã e (z) ¼ P e , we can set N Ã e (z) ¼ N 0 sech 2 (bz) þ n e (z) and P Ã e (z) ¼ P e þ p e (z) for small value of e. The linearized approximations of Equations 19a and 19b around (N Ã e ; P Ã e ; a Ã 2 ) ¼ (N 0 sech 2 (bz), P e , a 2 ) give us the velocity of the peak at t ¼ 0 (see Materials and Methods): where W 1 (x) ¼À2tanh(bx)sech 2 (bx), and Z is a constant. G a (x) is defined by G a (x) ¼ (4b 2 N 0 /a 2 )sech 2 (bx)a e (x), which reflects the position dependence of a Ã 2 (x). Equation 20 indicates that the velocity of the peak is determined by the integral of W 1 (x)eG a (x) with respect to x. Taking into consideration that eG a (x) represents the position dependence of a Ã 2 (x) and that the position dependence at the site where the value of W 1 (x)is trivial has little influence on the velocity, we can regard W 1 (x), which is a kind of weight function, as a ''sensing window.'' As shown in Figure 6, this window has significant value only at the site of the concentration peak. Indeed, as for the mass conserved reaction-diffusion system, W 1 (x) is proportional to dN Ã e /dx, and thus the sensing window is locally open at the preexisting peak, especially at the peripheries of the peak. Note that the integral of W 1 (x)eG a (x) is zero when a e (x) is an even function; therefore, the sensing window can detect a gradient (or slope) of a Ã 2 (x). Based on this analysis, the polarized peak can detect the position dependence of a parameter and move depending on the gradient at the site of the peak. This property can be called ''localized sensitivity.'' Analysis (4): Verification of analysis by computations. Finally, in Analyses (4.1)-(4.3), we verified our Analyses (1)-(3), respectively, by comparing analytical results with the values obtained by numerical simulations.
Analysis (4.1): Approximation of the one-peak solution. First, we verify the results of Analysis (1): Existence of a one-peak stationary solution. We obtain analytical approximations of the one-peak solution as Equations 8a and 8b, and we compare this approximation with the final profile of the numerical simulation of Model II. The left panel of Figure 7A is the final profile of numerical simulation (see Materials and Methods). The solid line indicates the profile of N, and the dashed line indicates P. The right panel of Figure 7A shows the approximations given by Equations 8a and 8b, taking the center of the peak as x p ¼ À2. 15. The results of the computations show good agreement with our analytical results. The approximations were also sufficient when we set L ¼ 20, 40, and 80 (unpublished data).
Analysis (4.2): Instability of a two-peak solution. Next, we verify the results of Analysis (2.2): Instability of multiple-peak solutions. According to our analysis, a two-peak solution is unstable, and perturbations will grow exponentially with a growth rate l anl given by Equation 16b with n ¼ 2 and k ¼ 1, that is, l anl ¼ (96D v a 2 )/(L 3 N 2 b). We numerically examined the growth rate of the perturbation given to two-peak solutions with D v ¼ 1, 2, and L ¼ 20, 30, 40 (see Materials and Methods). We logarithmically plotted the growth rates estimated by analytical results (l anl ; Figure 7B, solid and dashed lines) and those obtained by numerical simulations (l sml ; Figure 7B, filled circles and squares) against the system size L. The results of the analyses and computations were nearly identical.
Analysis According to our analysis, the polarized peak will move when a parameter gains position dependence. A concentration peak formed in Model II (Equation 3) with uniform a 2 will move when a 2 is replaced by a Ã 2 (x) ¼ a 2 f1 þ (e/2)sin[2p(x/ L)]g with a velocity v anl obtained by Equation 20 as: We numerically examined the velocity of movement of the polarized peak in response to a parameter gradient with D v ¼ 1, 2, and e ¼ 0.02, 0.04, 0.06 (see Materials and Methods). We plotted the velocities estimated by analytical results (v anl ; Figure 7C, solid and dashed lines) and those obtained by numerical simulations (v sml ; Figure 7C, filled circles and squares) against the parameter gradients e. The results of the analyses and computations were nearly identical.

Rho GTPases and Conceptual Model
In this study, we used a mathematical model to clarify the role of the interaction between the Rho GTPases in cell polarity and developed a conceptual model of cell polarity to glean a theoretical understanding of the unique behaviors of cell polarity.
The Rho GTPases regulate the remodeling of the actin cytoskeleton via actin polymerization, depolymerization, and myosin activity [15,17,26], which ultimately establishes cell polarity. Then, what regulates the spatial activity of the Rho GTPases? The model proposed by Sakumura et al. indicates that the interaction of the Rho GTPases can regulate their own temporal activities [25]. We demonstrated that the interaction of the Rho GTPases can regulate their own spatial activities. In reality, the interaction of the Rho GTPases can provide more complicated temporal and spatial regulation of their activities. Further study, including the determination of kinetic parameters of the interaction, is necessary to develop a more realistic model of the Rho GTPases.
The activator-inhibitor model for pattern formation [38] and the local excitation and global inhibition model for directional sensing in chemotaxis [6] belong to conceptual models, rather than to detailed biological models. Such conceptual models with reduction of equations and parameters make analysis simpler and clearer. We identified a mass conserved reaction-diffusion system with instability as common properties between the cell polarity models. The model belonging to this system can sufficiently reproduce the important behaviors of cell polarity, such as uniqueness of axis and localization of sensitivity, and enabled us to theoretically understand such behaviors, which are difficult to examine without the models.

Uniqueness of Axis and Localization of Sensitivity
When interleukin-8, a chemoattractant, is applied simultaneously from two directions at a 458 angle, normal neutrophils choose one direction for migration instead of responding to both sources [11]. Neutrophils with multiple leading edges are rarely observed under normal conditions [39]. When HL60 cells are transfected with a dominantnegative Rho construct or treated with Rho-kinase inhibitors, many cells exhibit the multiple pseudopods, where, in some cells, protrusions gradually withdraw, leaving a single, prominent pseudopod [10]. In addition, inhibitions of PI3Ks cause HL-60 cells to form multiple pseudopods, which are weak and transient [39]. These results suggest the instability of multiple leading edges, which may make the front of a migrating cell single and stable. Chemotactic cells must have only one front-back axis because multiple fronts would prevent fine migration. Subramanian and Narang investigated the response of their model to two almost identical stimulations [23] and showed that only one of the two peaks that arise persists, which agrees with our results. Here we show that uniqueness of axis emerges from instability of multiple peak solutions in a mass conserved reactiondiffusion system (Figures 4C and 5).
In neutrophils [40], HL-60 cells [10], and Dictyostelium cells [41], the polarized cells respond to changes in direction of a gradient by performing U-turns rather than by simply reversing polarity. In addition, polarized migrating cells move forward without responding to the chemoattractant source near their rears [12]. These experimental findings indicate that the sensitivity to chemoattractants is localized at the leading edge of polarized cells. The localized sensitivity focuses the activity of the actin cytoskeleton at the leading edge, resulting in faster movement toward a chemoattractant source [12]. Few mathematical theories, however, have been proposed to explain the localization of sensitivity. Here we show that localization of sensitivity depends on the specific localization of a sensing window at the polarized peak in a mass conserved reaction-diffusion system ( Figure 4D and 4E and Figure 6). It should be added that many other systems can also exhibit a localized sensing window.

Biological Interpretation of Mass Conservation
Consider a molecule that satisfies the following conditions: (1) the molecule (X) has two states (Xm and Xc); (2) the total amount of X is conserved; and (3) the diffusion coefficient of Xc is larger than that of Xm. Two states of this molecule can be treated as components of a mass conserved system. Some kinds of small GTPases, such as those of the Rho family, have two forms, active and inactive forms; the Rho GTPases in the active forms are located in the membrane, and those in the inactive forms are in the cytosol [26]. Some enzymes involved in the cell polarity of chemotactic cells, such as PI3K and PTEN, are also reported to show a relationship between their activity and membrane binding [42][43][44][45][46]. Molecules in the cytosol may well diffuse faster than those in the plasma membrane. Thus, these molecules can be considered as components of mass conserved systems.
Chemotactic cells, such as Dictyostelium cells and neutrophils, polarize within a few minutes (30 s to 3 min) after they are exposed to chemoattractants [6,10,47]. Because the time scale of cell polarity is likely to be much shorter than that of gene expression and protein synthesis, we can assume that the masses of molecules are constant during the polarization of chemotaxis.

Instability of Multiple-Peak Solutions
We numerically and analytically show that multiple-peak solutions are unstable. To facilitate an understanding of the physics of this instability, we attempt to give an intuitive physical explanation of the behavior of molecules in the case where there exist two peaks (Figure 8), just as in Figure 4C. We simplify the situation as follows. (1) There are two spaces (each space has one peak). (2) The molecules have two forms, u and v, which have small and large diffusivities, respectively.
(3) No molecule is generated or degraded in the spaces. (4) The molecules move between spaces, mainly in v-form, depending on the concentration gradient of v-form molecules. (5) The u-form molecules convert v-form molecules into u-form, and this positive feedback is so strong that infusion of molecules into the space causes a decrease in vform molecules. Here, consider that a few molecules move from one space (S1) to the other (S2). According to (5), because of the nature of the positive feedback, the number of v-form molecules in S2 declines as the total number of molecules in S2 increases. In turn, according to (4), the declining of the number of v-form molecules in S2 successively facilitates further transfer of v-form molecules from S1 to S2, resulting in the further increase of the total number of molecules in S2. This flux is never disrupted by the generation or degradation of molecules because of the mass conservation. Therefore, two peaks in such a system are unstable. The condition (5) seems critical for instability, and we analyzed such conditions mathematically elsewhere [48].

Specific Feature of the Mass Conserved Reaction-Diffusion System
One of the most extensively studied reaction-diffusion models is the Turing model, in which robust spatial patterns, such as stripes or spots, emerge via a diffusion-driven instability [38,49]. An ordinary Turing pattern in 1-D space is stripes with an intrinsic scale length [50]. Mass conserved models also generate multiple peaks from the homogenous state during the early phase, which is explained by diffusiondriven instability. However, they exhibit characteristic transitions after the initial peaks arise: most peaks become smaller and eventually disappear, and only one peak persists ( Figure 4A and 4B). Why is the behavior of the mass conserved system so different from that of ordinary Turing models? Consider a reaction-diffusion system with vast size (L ! ') and interval I [x 1 , x 2 ] within the system, where x 1 and x 2 are arbitrary but far apart. Can we predict what will happen to interval I? For an ordinary Turing model, the linearization analysis around the homogenous solution gives us sufficient information [50]. The mass conserved system is more complex, however, because the behavior differs between the case where the components flow into interval I versus the case where they flow out, and we cannot predict which will occur. This difference in predictability seems to be fundamentally linked to the different behavior and the specificity of the mass conserved system. Mass conserved models have multiple stationary states that are spatially homogenous or periodic, including a one-peak state and multiple-peak states. We show that the multiple-peak stationary states are unstable, resulting in a single stable peak.
In some reaction-diffusion systems, such as the activatorinhibitor model (or substrate-depleted model), high diffusivity of inhibitor (or substrate) make multiple peaks unstable, resulting in a single stable peak [51]. An intuitive explanation for this instability is that an inhibitor, which is generated at the peak, rapidly diffuses throughout the cell. In the cell where the inhibitor can be generated and degraded, the spread of the inhibitor requires the large diffusivity to overwhelm the inhibitor degradation. On the other hand, in a mass conserved reaction-diffusion system, where no molecule is generated or degraded, a peak takes up molecules from its surroundings to grow. That is, the growing peak inhibits the system not by spread of inhibitor but by deprivation of molecules. In this case, large diffusivity is not required because there is no competitor to overcome, such as generation of molecules, and the inhibition can eventually spread throughout the cell. Indeed, Equation 16b clearly indicates that any D v can make multiple peaks unstable (l . 0), at least in Model II.
It may be counterintuitive that any mass conserved system finally exhibits a one-peak pattern. For Model II, the final steady state was a one-peak solution regardless of the system size, even when L was infinitely large. But for Model I, the final steady state had two peaks when we set L ¼ 80 (unpublished data). Some mass conserved models probably have a maximum size to have a unique peak. However, this maximum size is independent of the linearization analysis. The conditions for the uniqueness of concentration peak will be elucidated in future analyses.

Future Perspective and Conclusion
Because properties observed in simple models are expected to be conserved in more detailed models, we assume that movement of molecules follows a simple diffusion equation in our conceptual model. Active transport systems, such as actomyosin-and microtubule-based active transports, regulate cell polarity in various cellular processes [52]. Such active We set the two-peak solution (S1 and S2, top panel), which is similar to Figure 4C. We consider the following conditions. There are two spaces (each space has one peak). The molecules have two forms, u and v, which have small and large diffusivities, respectively. No molecule is generated or degraded in the spaces. The molecules move between spaces, mainly in v-form, depending on the concentration gradient of vform molecules. The u-form molecules convert v-form molecules into uform, and this positive feedback is so strong under the condition where @P e /@N , 0 [48] that infusion of molecules into the space causes a decrease in v-form molecules. Infusion of molecules to S2 accompanies the removal of the same number of molecules from S1 because of the mass conservation (bottom panel). doi:10.1371/journal.pcbi.0030108.g008 transports are likely to follow the formation of intracellular asymmetry, which takes place under the resting condition where the cell polarity is yet to be generated. Under such conditions, the diffusion of the Rho small GTPases has been measured and shown to be approximated by an apparent simple diffusion, if viewed on the order of seconds or tens of seconds [37]. Because our concern in this study is the earlier asymmetry formation rather than the completion of the polarity, which includes active transport, we here assumed the simple diffusion of the Rho GTPases. However, we will readily incorporate the detailed mechanism of the transport system of the Rho-GTPases in a future model.
In this study, we focused on the stationary state, but not on the transient state, which involves adaptation in response to a transient signal [44,[53][54][55]. Such properties, as well as high dimensionality and multiple components, should be incorporated into a future model.
Although our model is rather simple, it shows the important properties of cell polarity such as switch-like reversible response, uniqueness of axis, and localization of sensitivity. We further demonstrated that the instability of multiple-peak solutions and the specific localization of a sensing window at the polarized peak in a mass conserved reaction-diffusion system are responsible for uniqueness of axis and localization of sensitivity, respectively. One remarkable feature of a mass conserved reaction-diffusion system compared with other models is that a mass conserved reaction-diffusion system does not require strict assumptions for diffusion coefficients, such as smaller and much larger diffusivities of excitatory and inhibitory molecules, respectively, in the local excitation and global inhibition models [6,18], or an extraordinary large diffusivity of the inhibitor in the Gierer-Meinhardt model [51]. Since the Rho GTPases, PI3K, or PTEN have thus far not been demonstrated to involve such ad hoc assumptions of diffusivity, a mass conserved reaction-diffusion system is more likely to explain cell polarity where these molecules are involved, and to be adapted to a wide range of cell polarity. Taking into consideration that the Rho GTPases system satisfies conditions of a mass conserved reaction-diffusion system, it is likely that this system is one of the fundamental principles of cell polarity.

Materials and Methods
Simulation. We considered a one-dimensional circular system with circumference L. The position is represented by x (ÀL/2 x L/2). We applied the periodic boundary condition, which is used in some models that explain cell polarity [23,24]. We used explicit difference methods to perform simulations. The difference intervals for calculations are shown in the following text.
Simulation of the Rho GTPases model. The parameter values were set as follows: L ¼ 10, The difference intervals for calculations were taken to be Dt ¼ 0.01 and Dx ¼ 0.33. We set X m (x) ¼ 0.3, X c (x) ¼ 0.7 (X ¼ Rac, Cdc, Rho), unless specified.
Simulation of Model I. The parameter values were set as follows: a 1 ¼ 2.5, a 2 ¼ 0.7, D u ¼ 0.01, D v ¼ 1, and L ¼ 10. The difference intervals for calculations were taken to be Dt ¼ 0.01 and Dx ¼ 0.2. We set u ¼ 1 and v ¼ 1, unless specified.
Simulation of Model II for Analysis (4). The parameter values were set as follows: The difference intervals for calculations were taken to be Dt ¼ 0.005 and Approximation of a one-peak solution. The computation was performed by setting L ¼ 10 and D v ¼ 1 and taking the initial state as u ¼ 1 and v ¼ 1 with small perturbation (60.01). The final profile (t ¼ 200) is shown in the left panel of Figure 7A. The solid line indicates the profile of N, and the dashed line indicates P.
Instability of a two-peak solution. We examined the instability of twopeak solutions. First, we obtained a stable one-peak pattern in Model II (Equation 3) by taking the size to be L/2, where L ¼ 20, 30, 40, and taking the initial state as u ¼ 1 and v ¼ 1 with small perturbation (60.01). Because we applied the periodic boundary condition to this system, we could set the center of the concentration peak at x ¼ 0 by translation. Next, by duplicating and coupling this profile (L/2), we obtained a new profile (L) with two peaks. We used this profile (L) with small perturbations (60.01) as the initial state of the following simulation. All trials (D v ¼ 1, 2 and L ¼ 20, 30, 40) showed instabilities of two-peak profiles, and we obtained the growth rate, l sml , from the change in peak height.
Movement of a polarized peak in response to the parameter gradient. We examined the response to the parameter gradient. First, we obtained a stable one-peak pattern in Model II (Equation 3) by setting L ¼ 10 and taking the initial state as u ¼1 and v ¼ 1. We set the center of the concentration peak at x ¼ 0 by translation. Next, we substituted a Ã 2 (x) ¼ a 2 f1 þ (e/2)sin[2p(x/L)]g for a 2 in Equation 3. All trials (D v ¼ 1, 2 and e ¼ 0.02, 0.04, 0.06) showed movement of the polarized peaks, and we obtained the velocities, v sml , from the results.
Linearization analysis around the homogenous state. In the homogenous stationary state, the Jacobian matrix for the reaction terms is given by where f u and f v denotes the partial derivatives of f by u and v, respectively, at a homogenous stationary state of Equations 1a and 1b. We obtained a condition for instability of the homogenous solution: For example, the homogenous solution is stable with S ¼ 0.2 in Model I, whereas unstable with S ¼ 1, under the following conditions: The homogenous solution in Model II is always unstable. Through the stability analysis using J, the range of wave numbers (k h ) that have positive eigenvalues is obtained as 0 , k h , [(D u f v À D v f u )/ (D u D v )] 1/2 , and the wave number that has the largest eigenvalue and grows most rapidly from the homogenous state, k Ã h , is obtained as follows: For Model I (Equation 2), we obtain k Ã h ¼ 1.32 under the following conditions: Periodic solution of Equation 6b. Equation 6b has the same formulation as classical Newton mechanics. We define V(N e ) as VðN e ; P e Þ ¼ À and Equation 6b implies where E is a constant value, corresponding to period and total mass of N e (x). The period k and the average mass N ¼ (1/k) R k 0 N e (x)dx satisfy the following equations: where N min and N max are minimum and maximum levels of N e (x), respectively, and are derived from V(N min ) ¼ V(N max ) ¼ E (0 , N min , N , N max ). Equations 26a and 26b give the relationship among P e , k, and N, where N is straightforwardly derived from the initial condition of (u, v). For Model II (Equation 3), and E can range between E * , E , 0 for N e (x) to be a periodic solution. Here E * ¼ Àf[(D v À D u )/D u D v ][(D 2 v a 1 a 3 2 )/6P 2 e ]g. As E becomes smaller (E ! E * ), the period k converges to k min , which is the shortest wavelength in the periodic solutions. As E becomes larger (E ! 0), the period k diverges. The solution of N e (x) at E ¼ 0 corresponds to the separatrix of Equation 6b, indicating an infinite period (k ! '). The explicit form of N e (x) for E ¼ 0 can be obtained by which has the sole peak at x ¼ x p and decays to zero as x ! 6'. For a sufficiently large system, Equation 28 is a good approximation of the solution for ÀL/2 , x , L/2. Equations 28 and 7 lead to Equations 8a and 8b, and 9a-9c. Solution of Equations 11a and 11b. If there is nontrivial (n l (x), p l (x)) that satisfies Equations 11a and 11b for l with a positive real part, the solution is unstable. Note that (n l0 , p l0 ) ¼ (n 0 sech 2 (bx), À n 0 P e /N 0 ) satisfies Equations 11a and 11b for l ¼ 0 under periodic boundary conditions. Here n 0 is an arbitrary factor, originated from the linearity of equations, and we can set n 0 ¼ 1. For l with an absolute value near zero, we can obtain (n l , p l ) by the expansion from (n l0 , p l0 ) with regard to l. To do this, we take n l ¼ n l0 þ ln l1 þ . . . and p l ¼ p l0 þ lp l1 þ . . . . In the first order of the expansion, (n l1 , p l1 ) obeys the following equations: Thus, p l is immediately derived from Equation 29a as: where C 1 and C 2 are integral constants. We can obtain n l by solving Equation 29b with substitution of p l1 . C 2 is determined by the mathematical condition that n l1 should be orthogonal to n l0 . In practice, C 2 has little influence on (n l , p u ), and we set C 2 ¼ 0 in the analysis.
Derivation of v p from Equations 19a and 19b. The linearized approximations of Equations 19a and 19b are given as follows: 2v p N 0 tanhðbzÞsech 2 ðbzÞ ¼ @ 2 p e ðzÞ @z 2 ; ð31aÞ 0 ¼ ðD v þ D u Þ½2v p N 0 tanhðbzÞsech 2 ðbzÞ À D u D v @ 2 n e ðzÞ @z 2 þ ðD v À D u Þ½h N ðzÞn e ðzÞ þ h P ðzÞp e ðzÞ þ h a ðzÞea Ã 2 ðxÞ: ð31bÞ where h N (z) ¼ @f*(N e (z), P e ,a 2 )/@N, h P (z) ¼ @f*(N e (z), P e ,a 2 )/@P, h a (z) ¼ @f*(N e (z), P e ,a 2 )/@a Ã 2 . Because t is no longer a variable, we set t ¼ 0 without loss of generality and replace z with x in the following analysis. Equation 31a under the periodic boundary condition leads to p e as follows: where C 3 is an integral constant. By substituting Equation 32 into Equation 31b, we obtain the following: d 2 n e dx 2 ¼ 4b 2 ½1 À 3sech 2 ðbxÞn e þ v p bN 0 G n ðxÞ þ eG a ðxÞ; ð33Þ where G n (x) and G a (x) are defined by: Considering the periodic boundary condition for n e (x), we obtain the following equation for sufficiently large L (solvable condition): ÀL=2 W 1 ðxÞ½v p bN 0 G n ðxÞ þ eG a ðxÞ dx ¼ 0: ð37Þ This leads to the velocity of the peak: where Z is given as follows: