Necessary and sufficient condition for hysteresis in the mathematical model of the cell type regulation of Bacillus subtilis
aa r X i v : . [ q - b i o . CB ] S e p Necessary and sufficient condition for hysteresis in themathematical model of the cell type regulation of
Bacillus subtilis
Sohei Tasaki a,b, ∗ , Madoka Nakayama c , Izumi Takagi d,e , Wataru Shoji f a Institute for the Advanced Study of Human Biology (WPI-ASHBi), Kyoto University Institute forAdvanced Study (KUIAS), Kyoto University, Japan b RIKEN Center for Biosystems Dynamics Research (RIKEN BDR), Japan c Department of General Engineering, Sendai National College of Technology, Japan d Mathematical Institute, Tohoku University, Japan e Institute for Mathematical Sciences, Renmin University of China, China f Frontier Research Institute for Interdisciplinary Sciences (FRIS), Tohoku University, Japan
Abstract
The key to a robust life system is to ensure that each cell population is maintainedin an appropriate state. In this work, a mathematical model was used to investi-gate the control of the switching between the migrating and non-migrating statesof the
Bacillus subtilis cell population. In this case, the motile cells and matrixproducers were the predominant cell types in the migrating cell population andnon-migrating state, respectively, and could be suitably controlled according tothe environmental conditions and cell density information. A minimal smoothmodel consisting of four ordinary differential equations was used as the mathe-matical model to control the
B. subtilis cell types. Furthermore, the necessaryand sufficient conditions for the hysteresis, which pertains to the change in thepheromone concentration, were clarified. In general, the hysteretic control of thecell state enables stable switching between the migrating and growth states of the
B. subtilis cell population, thereby facilitating the biofilm life cycle. The resultsof corresponding culture experiments were examined, and the obtained corollarieswere used to develop a model to input environmental conditions, especially, theexternal pH. On this basis, the environmental conditions were incorporated in asimulation model for the cell type control. In combination with a mathematical ∗ Corresponding author.
Email address: [email protected] (Sohei Tasaki)
Preprint submitted to arXiv September 15, 2020 odel of the cell population dynamics, a prediction model for colony growth in-volving multiple cell states, including concentric circular colonies of
B. subtilis ,could be established.
Keywords: Bacillus subtilis , Cell type, Hysteresis
1. Introduction
Cellular state diversity is the source of the morphology and function of lifesystems. Even in prokaryotes, the robust structure of bacterial biofilms can beattributed to the heterogeneous presence of various types of cells (Stoodley et al.,2002; Donlan, 2002; Hall-Stoodley et al., 2004; Branda et al., 2005; Kobayashi and Iwano,2012; Vlamakis et al., 2013; Hobley et al., 2015; Flemming et al., 2016). Amongsuch cells,
Bacillus subtilis is considered to be the master of differentiation, as itcan exhibit extremely diverse cell types (L´opez et al., 2009). The cellular diver-sity helps develop a diverse colony morphology (Wakita et al., 1994) and com-plex biofilm structure, supporting long-term survival and growth in response toenvironmental variations (Tasaki et al., 2017a). In this study, we considered thefollowing mathematical model that describes the cell type regulation of
B. subtilis (Tasaki et al., 2020): dSdt = c S Ha S + b S C − d S SdHdt = c H a H + b H A − d H HdAdt = c A a A + b A S − d A AdCdt = c C X A − d C C (1)In this model, the cell states are described by four variables S = S ( t ) , H = H ( t ) , A = A ( t ) and C = C ( t ) , each of which represents a group of cooperating genesand their products (Fig. 1A). Specifically, S is a group represented by Spo0A ∼ P(including the phosphorelay of Spo0F, Spo0B, and Spo0A); H and A correspondto SigH and AbrB, respectively; and C is a group represented by ComK (drivenby the ComX-ComP-ComA pathway).The output of this system is the cell type, which corresponds to a matrix pro-ducer and motile cell when S is high and low, respectively (Fig. 1A). Sporula-tion is initiated when S exhibits continuously high values; however, this aspectwas not considered in this work. Instead, this work was focused on examining2 H SA X
External signal
Cell
Low pH(7.4 - 5.0) S : high S : low Matrix producer Motile cell
Anhysteretic
S X
Hysteretic
Matrix producerMotile cell
S X X X A B C
Figure 1: Model to determine the response of the cells and cell populations to the environmentalpH. (A) Model for the cell type selection against the environmental pH and cell density. (B and C)Two types of cell type controls. Anhysteretic (B) and hysteretic (C). the switch between the migrating (planktonic) and non-migrating (biofilm) states(Kearns and Losick, 2005; Kobayashi, 2007; Chai et al., 2008; L´opez et al., 2009;Cairns et al., 2014). The inputs of the cell state control system include the ex-ternal environmental conditions and auto-inducing signals that represent the celldensity. Among such signals, one corresponds to a small peptide ComX secretedby
B. subtilis cells. In the following text, the concentration of this peptide is de-noted as X , and it is an external parameter input from C to the cell type regulationsystem (1). The curve of the set of equilibrium points of (1) can be divided intotwo types pertaining to the increase and decrease in X , indicating that the cellstate (for example, S ) is monotonic and non-monotonic (Figs. 1B and C), respec-tively. In the latter case, the curve is an S-shaped curve with two turning points.Specifically, there exist two cases in which the cell type control is not hysteretic(Fig. 1B) and hysteretic (Fig. 1C), which pertains to the increase/decrease in thecell density signal X , respectively. In general, the hysteretic control is necessaryto facilitate the biofilm life cycle or concentric colony formation. When hysteresisexists in the choice of cell state, in the context of an increase and decrease in X ,a life cycle for the cell population occurs as follows: (Migration phase) The celldensity information X at the growth front of the cell population decreases as it dis-perses in the motility state, and the state switches to the matrix-production statewhen the concentration falls below a certain threshold X . (Growth phase) As thecell population grows and matures in the matrix production state, X increases, and3hen it exceeds a certain threshold X ( > X ) , the state switches to the motilitystate. In this regard, the objective of this study was to classify the presence or ab-sence of hysteresis in the selection of such cell states through certain parameters.Furthermore, the mechanism to control the state of the cells and cell populationsin response to environmental conditions was discussed.
2. Results
The main theorem described herein is a mathematical and formal claim. Inthis context, the meaning of the parameters may be difficult to understand. Thistype of unbiased form of writing is intended to facilitate the subsequent testing oftwo different interpretations.To describe the results, first, the steady state hysteresis is defined. A set ofsteady states is considered to be anhysteretic if the steady state is unique to X , andthe steady state S decreases monotonically with respect to X (Fig. 1B). In contrast,a set of steady states is considered to be hysteretic if the steady state is not uniqueto X , and the steady state S is an (inverse) S-shaped curve (Fig. 1C). When thesteady state is anhysteretic, any equilibrium point on the curve X = X ( S ) of theset of steady states is stable. In comparison, when the steady state is hysteretic,any equilibrium point between the two folding points ( X ′ ( S ) =
0) is unstable, andthe outer equilibrium point is stable (Appendix B).The parameters for the classification can be defined as follows: D = a A − b H c A a H d A , ˆ P = a H b H a A c A d A ( a H a A d A − b H c A ) , ˆ Q = (cid:26) c S a H c H b A d A a S d S a H d H ( b H c A − a H a A d A ) − (cid:27) . (2)In this case, the following holds true. Theorem 1.
If D ≥ , the set of steady states is anhysteretic (Fig. 1B). If D < ,any ˆ P admits ˆ Q c ( ˆ P ) such that If ˆ Q ≤ ˆ Q c ( ˆ P ) , the set of steady states is anhysteretic (Fig. 1B). If ˆ Q > ˆ Q c ( ˆ P ) , the set of steady states is hysteretic (Fig. 1C).Moreover, the threshold ˆ Q c = ˆ Q c ( ˆ P ) satisfies the following (Fig. 2). ˆ Q ′ c ( ˆ P ) > , ˆ Q ′′ c ( ˆ P ) > , lim ˆ P → + ˆ Q c ( ˆ P ) ˆ P / = , lim ˆ P → + ∞ ˆ Q c ( ˆ P ) ˆ P = . (3)4 AnhystereticHysteretic -3-2-1 0 1 2 3 4 5 6 7 8 9 10-4 -3 -2 -1 0 1 2 3 4
AnhystereticHysteretic
A B
Figure 2: Range of hysteresis in the ˆ P − ˆ Q plane. (A) Normal graph. (B) Log-log graph. A focus of this study was to examine the relation between the cell type cyclegeneration conditions and the environmental factors, specifically, to examine themechanism using which each cell of
B. subtilis reflects the environmental changesin the control pattern of the cell type. In this regard, we considered the environ-mental pH as a sample environmental factor. As mentioned previously, when thecontrol of two cell types, motile cells and matrix producers, is hysteretic, a life cy-cle occurs in the cell population growth process, known as the biofilm life cycle.One of the simplest observations of the cell population life cycle is the forma-tion of concentric colonies (Fujikawa, 1992; Itoh et al., 1999; Wakita et al., 2001;Shimada et al., 2004; Yamazaki et al., 2005). This colony growth pattern alter-nates between growth and migration phases. The dominant cell type for colonygrowth periodically switches between matrix producers and motile cells. In otherwords, the formation of concentric colonies depends on the presence or absenceof hysteresis in the cell type control. In a recent study, the relationship betweenconcentric colony formation and environmental pH was clarified (Tasaki et al.,2020). It was noted that under appropriate conditions, concentric colonies areformed on a solid nutrient medium containing approximately 0.7% agar. Initially,in the neutral region (pH 6.8–8.0) with an intracellular pH of 7.4 (Shioi et al.,1980), concentric colonies are not formed, and only the growth phase through thematrix production cells is observed (Fig. 3A). As the environmental pH decreasesto approximately 6.8, many extremely short migration phases appear. When thepH is close to this transition point, there exists a considerable variation in space,the periodicity is not clear, and the pattern is considerably different from concen-5
H6.86.55.35.1
No colonies Hysteretic Anhysteretic
Matrix producerMotile cell Matrix producer α SH Anhysteretic Hysteretic
Motile cellMatrix producer
Low pH (7.4 - 5.0) ? α A α A, L α A, M α A, U Anhysteretic Hysteretic Anhysteretic
Motile cellMatrix producerMatrix producer
Low pH (7.4 - 5.0)
Motile cell α A α A, L α A, U BD CA
Figure 3: Hysteresis and parameter change. (A) Hysteresis and environmental pH. The schematicpertains to concentric colony formation experiments (Tasaki et al., 2020). (B) Hysteresis and ex-ternal activation of Spo0A ∼ P and SigH, α SH . (C) Curve of ˆ Q − ˆ Q c ( ˆ P ) vs. α A . (D) Hysteresis andexternal activation of AbrB, α A .tric circles. When the pH reduces to less than 6.5, concentric circular colonies thatexpand periodically are formed. When the pH is less than 5.3, colony formationbecomes unstable and stops halfway, or no colony is formed. Moreover, at a pHless than 5.1, colonies were never formed.To consider the environmental pH sensitivity, the classification parameters D ,ˆ P , ˆ Q and ˆ Q c are associated with S , H , A and C , respectively. The sign of D = p − q is equivalent to whether α A = q / p is less than or greater than 1. This parameter α A can be expressed as α A = b H a H × c A / a A d A = (rate of suppression of H by A ) × (basic production rate of A )(decay rate of A )6hich indicates the rate at which A , AbrB, is activated from outside the modelsystem of the cell type regulation (1). Similarly, we can define the rate at which S and H , Spo0A ∼ P and SigH, are activated from outside the model system as α SH = b A a A × c S / a S d S × c H / a H d H = (rate of suppression of A by S ) × (basic production rate of S activated by H )(decay rate of S ) × (basic production rate of H )(decay rate of H ) . The classification parameters (2) can be expressed asˆ P = α A ( α A − ) , ˆ Q = (cid:18) α SH α A − − (cid:19) , (4)and the following holds from Theorem 1. Corollary 1. If α A ≤ , the set of steady states is anhysteretic. If α A > , thereexists α SH , c ( α A ) such that If α SH ≤ α SH , c ( α A ) , the set of steady states is anhysteretic. If α SH > α SH , c ( α A ) , the set of steady states is hysteretic. According to Corollary 1, the presence or absence of hysteresis can be con-trolled by Spo0A ∼ P or SigH when AbrB is functioning to a certain extent (Fig. 3B).Therefore, the question is whether the activity of Spo0A ∼ P or SigH controlsthe environmental pH-dependent cell type hysteresis. It is known that SigH ispositively regulated with an increase in the environmental pH (Cosby and Zuber,1997; Wilks et al., 2009). Therefore, according to Corollary 1, if SigH is the inputpoint for pH-dependent control, the hysteresis disappears at a low pH (gray dot-ted arrow in Fig. 3B). However, this finding contradicts the pH-dependent controlphenomenon in actual colony observation (left part in Fig. 3A).Note that AbrB has a stronger pH dependence than SigH (Wilks et al., 2009).Moreover, AbrB is upregulated as the environmental pH decreases in this range.Therefore, we consider the possibility that the environmental-pH-dependent AbrBactivity controls the cell type hysteresis. To examine this aspect, we consider the7apping α A ˆ Q − ˆ Q c ( ˆ P ) ( α A ∈ ( , ∞ )) . If α A →
1, then ˆ P → ∞ , and it followsfrom Theorem 1 thatˆ Q − ˆ Q c ( ˆ P ) ∼ (cid:18) α SH α A − − (cid:19) − α A ( α A − ) → − ∞ . Similarly, we see that ˆ P → + α A → ∞ . Theorem 1 indicates thatˆ Q − ˆ Q c ( ˆ P ) ∼ (cid:18) α SH α A − − (cid:19) − α / A ( α A − ) / → − . Furthermore, according to this theorem ∂∂α A ˆ Q ′ c ( ˆ P ) = − ˆ Q ′′ c ( ˆ P ) α A + ( α A − ) < . Therefore, when α A ∈ ( , ∞ ) increases, ˆ Q ′ c ( ˆ P ) is positive and decreases monoton-ically. Consequently, ∂∂α A (cid:0) ˆ Q − ˆ Q c ( ˆ P ) (cid:1) = ( α A − ) (cid:26)(cid:18) + α A − (cid:19) ˆ Q ′ c ( ˆ P ) − α SH (cid:27) changes its sign only once from positive to negative. In other words, there exists an α A , M ∈ ( , ∞ ) such that α A ˆ Q − ˆ Q c ( ˆ P ) increases when α A < α A , M and decreaseswhen α A > α A , M (Fig. 3C).Accordingly, Theorem 1 can be expressed as follows. Corollary 2.
There exists an α ∗ SH > such that if α SH ≤ α ∗ SH , the set of steady states is hysteretic; if α SH > α ∗ SH , < α A , L < α A , U such that (i) if < α A ≤ α A , L , the set of steady states is anhysteretic. (ii) if α A , L < α A < α A , U , the set of steady states is hysteretic. (iii) if α A , U ≤ α A , the set of steady states is anhysteretic. This indicates that when α SH is sufficiently large, that is, when Spo0A ∼ P andSigH function to a reasonable extent, α A or AbrB can control the hysteresis ofthe cell type selection (Fig. 3D). Moreover, the cell type control is (i) anhysteretic(matrix producers dominant) when α A is small; (ii) hysteretic (periodic) when α A is intermediate; (iii) anhysteretic (motile cells dominant) when α A is large. In8erms of the effect of the environmental pH, α A exhibits a negative correlationin the neutral range, that is, AbrB is upregulated at a low pH. In this case, thechange in the cell type control from anhysteretic (matrix producer) to hysteric(periodic) with decreasing pH (Fig. 3A) can be explained by the upregulation of α A (Fig. 3D). Therefore, it is suggested that AbrB plays a central role in cell typeregulation in response to environmental pH changes. In colony formation, it isexpected that the activity of AbrB is the key to the selection of the concentricpattern.
3. Discussion
Flexible and stable control of the cell state according to the environmental con-ditions is the basis for realizing a robust life system. In general, hysteretic controlis one of the methods that exhibits a prompt and stable response to environmentalchanges. In this work, we clarified the necessary and sufficient conditions for thehysteretic control of the cell state, considering the case of the bacterial cell typeregulation as an example. By incorporating this cell type regulation model intoa model of the cell population dynamics (Ben-Jacob et al., 1994; Wakita et al.,1994; Kitsunezaki, 1997; Kawasaki et al., 1997; Golding et al., 1998; Mimura et al.,2000; Tasaki et al., 2017b), the colony morphology can be predicted under a widerange of environmental conditions. In particular, it is possible to correctly repro-duce the formation of concentric colonies that expand periodically, which has notbeen realized so far (Mimura et al., 2000). Moreover, the effect of the environ-mental conditions can be compared with the experimental findings.In addition, we developed a model of the cell type regulation influenced bythe environmental pH changes, and the findings were noted to be consistent withthose of the concentric colony formation experiment. In the existing studies, theconcentration of agar in the medium, which is a control parameter of the cellpopulation motility, has been widely examined as an environmental factor. It isexpected that the growth dynamics of concentric circle colonies depending on theagar concentration can be discussed in combination with the model of the cellpopulation dynamics.Furthermore, the structure of the hysteresis for external signals presented hereinis not limited to
B. subtilis cell type selection. Specifically, a similar structure oc-curs for phenomena that can be expressed in the form of model (1). Furthermore,similar properties are expected if the regulatory network is similar to that shownin Fig. 1A. In addition, as indicated previously, the main result, Theorem 1, can beexpressed in terms of only the indices α A and α SH , which represent the activation9ate parameters from outside the system, as in Corollaries 1 and 2, respectively.Such expressions are universal and can help elucidate the regulatory mechanismsof cell populations. This analysis methodology, which focuses on the presenceor absence of hysteresis related to external signals, can provide a basis for thecomprehensive understanding of other control systems, among other applications. Acknowledgments
This work was supported by JSPS KAKENHI, Grant Number 19K03645 (S.T.),and MEXT KAKENHI, Grant Number 17H06327 (S.T.).
Appendix A. Hysteretic and anhysteretic curves: Proof of Theorem 1
Herein, we present the proof of Theorem 1. The equation for the steady state(equilibrium point) ( S , H , A , C ) of (1) is as follows: = c S Ha S + b S C − d S S = c H a H + b H A − d H H = c A a A + b A S − d A A = c C X A − d C C (A.1)If we eliminate H , A and C using these expressions, X in S can be defined as X = λσ ( σ − p ) ( σ + q ) − µσ = − σ ( σ − p ) ( σ + q ) (cid:2) µσ − { λ + µ ( p − q ) } σ − µ pq (cid:3) (A.2)where σ = a A + b A S , p = a A , q = b H c A a H d A , λ = c S c H b A d A d C b S d S a H d H c A c C , µ = a S d A d C b S c A c C . (A.3)Furthermore, by setting µσ − { λ + µ ( p − q ) } σ − µ pq = µ ( σ + a ) ( σ − b ) , (A.4)10 can be expressed as X = − µσ ( σ + a ) ( σ − b )( σ − p ) ( σ + q ) , (A.5)where, κ = λµ = c S c H b A a S d S a H d H , − a = (cid:26) κ + p − q − q ( κ + p − q ) + pq (cid:27) < , b = (cid:26) κ + p − q + q ( κ + p − q ) + pq (cid:27) > . (A.6)In addition, the range of σ for which X > p < σ < b . Differentiating (A.5)with respect to σ yields dXd σ = µ ( σ − p ) ( σ + q ) ϕ ( σ ) , (A.7)where ϕ ( σ ) = σ ( σ + a ) ( σ − b ) ( σ + q − p ) − ( σ − p ) ( σ + q ) { ( σ + a ) ( σ − b ) + σ ( σ − b ) + σ ( σ + a ) } . (A.8)This indicates that ϕ ( p ) = p ( p + a )( p − b )( p + q ) < , ϕ ( b ) = − b ( b − p )( b + q )( b + a ) < . (A.9)According to this expression, and because ϕ ′ ( p ) = − p κ <
0, the quartic func-tion ϕ = ϕ ( σ ) can be classified into the following two types:1. ϕ ( σ ) ≤ p < σ < b .2. There exist σ − and σ + such that p < σ − < σ + < b , ϕ ( σ ) > σ − < σ < σ + , ϕ ( σ ) ≤ σ ≤ σ − or σ ≥ σ + ).Each case corresponds to an anhysteretic (Fig. 1B) and hysteretic case (Fig. 1C).Expanding (A.8) yields ϕ ( σ ) = − σ + ( p − q ) σ + { pq + ( p − q )( a − b ) − ab } σ + ( a − b ) pq σ − abpq . (A.10)11ecause b − a = p − q + κ and ab = pq , ϕ ( σ ) = − σ + ( p − q ) σ + { pq − ( p − q )( p − q + κ ) } σ − ( p − q + κ ) pq σ − p q . (A.11)Furthermore, because D = p − q and P = pq , ϕ ( σ ) = − σ + D σ + { P − D ( D + κ ) } σ − ( D + κ ) P σ − P . (A.12)As p < σ < b , the range of the variable σ is12 (cid:16) D + p D + P (cid:17) < σ < (cid:18) D + κ + q ( D + κ ) + P (cid:19) . (A.13)This treatment indicates that the hysteretic or anhysteretic nature of the set ofsteady states depends on whether the function ϕ ( σ ) , defined by (A.12) in therange (A.13), does or does not (non-positive) undergo a sign change, respectively.In the following text, we examine the cases of D ≥ D < Appendix A.1. No hysteresis under low influence of AbrB
First, we demonstrate that hysteresis does not occur when AbrB has a lowinfluence, that is, D ≥ α A ≤ D = α A = √ P < σ < (cid:16) κ + p κ + P (cid:17) , (A.14)and ϕ ( σ ) = − (cid:0) σ − P (cid:1) − κ P σ < . (A.15)When D > α A < P , κ , σ by D = − η so that P = η ˆ P , κ = η ˆ κ , σ = η ˆ σ , then η <
0, ˆ P >
0, ˆ κ <
0, ˆ σ <
0, and (A.13) becomes12 (cid:18) − η + q η + η ˆ P (cid:19) < η ˆ σ < (cid:26) − η + η ˆ κ + q ( − η + η ˆ κ ) + η ˆ P (cid:27) . (A.16)Dividing this expression by η ( < ) yields12 (cid:26) − + ˆ κ − q ( − ˆ κ ) + P (cid:27) < ˆ σ < (cid:16) − − p + P (cid:17) . (A.17)By setting ˆ Q = ( ˆ κ − ) < − /
2, the domain can be expressed asˆ Q − q ˆ Q + P < ˆ σ < (cid:16) − − p + P (cid:17) . (A.18)12urthermore, by considering h ( z ) = z − z − (cid:0) ˆ P + ˆ Q (cid:1) z − P ˆ Qz + ˆ P , (A.19)we can obtain ϕ ( σ ) = − η h ( − ˆ σ ) . As p < σ < b , the range of the variable z = − ˆ σ is 12 (cid:16) + p + P (cid:17) < z < − ˆ Q + q ˆ Q + P . (A.20)Therefore, it is sufficient to show that the function h ( z ) , defined by (A.19) in(A.20), is positive. From − Q >
1, it follows that h ( z ) = z (cid:8) z − z + (cid:0) − Q (cid:1)(cid:9) + Pz (cid:8) − z + (cid:0) − Q (cid:1)(cid:9) + ˆ P > z (cid:0) z − z + (cid:1) + Pz ( − z + ) + ˆ P = z ( z − ) − Pz ( z − ) + ˆ P = (cid:0) z − z − ˆ P (cid:1) > . (A.21)Therefore, no hysteresis occurs. Appendix A.2. Hysteresis condition under high influence of AbrB
Herein, we consider in detail the conditions under which the steady state be-comes hysteretic when AbrB has a strong influence, that is, when D < α A > D = − η , P = η ˆ P , κ = η ˆ κ , σ = η ˆ σ yields ϕ ( σ ) = η (cid:8) − ˆ σ + σ + (cid:0) ˆ P + ˆ Q (cid:1) ˆ σ − P ˆ Q ˆ σ − ˆ P (cid:9) , (A.22)where ˆ Q = ( ˆ κ − ) > − /
2. Hence, considering h ( z ) = − z + z + (cid:0) ˆ P + ˆ Q (cid:1) z − P ˆ Qz − ˆ P (A.23)yields ϕ ( σ ) = − η h ( ˆ σ ) . Because p < σ < b , the range of the variable z = ˆ σ is12 (cid:16)p + P − (cid:17) < z < q ˆ Q + ˆ P + ˆ Q . (A.24)Thus, the hysteretic or anhysteretic nature of the set of steady states depends onwhether the function h ( z ) , defined as in (A.23) in the domain (A.24), does or doesnot (non-negative) undergo a sign change.Because h ( z ) = Qz (cid:0) P − z (cid:1) + (cid:0) z − ˆ P (cid:1) + z , (A.25)13 ( z ) > z ≤ P . In particular, h has a fixed value independent of ˆ Q at z = P : h ( P ) = ˆ P ( P + ) . Conversely, in view of h ( z ) < ⇔ z + z − Pz + ˆ P < z (cid:0) z − P (cid:1) ˆ Q , (A.26)for each z > P , we see that h ( z ) < ⇔ ˆ Q > (cid:0) z − ˆ P (cid:1) + z z (cid:0) z − P (cid:1) . (A.27)Hence, h ( z ) is negative if ˆ Q is sufficiently large, and the set of steady states ishysteretic. In contrast, if ˆ Q ≤ h ( z ) is always non-negative, and the set of steadystates is anhysteretic. Furthermore, ∂ h ∂ ˆ Q = − z (cid:0) z − P (cid:1) (A.28)is negative if z > P . Therefore, for each ˆ P , there exists a threshold ˆ Q c ( ˆ P ) > Q ≤ ˆ Q c ( ˆ P ) , the set of steady states is anhysteretic;2. if ˆ Q > ˆ Q c ( ˆ P ) , the set of steady states is hysteretic. Appendix A.3. Profile of the hysteresis threshold curve
We present the profile (3) of the threshold function ˆ Q c = ˆ Q c ( ˆ P ) to completethe proof of Theorem 1. First, the threshold curve is characterized as (cid:0) z , ˆ Q (cid:1) = (cid:0) z c ( ˆ P ) , ˆ Q c ( ˆ P ) (cid:1) , which satisfies the following two equations for ˆ P : F (cid:0) z , ˆ P , ˆ Q (cid:1) = z + z − (cid:0) ˆ P + ˆ Q (cid:1) z + P ˆ Qz + ˆ P = , (A.29) F z (cid:0) z , ˆ P , ˆ Q (cid:1) = z + z − (cid:0) ˆ P + ˆ Q (cid:1) z + P ˆ Q = . (A.30)Then, (3) follows from the following three lemmas. Lemma 1. ˆ Q ′ c ( ˆ P ) = − F ˆ P F ˆ Q > , ˆ Q ′′ c ( ˆ P ) = − F ˆ P ˆ P F ˆ Q − F ˆ P F ˆ Q ˆ P F Q > . (A.31) Here, (cid:0) z c ( ˆ P ) , ˆ P , ˆ Q c ( ˆ P ) (cid:1) was omitted on the right side, for example, F ˆ P = F ˆ P (cid:0) z c ( ˆ P ) , ˆ P , ˆ Q c ( ˆ P ) (cid:1) . emma 2. As ˆ P → , the following holds.z c ( ˆ P ) ˆ P / → , ˆ Q c ( ˆ P ) ˆ P / → . (A.32) Lemma 3. As ˆ P → ∞ , the following holds.z c ( ˆ P ) ˆ P → , ˆ Q c ( ˆ P ) ˆ P → . (A.33) Proof of Lemma 1.
By differentiating F (cid:0) z c ( ˆ P ) , ˆ P , ˆ Q c ( ˆ P ) (cid:1) = P ,(A.31) can be obtained. In particular, we can show that F ˆ P = − z + Qz + P > , F ˆ P ˆ P = > , (A.34) F ˆ Q = − z (cid:0) z − P (cid:1) < , F ˆ Q ˆ P = z > (cid:0) z , ˆ Q (cid:1) = (cid:0) z c ( ˆ P ) , ˆ Q c ( ˆ P ) (cid:1) . z − P > z − Qz − ˆ P <
0, as follows.Assume that z ≥ ˆ Q . According to (A.24), z − ˆ Q < p ˆ Q + ˆ P . If both sidesare squared, we obtain (cid:0) z − ˆ Q (cid:1) < ˆ Q + ˆ P as z − ˆ Q ≥
0. Rearranging this ex-pression yields z − Qz − ˆ P <
0. However, when z < ˆ Q , we have ˆ Q > z − Qz − ˆ P = − z (cid:0) ˆ Q − z (cid:1) − ˆ Qz − ˆ P < Proof of Lemma 2.
Substituting z = α ˆ P / into (A.29) and (A.30) yieldsˆ P / (cid:16) α ˆ P / + α ˆ P / − α ˆ P − α ˆ Q + α ˆ P / ˆ Q + ˆ P / (cid:17) = , (A.36)2 ˆ P / (cid:16) α ˆ P / + α ˆ P / − α ˆ P − α ˆ Q + P / ˆ Q (cid:17) = , (A.37)respectively. Because ˆ P > α > Q = α ˆ P / + α ˆ P / − α ˆ P + ˆ P / α (cid:0) α − P / (cid:1) , (A.38)ˆ Q = α ˆ P / + α ˆ P / − α ˆ P (cid:0) α − ˆ P / (cid:1) . (A.39)In this case, we assumed that α > P / , which is true for small ˆ P if α → α > P →
0. Therefore, α = α ˆ P / + α ˆ P / − α ˆ P + ˆ P / α ˆ P / + α ˆ P / − α ˆ P · α − ˆ P / α − P / . (A.40)15f we assume that α → α as ˆ P →
0, then α = α + α . (A.41)Hence α =
1, and (A.32) can be obtained.The above arguments are based on the assumption that there exists a limit of α ( ˆ P ) . The implicit function theorem can be used to justify these arguments. Proof of Lemma 3.
Substituting z = α ˆ P into (A.29) and (A.30) yieldsˆ P (cid:0) α ˆ P + α ˆ P − α ˆ P − α ˆ Q + α ˆ Q + (cid:1) = , (A.42)2 α ˆ P + α ˆ P − α ˆ P − ( α − ) ˆ P ˆ Q = , (A.43)respectively. Because ˆ P > α > α ˆ P + α ˆ P − α ˆ P + = α ( α − ) ˆ Q , (A.44)2 α ˆ P + α ( α − ) ˆ P = ( α − ) ˆ Q . (A.45)If we assume that α > Q = α (cid:8) α ˆ P + ( α − ) ˆ P (cid:9) + α ( α − ) , (A.46)ˆ Q = α (cid:8) α ˆ P + ( α − ) ˆ P (cid:9) ( α − ) . (A.47)Therefore, α ( α − ) α − = α (cid:8) α ˆ P + ( α − ) ˆ P (cid:9) + α (cid:8) α ˆ P + ( α − ) ˆ P (cid:9) . (A.48)If we assume that α → α ∞ as ˆ P → ∞ , then α ∞ ( α ∞ − ) α ∞ − = α ∞ . (A.49)Since we have assumed that α >
2, necessarily we have α ∞ ≥
2. In this case, α ∞ =
3, and (A.33) can be obtained.Similar to Lemma 2, the implicit function theorem can be used to justify thesearguments. 16 ppendix B. Stability analysis
We present the proofs pertaining to the stability and instability of the equilib-rium points. The stability of an equilibrium point ( S , H , A , C ) in (1) indicates thatthe real parts of all the roots of the following characteristic equation are negative:det − d S − λ p S − n S − d H − λ − n H − n A − d A − λ
00 0 p C − d C − λ = , (B.1)where n S = b S c S H ( a S + b S C ) , n H = b H c H ( a H + b H A ) , n A = b A c A ( a A + b A S ) , (B.2) p S = c S a S + b S C , p C = c C X (B.3)are all positive constants (except p C = X = λ + a λ + a λ + a λ + a = , (B.4)the coefficients are as follows: a = d S d H d A d C − ( p S n A n H d C + p C n A n S d H ) , a = d S d H ( d A + d C ) + d A d C ( d S + d H ) − ( p S n A n H + p C n A n S ) , a = d S d H + d A d C + ( d S + d H ) ( d A + d C ) , a = d S + d H + d A + d C . (B.5)Then, according to the Routh–Hurwitz criterion, the necessary and sufficient con-dition for the real parts of λ , λ , λ , and λ to be negative is for the following sixinequalities to be satisfied: a > a > a > a > (cid:18) a a a (cid:19) = a a − a > , det a a a a a a = a a a − a a − a > . Among these relations, a > a > a a − a > emma 4. If a ≥ , then a a a − a a − a > . Overall, only two conditions remain for stability: a > a >
0. On this basis,we can state that1. Any equilibrium point is stable when the equilibrium curve X = X ( S ) isanhysteretic (Fig. 1B).2. When the curve of the equilibrium point is hysteretic, the equilibrium pointbetween the two folding points (point at X ′ ( S ) =
0) is unstable, and the outerequilibrium point is stable (Fig. 1C).To prove this aspect, it is sufficient to demonstrate the following four lemmas.
Lemma 5. sgn λ λ λ λ = − sgn ϕ ( σ ) . Here, ϕ ( σ ) ( σ = a A + b A S ) is the function defined by (A.8) representing the in-crease/decrease in the curve X = X ( S ) of the equilibrium points, and sgn X ′ ( S ) = sgn ϕ ( σ ) . In other words, according to Lemma 5, the zero eigenvalue appearsat the turning point of the curve. Moreover, the sign of a = λ λ λ λ changesthrough the turning point. Specifically, if the system is stable from a certain pointto the turning point, it becomes unstable at the turning point. Lemma 6.
If the real parts of all the eigenvalues λ , λ , λ , and λ are non-positive, none of them are pure imaginary numbers. This Lemma indicates that when the stability changes along the curve of theequilibrium points, it must pass through the zero eigenvalue and not the pure imag-inary number. Moreover, according to Lemma 5, the point at which the stabilitychanges is the turning point of the curve of the equilibrium points.
Lemma 7.
If X is sufficiently large, any equilibrium point is stable.
Lemma 8.
If X = , any equilibrium point is stable. Thus, the equilibrium point is always stable if no folding points occur. When thereare two turning points, all equilibrium points between them are unstable, and theother equilibrium points are stable. 18 roof of Lemma 4.
By substituting ˜ p S = p S n A n H and ˜ p C = p C n A n S , we have a a a − a a − a = (cid:8) ( d C + d A ) d H + (cid:0) d C + d A d C + d A (cid:1) d H + d A d C + d A d C (cid:9) d S + (cid:8) − ( d H + d A ) ˜ p S − ( d C + d A ) ˜ p C + ( d C + d A ) d H + (cid:0) d C + d A d C + d A (cid:1) d H + (cid:0) d C + d A d C + d A d C + d A (cid:1) d H + d A d C + d A d C + d A d C (cid:9) d S + h(cid:8) − d H + ( d C − d A ) d H + d C + d A d C − d A (cid:9) ˜ p S + (cid:8) d H + ( d A + d C ) d H − d A − d A d C − d C (cid:9) ˜ p C + (cid:0) d C + d A d C + d A (cid:1) d H + (cid:0) d C + d A d C + d A d C + d A (cid:1) d H + (cid:0) d A d C + d A d C + d A d C (cid:1) d H + d C d A + d C d A i d S − ˜ p S + (cid:8) − p C − d A d H + (cid:0) d C + d A d C − d A (cid:1) d H + d C + d A d C (cid:9) ˜ p S − ˜ p C + (cid:8) d H + ( d C + d A ) d H + d A d C d H − d A d C − d A d C (cid:9) ˜ p C + (cid:0) d A d C + d A d C (cid:1) d H + (cid:0) d A d C + d A d C + d A d C (cid:1) d H + (cid:0) d A d C + d A d C (cid:1) d H . As a ≥
0, ˜ p S ≤ d S d H d A , ˜ p C ≤ d S d A d C . Using these expressions for the negativeterms in the abovementioned equation, we can obtain the following expression: a a a − a a − a ≥ (cid:8) d C d H + (cid:0) d C + d A d C (cid:1) d H (cid:9) d S + (cid:8) d C d H + (cid:0) d C + d A d C (cid:1) d H + (cid:0) d C + d A d C + d A d C (cid:1) d H (cid:9) d S + h(cid:8) d C d H + d C + d A d C (cid:9) ˜ p S + (cid:8) d H + ( d A + d C ) d H (cid:9) ˜ p C + (cid:0) d C + d A d C (cid:1) d H + (cid:0) d C + d A d C + d A d C (cid:1) d H + (cid:0) d A d C + d A d C + d A d C (cid:1) d H i d S + (cid:8)(cid:0) d C + d A d C (cid:1) d H + d C + d A d C (cid:9) ˜ p S + (cid:8) d H + ( d C + d A ) d H + d A d C d H (cid:9) ˜ p C + (cid:0) d A d C + d A d C (cid:1) d H + (cid:0) d A d C + d A d C + d A d C (cid:1) d H + (cid:0) d A d C + d A d C (cid:1) d H , the right-hand side of which is positive. Proof of Lemma 5.
By using the equations of the equilibrium point, (A.1), in(B.5) to eliminate H , A , C and X , and by rearranging the equation by using σ (= a A + b A S ) , α A , α SH , we can obtain a = − d S d H d A d C a A σ ( σ + a A α A ) α SH ϕ ( σ ) . (B.6)19herefore, from a = λ λ λ λ , the claim is verified. Proof of Lemma 6.
The necessary and sufficient condition for the characteristicequations (B.4) to have a pure imaginary root λ = ρ i ( ρ = ) is a a a − a a − a = ρ − a ρ + a = ρ (cid:0) a ρ − a (cid:1) =
0. According to the as-sumption, because the Routh–Hurwitz conditions hold with equal signs, we have a ≥
0. Therefore, according to Lemma 4, a a a − a a − a >
0. Thus, thecharacteristic equation (B.4) does not have any pure imaginary root.
Proof of Lemma 7.
As mentioned previously, only two of the Routh–Hurwitz con-ditions must be examined: a > a >
0. When X is sufficiently large, S ∼ ( σ ∼ a A ) . Then, the equilibrium point satisfies X ′ ( S ) <
0, and therefore, ϕ ( σ ) <
0. Thus, a > a ,in the same manner as (B.6), by rearranging the expression using σ , α A and α SH ,we can obtain a = d S d H ( d A + d C ) + d A d C ( d S + d H ) − d S d H d A a A α A ( σ − a A ) σ ( σ + a A α A ) − d A d C d S σ − a A σ (cid:26) − ( σ − a A ) ( σ + a A α A ) α SH a A σ (cid:27) . (B.7)Because σ ∼ a A , we see that a ∼ d S d H ( d A + d C ) + d A d C ( d S + d H ) > Proof of Lemma 8.
As in the previous case, it is sufficient to consider only twoconditions: a > a >
0. When X =
0, we have σ = b . Then, the equilibriumpoint satisfies X ′ ( S ) <
0, and therefore, ϕ ( σ ) <
0. Thus, a > σ = b using α A and α SH , we obtain b = η (cid:18)q ˆ Q + ˆ P + ˆ Q (cid:19) = p (cid:18) qp − (cid:19) ( s(cid:18) α SH − α A + α A − (cid:19) + α A ( α A − ) + (cid:18) α SH − α A + α A − (cid:19)) = a A ( α A − ) ( s(cid:18) α SH − α A + α A − (cid:19) + α A ( α A − ) + (cid:18) α SH − α A + α A − (cid:19)) = : a A ( α A − ) b ∗ . (B.8)20ere, we can see that b ∗ satisfies the following condition b ∗ ≥ α A − , b ∗ ≥ α SH − α A + α A − b > p = a A . Then, from p S = c S a S , n H = b H c H a H + b H c A d A σ ! , n A = b A c A σ , (B.10)we have p S n H n A d A d S d H = α SH α A σ a A + α A ! = α SH α A { b ∗ ( α A − ) + α A } . (B.11)By using the two inequalities (B.9) at the square of the denominator on the right-hand side, one can obtain p S n H n A d A d S d H ≤ α SH + α SH α A + α A < . (B.12)Moreover, recall that p C = X =
0. In this case, a = d S d H d C + d A d C ( d S + d H ) + d S d H d A (cid:18) − p S n H n A d A d S d H (cid:19) > d S d H d C + d A d C ( d S + d H ) , (B.13)and, a > References
Ben-Jacob, E., Schochet, O., Tenenbaum, A., Cohen, I., Czirok, A., Vicsek, T.,1994. Generic modelling of cooperative growth patterns in bacterial colonies.Nature 368, 46–9. doi: .Branda, S.S., Vik, ˚A., Friedman, L., Kolter, R., 2005. Biofilms:the matrix revisited. Trends in Microbiology 13, 20–26.doi: https://doi.org/10.1016/j.tim.2004.11.006 .21airns, L.S., Hobley, L., Stanley-Wall, N.R., 2014. Biofilm formation by
Bacillussubtilis : new insights into regulatory strategies and assembly mechanisms. MolMicrobiol 93, 587–98. doi: .Chai, Y., Chu, F., Kolter, R., Losick, R., 2008. Bistability andbiofilm formation in
Bacillus subtilis . Mol Microbiol 67, 254–63.doi: .Cosby, W.M., Zuber, P., 1997. Regulation of
Bacillus subtilis σ H (Spo0H) andAbrB in response to changes in external pH. J Bacteriol 179, 6778–87.Donlan, R.M., 2002. Biofilms: microbial life on surfaces. Emerging infectiousdiseases 8, 881–890. doi: .Flemming, H.C., Wingender, J., Szewzyk, U., Steinberg, P., Rice, S.A., Kjelle-berg, S., 2016. Biofilms: an emergent form of bacterial life. Nature ReviewsMicrobiology 14, 563–575. doi: .Fujikawa, H., 1992. Periodic growth of Bacillus subtilis colonies on agarplates. Physica A: Statistical Mechanics and its Applications 189, 15–21.doi: http://dx.doi.org/10.1016/0378-4371(92)90123-8 .Golding, I., Kozlovsky, Y., Cohen, I., Ben-Jacob, E., 1998. Studies of bac-terial branching growth using reactiondiffusion models for colonial develop-ment. Physica A: Statistical Mechanics and its Applications 260, 510–554.doi: http://dx.doi.org/10.1016/S0378-4371(98)00345-8 .Hall-Stoodley, L., Costerton, J.W., Stoodley, P., 2004. Bacterial biofilms: fromthe natural environment to infectious diseases. Nat Rev Microbiol 2, 95–108.doi: .Hobley, L., Harkins, C., MacPhee, C.E., Stanley-Wall, N.R., 2015. Giv-ing structure to the biofilm matrix: an overview of individual strate-gies and emerging common themes. FEMS Microbiol Rev 39, 649–69.doi: .Itoh, H., Wakita, J.i., Matsuyama, T., Matsushita, M., 1999. Periodic patternformation of bacterial colonies. Journal of the Physical Society of Japan 68,1436–1443. doi: .22awasaki, K., Mochizuki, A., Matsushita, M., Umeda, T., Shigesada, N., 1997.Modeling spatio-temporal patterns generated by
Bacillus subtilis . J Theor Biol188, 177–85. doi: .Kearns, D.B., Losick, R., 2005. Cell population heterogeneity during growth of
Bacillus subtilis . Genes Dev 19, 3083–94. doi: .Kitsunezaki, S., 1997. Interface dynamics for bacterial colony for-mation. Journal of the Physical Society of Japan 66, 1544–1550.doi: .Kobayashi, K., 2007. Gradual activation of the response regulatorDegU controls serial expression of genes for flagellum formation andbiofilm formation in
Bacillus subtilis . Mol Microbiol 66, 395–409.doi: .Kobayashi, K., Iwano, M., 2012. BslA(YuaB) forms a hydrophobic layeron the surface of
Bacillus subtilis biofilms. Mol Microbiol 85, 51–66.doi: .L´opez, D., Vlamakis, H., Kolter, R., 2009. Generation of multiplecell types in
Bacillus subtilis . FEMS Microbiol Rev 33, 152–63.doi: .Mimura, M., Sakaguchi, H., Matsushita, M., 2000. Reac-tiondiffusion modelling of bacterial colony patterns. Phys-ica A: Statistical Mechanics and its Applications 282, 283–303.doi: http://dx.doi.org/10.1016/S0378-4371(00)00085-6 .Shimada, H., Ikeda, T., Wakita, J.i., Itoh, H., Kurosu, S., Hiramatsu, F., Nakat-suchi, M., Yamazaki, Y., Matsuyama, T., Matsushita, M., 2004. Dependenceof local cell density on concentric ring colony formation by bacterial species
Bacillus subtilis . Journal of the Physical Society of Japan 73, 1082–1089.doi: .Shioi, J.I., Matsuura, S., Imae, Y., 1980. Quantitative measurements of protonmotive force and motility in
Bacillus subtilis . J Bacteriol 144, 891–7.Stoodley, P., Sauer, K., Davies, D.G., Costerton, J.W., 2002. Biofilms as com-plex differentiated communities. Annual Review of Microbiology 56, 187–209.doi: .23asaki, S., Nakayama, M., Shoji, W., 2017a. Morphologies of
Bacillus subtilis communities responding to environmental variation. Dev Growth Differ 59,369–378. doi: .Tasaki, S., Nakayama, M., Shoji, W., 2017b. Self-organization of bacterial com-munities against environmental pH variation: Controlled chemotactic motil-ity arranges cell population structures in biofilms. PLoS One 12, e0173195.doi: .Tasaki, S., Nakayama, M., Takagi, I., Shoji, W., 2020. Life cycle of
Bacillussubtilis population: biofilm life cycle generation and response to environmentalpH. bioRxiv, doi: .Vlamakis, H., Chai, Y., Beauregard, P., Losick, R., Kolter, R., 2013. Stickingtogether: building a biofilm the
Bacillus subtilis way. Nat. Rev. Microbiol. 11,157–168. doi: .Wakita, J.i., Komatsu, K., Nakahara, A., Matsuyama, T., Matsushita, M., 1994.Experimental investigation on the validity of population dynamics approach tobacterial colony formation. Journal of the Physical Society of Japan 63, 1205–1211. doi: .Wakita, J.i., Shimada, H., Itoh, H., Matsuyama, T., Matsushita, M., 2001. Periodiccolony formation by bacterial species
Bacillus subtilis . Journal of the PhysicalSociety of Japan 70, 911–919. doi: .Wilks, J.C., Kitko, R.D., Cleeton, S.H., Lee, G.E., Ugwu, C.S., Jones, B.D.,BonDurant, S.S., Slonczewski, J.L., 2009. Acid and base stress and tran-scriptomic responses in
Bacillus subtilis . Appl Environ Microbiol 75, 981–90.doi: .Yamazaki, Y., Ikeda, T., Shimada, H., Hiramatsu, F., Kobayashi, N., Wakita, J.i.,Itoh, H., Kurosu, S., Nakatsuchi, M., Matsuyama, T., Matsushita, M., 2005.Periodic growth of bacterial colonies. Physica D: Nonlinear Phenomena 205,136–153. doi:10.1016/j.physd.2004.12.013