From Staphylococcus aureus gene regulation to its pattern formation

The focus of this paper is to develop a new partial differential equation model for the pattern formation of the human pathogen Staphylococcus aureus, starting from a newly developed model of selected gene regulation mechanisms. In our model, we do not only account for the bacteria densities and nutrient concentrations, but also for the quorum sensing and biofilm components, since they enable bacteria to coordinate their behavior and provide the environment in which the colony grows. To this end, we model the relevant gene regulation systems using ordinary differential equations and therefrom derive our evolution equations for quorum sensing and biofilm environment by time-scale arguments. Furthermore, we compare and validate our model and the corresponding simulation results with biological real data observations of Staphylococcus aureus mutant colony growth in the laboratory. We show that we are able to adequately display the qualitative biological features of pattern formation in selected mutants, using the parameter changes indicated by the gene regulation mechanisms.


Introduction
The bacterium Staphylococcus aureus (S. aureus) is a Gram-positive, non-sporulating bacterium found in nature and on mammal skin. Its pathogenic form is of specific interest due to its role in hospital-acquired infections from medical devices like prosthetics or surgical instruments (Bronner et al. 2004). In these infections the formation of biofilm is crucial to bacterial survival. Pattern formation in Gram-positive bacteria has been studied regarding biology and mathematical modeling, e.g. by Ben-Jacob et al. (2001), García-Betancur et al. (2017), Matsushita et al. (1999) and Mimura et al. (2000), often for the close S. aureus relative B. subtilis, e.g. by Dervaux et al. (2014), Kawasaki et al. (1997), Seminara et al. (2012) and Vlamakis et al. (2013). In this paper, we develop a new partial differential equation (PDE) model for the pattern formation in S. aureus, including nutrient concentration, replicative and non-replicative bacteria densities, quorum sensing substance and biofilm environment concentration. Therefore, we consider gene regulation mechanisms leading to the production of the representative substances autoinducing peptide (AIP) and polysaccharide intercellular adhesin (PIA) involved in S. aureus quorum sensing and biofilm formation. The complex regulation systems are described by ordinary differential equations (ODE) and the resulting model is subsequently reduced with the aim of using the derived equations in the description of bacterial pattern formation.
In this context, biofilm formation is crucial since in our modeling the biofilm constitutes the environment in which the bacterial colony grows and develops its characteristic features. We note that in our modeling the term biofilm does not denote the entire biomass as in the classical sense, it is only the environment in which the bacteria live and it is considered as a product of the replicative bacteria which does not grow on its own. In addition, quorum sensing (QS) is a cell-cell communication mechanism, which coordinates gene expression and behavior in bacterial populations, depending on population densities (Fuqua et al. 1994;Müller et al. 2006;Perez-Velazquez et al. 2016). In S. aureus, the QS signalling molecule AIP is produced by the agr operon (Boles and Horswill 2008;Yarwood and Schlievert 2003;Yarwood et al. 2004). AIP impedes the attachment and development of a biofilm (Boles and Horswill 2008). An increasing concentration of NaCl in the experimental set up leads to increased biofilm formation and lower growth rates and the availability of nutrients is thus linked to commitment towards biofilm dispersal (Mhatre et al. 2014). During population growth, on top of the biofilm a thin film of liquid is formed.
From a mathematical perspective, general quorum sensing and gene regulation at the agr operon are studied by Ward et al. (2001) and by Jabbari et al. (2010Jabbari et al. ( , 2012a, respectively. In the paper by Jabbari et al. (2010), an ODE model of the agr operon is established and equation systems are asymptotically analyzed on different time scales. A goal in the following is to derive PDE for the QS substance AIP and the biofilm substance PIA and terms for the diffusion coefficients for the nutrients and replicative bacteria. The modeling of the nutrient concentration and bacteria densities is based on the approach for pattern formation in S. aureus by García-Betancur et al. (2017). The rest of the paper is structured as follows: In Sect. 2, we introduce our ODE system for the gene regulation mechanisms in S. aureus and perform the model reduction, which yields the reaction terms for the quorum sensing substance and biofilm equations in the pattern formation PDE model. In Sect. 3, we verify that our model adequately predicts colony shapes for S. aureus mutant colonies by comparing our simulation results to biological observations of colony growth in the laboratory.

Mathematical modeling
In the paper by García-Betancur et al. (2017) a system of four equations was introduced, which takes into account nutrient concentration, replicative and non-replicative bacteria densities and quorum sensing substance concentration. Bacteria in the reproductive state are assumed to move according to a mixed diffusion dispersal term, while bacteria in the non-replicative state only change locations through pushing induced by the movement of the replicative bacteria. However, the model by García-Betancur et al. (2017) includes a generic quorum sensing term and does not take into account biofilm formation and the related effects, which are crucial to colony growth. Our newly developed model is based on the structure as introduced in García-Betancur et al. (2017), but with two major changes: there is a separate variable and equation included for the biofilm component, separate from the two bacterial components. Secondly, the details of the biologically known underlying gene regulation mechanisms are taken into account. With this improved understanding we aim for additional predictions, which associate parameter changes for a certain mutant to gene regulation effects. Thus, in the following we choose the terms unrelated to QS and biofilm formation as introduced by García- Betancur et al. (2017) and aim to determine where n, b, s, q and f denote the nutrient concentration, the replicative and nonreplicative bacteria densities, the AIP concentration and the biofilm environment concentration, which in our case corresponds to the concentration of PIA used as an indicator for biofilm. Note that we assume the diffusion of non-replicative bacteria to be driven by the movement of the replicative cells, which push the non-replicative bacteria in a limited process. In the case b b s , the diffusion of the non-replicative bacteria becomes approximately linear. However, there is only a fixed amount of nutrients in the system and if b b s , these nutrients are consumed very fast. This results in a fast conversion of replicative bacteria to non-replicative bacteria, such that the threshold b s holds again. For a detailed discussion of the remaining diffusion coefficients we refer to Sect. 2.6. Furthermore, nutrient availability influences Degradation and dilution κ Translation colony structure in that nutrient gradients drive colony growth and in that the availability of nutrients determines the transition of bacteria from the replicative to the non-replicative state. Additionally, it should be noted here that the spatial processes have been non-dimensionalized following the procedure introduced in Kawasaki et al. (1997). QS systems regulate the behavior of the bacterial population in response to the density. S. aureus performs QS via the agr operon, which is only found in the Staphylococcus genus (Gray et al. 2013), and the luxS system, which is common for many bacterial species (Xu et al. 2006). The agr system uses the signaling substance AIP. It represses biofilm formation and influences the structure of a biofilm through the regulation factors PSM-α, PSM-β and δ-toxin. The luxS system in S. aureus uses the signalling substance AI-2 and regulates the gene transcription of CP5, a capsular polysaccharide (CP) cell wall component responsible for interaction with the host immune system during the invasive process (Zhao et al. 2010). Since we consider bacterial growth in a laboratory setting, this pathway is not included into the model. Following the paper by Jabbari et al. (2010), we assume that the system is developed for the situation of a well-mixed environment with sufficient supply of ribosomes, such that at each pass the entire strand of mRNA is translated.
We assume that a proportion of the cells is up-regulated for the activity of a regulatory system, such as agr. The gene up-regulation at a locus is increased by the binding of external substances to the promoters. Down-regulation is proportional to the number of up-regulated cells and the gene locus down-regulation rate is increased by the binding of further substances. An overview of the recurring parameters for an example gene locus is found in Table 1. Furthermore, in Fig. 1 an overview of the gene regulation circuits in S. aureus is depicted, showing regulation mechanisms in the agr, the ica, the sarA and the sarA homologue regulation systems (Audretsch et al. 2013;Bronner et al. 2004;Jabbari et al. 2010;Periasamy et al. 2012).

The agr system
The agr system, depicted in the red area in Fig. 1, is considered in detail as an example of a subsystem of gene regulation. It is responsible for the production of the QS substance AIP. In Fig. 1   > 0 and no further down-regulation takes place. In the same way at P agr = 1 no further up-regulation is possible. This equation is of special interest since it represents the regulation of the activity of the agr subsystem, which is involved in the production of AIP. Within the agr subsystem further interactions take place, which are represented in the system of equations displayed in Fig. 2.
Our model of the agr operon depicted in Fig. 2 is based on the model developed by Jabbari et al. (2010), as we consider similar interactions but fewer equations. In the equation for mRNAII (M 2 ) in Fig. 2 we model that each of the N bacteria in the system produces M 2 at the basal production rate m 2 . Increased activity of the agr operon up-regulates this production at the velocity v 2 . Degradation and dilution take place at the rate δ M 2 . Furthermore, we see that the proteins AgrB (B), AgrD (D), AgrC (C) and AgrA (A) are produced by translation of mRNAII. While AgrB is a membrane protein, AgrA can be phosphorylated by AgrC with bound AIP and the QS substance AIP is produced from AgrD under the influence of AgrB. Since the entire strand of mRNA is translated at each pass, we take the translation rates to be of the same order. As in the paper by Jabbari et al. (2010), the variables a, A P and R denote the amounts of AIP, phosphorylated AgrA and AIP-bound AgrC, respectively. The binding of AIP to AgrC is modeled as ± βCa and this binding is assumed to be reversible (± γ R). The phosphorylation and dephosphorylation of AgrA are modeled as φ AR and μA P . Degradation of AIP takes place outside the cells with the corresponding degradation rate λ a . The P3 operon is responsible for the formation of the untranslated RNAIII, an intracellular effector, which up-regulates extracellular protein genes and downregulates cell wall colonisation factor genes (Jabbari et al. 2010). RNAIII experiences transcription but not translation, and can thus base pair with other mRNA strains in order to inhibit encoding of virulence factors (Waters and Storz 2009).
In a similar way, we model the production of the biofilm-promoting hydrophobic phenol soluble modulins (PSM), which are responsible for biofilm structuring (Schwartz et al. 2012) and dispersal and the production of amyloid fibrils (Schwartz et al. 2014). Due to their amphipathic α-helical structure, PSM can lyse eukaryotic cells, such as neutrophils, monocytes, and erythrocytes, by non-specific destruction of biological membranes . The regulation of P P SM takes place by direct binding of the AgrA response regulator (Peschel and Otto 2013;Queck et al. 2008). For P P SM , we obtain the equation While the proportions of up-regulated cells are dimensionless quantities, the dimensional model includes concentrations of molecules per volume unit (cm 3 or mm 3 ) for all proteins and RNAs and the rates are measured per second. We adjust for the additional units by choosing appropriate units for the corresponding rate constants, e.g. the constant k in the production of AIP from AgrD and AgrB has the unit [k] = cm 3 molecules·second . We non-dimensionalize the model for the simulations in order to reduce the number of parameters, using the timescale τ := δ M 2 t for all subsystems, where we assume that the degradation rates of the considered mRNAs are of similar order, i.e.
As an example of the non-dimensionalization we consider the agr subsystem. As in the paper by Jabbari et al. (2010), the nondimensionalization is performed with respect to the starting values calculated as the stationary points for k = 0 and P agr = P P SM = 0. Due to k = 0, we obtain a(0) = R(0) = A P (0) = 0, which allows to choose the non-dimensionalizations for these variables as a := Jabbari et al. (2010), we assume that dilution dominates degradation and thus the δ X are of the same order for all proteins X . We obtain the new parameters λ := δ X δ M 2 and η := δ D δ C ≈ 1. We furthermore reconsider the activation and deactivation factors as well since, due to the new timescale, the variables are rescaled. Thus we define with β agr := N m 3 δ R 3 , which yields b agr := b agr β agr and α := κ X δ X , assumed to be equal for X ∈ {[Sar A], [SarU ], A P }. The obtained non-dimensional equations are depicted in Fig. 3 and, for ease of notation, we omit the primes and directly state the equations for the remaining subsystems in their non-dimensional forms.

The sarA system
We see in Fig. 1 that regulation of the sarA system in the blue area depends on the environmental conditions of the colony, including the temperature T , the osmolarity and the pH level of the environment (Vuong et al. 2005). We summarize these factors as the stress on the system, denoted by [str] ∈ [0, 1], and include it into our model system as one single variable with the optimal reference temperature T ref and pH-value pH ref . In Fig. 1, we observe that the parameter [str] influences transcription of rsbU, rsbV, rsbW and sigB. As described above, the staphylococcal accessory regulator SarA activates transcription at the agr operon. In addition, also S. aureus biofilm formation depends on SarA and σ B (Valle et al. 2003).
We see in Fig. 1 that transcription at the sarA operon is positively regulated by σ B (Cheung and Manna 2005). Furthermore, SarA activates its own expression (Bronner et al. 2004), but this is neglected here for the sake of model simplification. We introduce proportions of up-regulated cells P RsbU for the RsbU-and P σ for the RsbV-, RsbW-and σ B -production as well as the proportion P sar A for the production of SarA. The resulting equations on the timescale τ are where the origin of the factor [C 1 ] is considered in detail in the following. An overview of the non-dimensional equations describing the interactions within the sarA subsystem is displayed in Fig. 4, where the equations are developed analogously to Sect. 2.1. Basal transcription of the genes rsbU, rsbV, rsbW and sigB takes place at the σ Adependent promoter sig B P1 , and the genes rsbV, rsbW and sigB are also transcribed from the σ B -dependent promoter sig B P3 (Senn et al. 2005). Both promoters show a rapid response to environmental stress (Senn et al. 2005). Therefore we distinguish between the variables P RsbU and P σ in (7)-(8). The rsbU gene product, in combination with the rsbV and rsbW gene products, influences the expression of the σ B -operon, responsible for the regulation of sarA transcription, in a competitively inhibitory process. In the regulation of rsbU activity in S. aureus, the RNA polymerase (RNAP) core enzyme and the σ B -factor can associate to form an RNAP holoenzyme (Senn et al. 2005), which recognizes promoter regions in the DNA, initiating transcription at the sig B P3 promoter (Dervaux et al. 2014). If the environmental stress on the system is low, RsbW binds to σ B at rate k 2 , keeping it from aggregating to the RNAP core enzyme and forming the complex C 2 . If the stress increases, RsbU dephosphorylates RsbV, such that RsbV can bind to RsbW in order to release activated σ B -factor, which then binds to the RNAP core enzyme and forms an RNAP holoenzyme (Knobloch et al. 2004). Thus σ B and RsbV compete for the binding with RsbW. The RsbV-RsbW complex C 1 forms at rate k 1 . All bindings occurring at rate k i , i ∈ {1, 2} are reversible at rate k −i . In this framework we formulate the dimensional equations With the stress-dependent rates k 1 := [str]k −1 and k 2 := (1 − [str])k −2 , the binding and unbinding are assumed to be in equilibrium.
Due to mass conservation for σ B , [RsbV ] and [RsbW ], it holds that where σ 0 , v 0 , r 0 ∈ R + are constant parameters. Furthermore, as both rsbV and sig B are transcribed from the P3-operon, we approximate using [RsbV ] ≈ [σ B ]. We obtain Here we concentrate on the dependence on the system stress and take r 0 = 1 to insert (8). We see in Fig. 4 that we assume the amount of dephosphorylated RsbV to be a direct consequence of the availability of RsbU. Furthermore, we denote the binding of RsbW to RsbV with a constant rate as −d[RsbW ] [RsbV ]. In the unstressed as well as in the stressed cell, the ratio of RsbV to phosphorylated RsbV (RsbV-P) is approximately 0.22 if RsbU is present. This ratio does not change significantly due to stress (Pane-Farre et al. 2009), and thus we do not model RsbV-P explicitly. Furthermore, the equations for RsbV, RsbW and σ B are simplified considerably by taking into account that the equations in (10) are assumed to be in equilibrium.

The ica system
The ica locus is a part of the accessory genes found in most S. aureus strains (Cue et al. 2012;Arciola et al. 2015) and consists of the icaR and the icaADBC locus, which are transcribed divergently (Cue et al. 2012). It is of special interest since it is involved in the synthesis of the main extracellular polymeric substance (EPS) component polysaccharide intercellular adhesin (PIA). PIA controls intercellular adhesion during the initial stages of biofilm formation (Cue et al. 2012). Its synthesis is achieved by the combination of all products of the icaADBC gene cluster and modeled as depicted in Fig. 5. While IcaA and IcaD are necessary for the expolysaccharide synthesis, IcaC translocates the polymer poly-N -acetylglucosamine to the bacterial surface and IcaB deacetylates the molecule, enabling fixation to the bacterial cell surface (Cue et al. 2012). The icaADBC cluster genes are therefore modeled together.
Transcription at the ica operon is a direct consequence of the level of Sar A, which binds to the icaADBC promoter region. The influence of σ B is indirect (Cue et al. 2012) since the ica locus does not have a corresponding binding site (Cerca et al. 2008), and thus not considered here. Furthermore, an increasing nutrient concentration inhibits repression of icaADBC by IcaR (Cue et al. 2012). Since the nutrients do not bind to the promoter, we assume logistic up-regulation by nutrients. The IcaR dimers produced at the icaR locus bind in cooperative pairs to the icaADBC locus, where they inhibit transcription (Cue et al. 2012). For the transcription at the icaR locus, we assume that, as in the close S. aureus relative Staphylococcus epidermidis, it is influenced by environmental stress. This is taken into account in the system of equations where we introduce the regulation terms The regulations within the ica system in non-dimensional form are depicted in Fig. 5. Experimental evidence suggests that an intermediate substance exists, whose production is up-regulated by σ B and down-regulated by SarA, and which either degrades PIA or represses PIA synthesis (Valle et al. 2003). We take the additional degradation term δ add [Sar A] for PIA in the dimensional equation setting, which is close to δ add if σ B and SarA are present at similar levels, greater than δ add if there is a higher level of σ B and less than δ add if there is more SarA. The amount of SarA is up-regulated by the amount of σ B , and thus the fraction approaches δ add as time increases, resulting in the additional degradation term δ add in the equation for PIA.

The sarA homologue system
There are at least nine major SarA homologues, including SarR, SarS, Rot, SarT, SarU, SarV, MgrA, SarX and SarZ, where the gene loci sarS and sarT are adjacent, but transcribed divergently . The sarA, sarT and sarU genes are expressed evenly among strains (Ballal and Manna 2009) and thus included into our model. Furthermore the genes mgrA, rot, sarZ, sarR and sarS are expressed differently among strains (Ballal and Manna 2009) and thus only included through the sarS gene as their representative.
The production of protein A and α-toxin through the sarA/agr network (Schmidt et al. 2003) is steered by delicate interactions between the SarA homologues. SarA directly regulates the expression of sarS and spa as well as sarT. The transcription of sarT is further reduced by a high activity of the agr system. SarT binds to the sarS promoter, increasing its activity and SarS influences the transcription of spa. Furthermore, the transcription of hla depends on the level of SarT. SarA binds to the Sar boxes within the promoter regions of the genes encoding protein A and α-toxin (Dunman et al. 2001) in the same way it binds to the promoter region for agr. In addition, mRNAIII regulates the transcription as well as the translation of α-toxin (Bronner et al. 2004). This yields the system of equations An overview of the processes and resulting non-dimensional equations in the sarA homologue system is depicted in Fig. 6.

Model reduction
In this section, we perform the model order reduction, which is aimed at the derivation of ODE for the concentrations of free AIP (a) and PIA ([P I A]). The model reduction includes two main steps: the choice of the timescale of interest, which is the time scale of the change in the proportions of up-regulated cells P act , where act ∈ {agr, P SM, RsbU , σ, sar A, ica1, ica2, sar T , sarU , sar S, spa, hla}, and the derivation of further simplifying parameter size considerations. A recurring structure in the regulation processes is the set of equations for P act , M act and a substance X .
The agr system constitutes a special case, since here mRNAII (M 2 ) and RNAIII (R 3 ) have to be considered. For details we refer to Fig. 3. Since the changes in M act and X take place much faster than the changes in P act , we introduce a parameterε 1. Then the parameters b act and u act are in fact of the form b act =εb act and u act =εû act , whereb act andû act are O(1), and we define the time scale of interest asτ :=ετ to obtain In the limitε → 0 we obtain a reduction of the system by two equations since X = M act = 1 + v act P act , where v act = const.. Then the evolution of X on the slow time scale is a direct consequence of the evolution of P act since We neglect differences among the speeds of those interactions that are faster than d P act dτ since on the slow timescale, for bothε → 0 andε 2 → 0 asymptotics, we approach the steady state. Only those concentrations which provide the connections between the different subsystems are required in an explicit form. Furthermore for effective QS, the induced regulation has to take place faster than the regular effects. Thus there are also significant differences among the remaining rate constants. In agreement with Jabbari et al. (2010), we use ε := m v < 1, the relationship between the basal mRNA transcription rate and the regulation-induced transcription for all subsystems. A possible value is ε = 10 −1 . We obtain for the agr system, where the parametersṽ,β,φ,μ,λ a are O(1).
The connecting concentrations to the agr subsystem are the concentrations a and A P of free AIP and phosphorylated AgrA. We apply the time-scale argument to the set of equations depicted in Fig. 3, calculate the concentration of AIP as and obtain that C a ≈ k D B+γ R β . Inserting this into the expression R = βCa γ +λ , we obtain that R = k D B λ and therefore with B = M 2 that a = k(1 + γ λ ) β (1 + v P agr ).

Furthermore, for A P we obtain that
kφ(1+v P agr ) 2 + 1 and for the concentrations of PSM, amyloid fibrils and mRNAIII we calculate that Using that d X dτ = d X d P act d P act dτ , we obtain the following equations for the concentrations of AIP, PSM, δ-toxin and PIA: Furthermore we use that with the parameters A 1 := k 2 (μ+λ)(λ+γ ) 3 φβ 3 λ 2 , A 2 := k 2 (μ+λ)(λ+γ ) 2 φβ 2 λ 2 and A 3 := k(λ+γ ) βλ . The general solution of (11) is and approaches the constant value exponentially fast. For the inter-cellular regulations to affect the proportions of up-regulated cells in another subsystem, a certain concentration of the substance has to be reached. These concentrations are thus taken as the steady state concentrations, which are obtained from (16)  With A P defined as in (15), the terms in (6) are determined by straight-forward calculations as where Table 2. Thus f ([str]) depends on the system stress in a limited growth process. By inserting the expressions from (15), (17)-(19) and (14), the differential equation for the amount of free AIP is given as with the parameters C i , i ∈ {1, . . . , 11} and B 5 as in Table 3. PSM molecules enable the lysis of cells to gain nutrients. Due to abundant nutrient supply in the medium and a homogeneous bacteria population, cell lysis only plays a minor role and the concentration of PSM is not included explicitly. We remark however that an increase in the PSM concentration is positively influenced by the concentration of free AIP, i.e. agr activity. Furthermore, the equation for δ-toxin only differs by a factor from that for free AIP. It inhibits spreading of the colony, which is used in the replicative bacteria diffusion term of the full model.
Inserting into (14) yields that the concentration of the main biofilm component PIA depends on the evolution of the proportion P ica2 as We insert the expression for [Sar A] from (18) and [str] . Using that v ica2 1, we simplify the equation to obtain with the model parameters P i , i ∈ {1, . . . , 6} described in Table 4.   . 7 Concentrations of AIP and PIA from the full and the reduced ODE system (colour figure online)

Inclusion into the PDE model
In the following, we combine the derived ODE with adequate terms for the spread in space to obtain a PDE model for bacterial pattern formation. The concentration of the S. aureus QS substance AIP is denoted by q in the PDE model. The evolution of the QS substance is obtained by multiplying the evolution of the AIP concentration with the density b of replicative bacteria, which can produce AIP. We derive f 2 (b, q) in (4) from (20) as Biofilm growth is determined by the principal biofilm component PIA, which is included explicitly with the concentration f . Although the production of PIA is not under the control of agr (Vuong et al. 2005), the agr system influences biofilm formation via other mechanisms since repression of agr is necessary to form biofilm and the reactivation of agr in established biofilms through AIP addition or glucose depletion triggers detachment (Boles and Horswill 2008). We include the agr-mediated biofilm detachment depending on the extracellular presence of AIP and f 3 (n, b, q, f ) in (5) is obtained with (21) as (23) Note that the parameters in (22) and (23) are as stated in Table 5. Furthermore the presence of the biofilm growth environment facilitates diffusion and especially nutrient transport in comparison to the situation where only bacteria are present, which is included into the nutrient diffusion term as Although S. aureus bacteria do not have a flagella, they are able to spread on soft agar with a velocity of approximately 100 µm min . This spread is inhibited by the secretion of inhibitors against colony-spreading such as δ-toxin (Omae et al. 2012). Due to (14) and the fact that both concentrations are zero initially, we assume that [δ-toxin] = α q a. Note that δ-toxin belongs to the family of PSM and inhibits colony spreading, but not the growth rate of S. aureus (Omae et al. 2012). Therefore it is included as a factor 1 1+α q q into the diffusion term of the bacteria and we obtain Furthermore, we use Fickian diffusion for f as we consider the biofilm indicator molecule PIA. However, the involved modeling of the replicative bacteria diffusion also influences f as the presence of b is necessary for the production of PIA. We note that the choice of linear diffusion for the biofilm growth environment f allows this component to diffuse outside the colony, such that the bacteria can grow and spread in the growth environment.
Other PSM than δ-toxin, such as PSM-α and PSM-β, are stimulants of biofilm spreading (Omae et al. 2012), which are able to lyse eukaryotic cells, especially in competitive settings with several bacteria types (Cheung and Manna 2005). In PSM mutants, PSM-induced lysing is not available and thus the bacteria growth parameters G 1 , G 2 in (1)-(2) are smaller than in the wildtype bacterium. Note that, since PSM production depends on the availability of phosphorylated AgrA, also agr mutants do not produce PSM and the same changes apply.
Amyloid fibrils provide structural integrity to biofilms (Romero et al. 2010) and increase resistance to degradation (Schwartz et al. 2012). The production of amyloid fibrils depends on the agr locus, such that in agr mutants, the diffusion coefficient σ of the replicative bacteria increases. The substance α-toxin is necessary for cell-to-cell interaction during biofilm formation. Colonies of bacteria lacking α-toxin are unable to adhere to plastic surfaces under static or flow conditions (Caiazza and O'Toole 2003). Since the initial microbial adhesion to surfaces is a complex process involving bacterial factors as well as physical interactions like Lifshitz-van der Waals forces, electrostatic forces, acid-base interactions and Brownian motion forces (Cerca et al. 2005), we do not consider α-toxin effects. Protein A is most important in the context of invasion of a biological host (DeDent et al. 2007), where it facilitates colony spreading and virulence. In attachment or spa mutants the parameter δ in (1)-(2) is increased since a decrease in spa activity is linked to an increase in agr activity. Furthermore, the decrease in cell wall anchoring of protein A in the spa mutant decreases the stability of the biofilm, which results in an increased nutrient diffusion coefficient. The detailed mutant parameter changes are described in Table 7 in the following section.

Numerical simulation results and real data comparison
In this section we compare numerical simulations of our derived partial differential equations to biological observations of S. aureus pattern formation in the laboratory. The experimental observations are extracted from the published article by García-Betancur et al. (2017) and show S. aureus aggregates growing on top of an Mg 2+ -enriched TSB medium (TSBMg), where clinical isolates can develop robust multicellular aggregates.
For the simulations of our system, we use a time-adaptive finite element method implemented in MATLAB. In this finite element method we use first order linear finite elements with mass lumping on a triangular unstructured grid (Braess 2003;Brenner and Scott 2008). These are then coupled to an implicit Euler method for the time discretization. Additionally, for the time discretization, the size of the time steps used here is steered by the amount of Newton steps needed to solve the linearized equation system. The time step size in our method takes values between an upper and a lower bound, which both ensure convergence of the method, and the adaptive choice of the time step speeds up the computations. Also a constant time step could be used for the computations. Furthermore, it should be noted that for all the following simulations convergence studies have been performed to avoid an influence of the grid on the resulting patterns.
We further use a random diffusion component from a triangular distribution and initial conditions as introduced by Kawasaki et al. (1997) to reflect the stochastic fluctuation of the random movement of bacteria and not fully homogeneous material properties. This means that the diffusion parameter σ has the form σ = σ 0 (1 + Δ), where σ 0 = 0.5 and Δ is taken from a triangular distribution supported by [−1, 1]. This fluctuation is applied for every grid unit once in the initialization. Thus one parameter value is chosen randomly in order to break the symmetry of the solution. We note that simulation results without the random diffusion component show symmetric patterns as depicted in the paper by Horger et al. (2015).
As boundary conditions, we choose no-flux conditions as the colony grows on a nutrient medium dish which neither component can leave, i.e.
where Ω ⊂ R 2 denotes the area of the dish, which corresponds to a rectangular computation domain. As the initial concentration of nutrients is homogeneous throughout the medium, we define the initial nutrient concentration as n(x, 0) = n 0 for x ∈ Ω and n 0 = 1.11. Furthermore, the initial concentrations of non-replicative bacteria, quorum sensing and biofilm component are set to zero as these substances are not present initially. The replicative bacteria are inoculated as a drop in the center of the dish, which is reflected by the initial concentration We note that in the region outside the colony we find the TSB medium, which is a liquid enrichment medium. It shows only very little resistance against spatial expansion, such that these effects can be neglected in our model, especially since the colony is growing on top of the medium. Nutrients are present in the entire computation domain Ω and diffusion takes place inside and outside the colony.
Furthermore, the experimentally observed colonies are three-dimensional structures, whose height is small compared to their radius. We thus consider a twodimensional simulation domain, on which we add the densities of replicative (b) and non-replicative (s) bacteria as well as the concentration of PIA as a representative of biofilm environment density ( f ). In the experimental results, it is not possible to distinguish between the single components. Thus, for the comparison to the experimental observations, we consider b + s + f , as they together correspond to the visible and comparable part of the experimental microcolony. At this point it should also be noted that here the height of a structure can only be interpreted as a concentration level. Since the dimensionless concentration of PIA is in the range of 0-4.5, these values are indicated by a light blue color and the larger density values, which range up to 14 in areas where bacteria are present, are indicated by shades of yellow and red.
We observe that many colonies show a concentric ring structure. Outside this ring structure, fingering behavior takes place, which is mainly influenced by the replicative bacteria. The structure inside the ring is due to the accumulation of non-replicative bacteria and varies for different experiments, even for the same mutant. Thus the features of the experimental observations we focus on in the numerical simulations are the structure of the fingering pattern in the outer areas of the colonies, the presence of a ring, but not its explicit structure, and the prominence of the surrounding environment layer. In order to investigate the structure inside the ring further, the biological  Table 6 in (b) in comparison to experimental observation by García-Betancur et al. (2017) in (a) (colour figure online)    processes underlying the transitions between the replicative and the non-replicative state as well as the interactions of the non-replicative bacteria with the environment could be studied in detail from a biological point of view and included into the dif-ferential equations. The real data image for the wildtype bacterium in Fig. 8a shows a distinct structure with a concentric ring of bacteria and a narrow fingering or wrinkling structure towards the outer area of the colony. The colony is mostly round with only few very shallow dents. Furthermore, there is a layer surrounding the bacteria. For the wildtype bacterium we choose moderate parameter values as indicated in Table 6, since all regulation subsystems are active and contribute to the growth of the colony. The real data features are reproduced in the simulation result depicted in Fig. 8b, where a distinct ring structure is observed, the colony fingers are close and a layer of biofilm environment surrounds the colony.
The changes in the parameters for the mutant colonies are indicated in Table 7. In Fig. 9 the ica, spa and ica/spa mutant simulation results are depicted. In the real data in Fig. 9a, we observe that the extracellular matrix or ica mutant shows a less pronounced colony structure, especially in the inner area where the ring structure is weaker than in the wildtype colony. The colony is round, very narrow wrinkles are observed in the outer parts and the colony is not surrounded by a layer. The parameters G f = μ f = 0 are chosen since the ica locus is deactivated. Furthermore the lack of biofilm environment leads to a slower availability of nutrients for colony growth, such that the growth rates G 1 and G 2 are decreased, and without the added structure of the biofilm, the bacteria move faster. In the simulation results depicted in Fig. 9b we observe that a good agreement is reached as the ring structure is less pronounced and the branches are very close.
In the attachment mutant, depicted in Fig. 9c, the spa locus is disabled. A very pronounced ring structure is observed and, in comparison to the wildtype, the spaces between the wrinkles in the outer part of the colony are larger. The colony is surrounded by a layer. Since a decrease in the concentration of protein A is linked to an increase in the AIP concentration, we increase the parameter δ. Protein A binds to the cell wall envelope (DeDent et al. 2007), thus providing stability to the biofilm. Due to the decreased stability, we slightly increase the nutrient diffusion coefficient d n . In the simulation result in Fig. 9d we observe a good agreement of the simulations as the distances between the fingers are increased, while a very distinct ring structure is observed.
The combined spa and ica mutant in Fig. 9e shows features from both mutations. Here the concentric ring structure is more pronounced than in the ica mutant and the outer areas of the colony show wider gaps between the wrinkles in the colony structure. As above, there is no biofilm formation due to the ica locus mutation. In the simulation we combine the parameters G 1 , G 2 and σ from the ica mutant with the parameters δ and d n from the spa mutant to obtain the simulation result depicted in Fig. 9f. In comparison to Fig. 9b, it shows a stronger ring structure and wider gaps between the fingers, and thus a good agreement of the real data and the simulation. Figure 10 shows the simulation results for the psm-α and psm-β mutants. The amyloid type 1 or psm-α mutant real data in Fig. 10a shows a pronounced ring structure and very narrow wrinkling in the outer areas of the round colony. In comparison to the wildtype colony, we observe a more prominent layer of biofilm environment surrounding the colony. Due to the lack of PSM, the POPC vesicle lysing capacity decreases. Since this capacity is smaller for PSM-α than for PSM-β, also the decrease is smaller. Due to the structuring effect of PSM, we decrease the parameter d n . Since   psm-β mutant depicted in Fig. 10c, we also slightly increase the replicative bacteria diffusion parameter since the lack of biofilm structure facilitates bacteria movement. We again decrease d n and, in contrast to the psm-α mutant, we increase d f to obtain the simulation results depicted in Fig. 10d, where again a good agreement of the simulation result is observed as very little structure and very little wrinkling are observed in a round colony shape. The QS or agr mutant real data displayed in Fig. 11a shows a distinct structure in the middle of the colony, while towards the outer areas, the colony flattens and shows more pronounced dents. The colony is surrounded by a slightly uneven layer. Since the agr locus is not active, we take G q = μ q = 0. No agr activity also means that there are no PSM, and thus the growth parameters are decreased as in the psm mutant colony. Furthermore, since the biofilm is less structured without PSM, the diffusion coefficient of the replicative bacteria is increased and the diffusion parameters for nutrients and biofilm are chosen between the values for the psm-α and psm-β mutant colonies. This is due to the fact that biofilm-enhancing effects in psm-α and psm-β mutants do not seem to be additive, as in S. epidermidis (Le et al. 2014;Wang et al. (a) Agr mutant colony evolution.
(b) Agr mutant simulation result. ). We observe a good agreement of the simulation result depicted in Fig. 11b, where the structure in the middle of the colony as well as the slightly uneven layer are reproduced.

Conclusion
We investigated pattern formation in several mutant colonies of the bacterium S. aureus. Starting from the mathematical modeling of the agr, sarA, ica and sarA homologue cellular gene regulation systems, we derived a system of non-dimensional ordinary differential equations. Time-scale arguments and parameter size considerations allowed us to significantly reduce the system complexity, yielding differential equations and diffusion functions for the effects of quorum sensing and biofilm, where the parameters for the mutant colonies were derived from the gene regulations. We developed a new model for pattern formation in S. aureus and, in order to compare our results to experimental observations, we performed finite element simulations for several mutant colonies. Using the parameter changes indicated by the gene regulation mechanisms, we were able to adequately display the qualitative biological features of pattern formation in the simulations, thus establishing a basis for qualitative predictions of S. aureus pattern formation, which improves existing models by a reliable prediction of mutant colony growth. This work therefore constitutes a foundation for the further study of complex mutations in multiple gene loci.