Chronic Inflammation in the Epidermis: a Mathematical Model

The epidermal tissue is the outmost component of the skin that plays an important role as a first barrier system in preventing the invasion of various environmental agents, such as bacteria. Recent studies have identified the importance of microbial competition between harmful and beneficial bacteria and the diversity of the skin surface on our health. We develop mathematical models (M1 and M2 models) for the inflammation process using ordinary differential equations and delay differential equations. In this paper, we study microbial community dynamics via transcription factors, protease and extracellular cytokines. We investigate possible mechanisms to induce community composition shift and analyze the vigorous competition dynamics between harmful and beneficial bacteria through immune activities. We found that the activation of proteases from the transcription factor within a cell plays a significant role in the regulation of bacterial persistence in the M1 model. The competition model (M2) predicts that different cytokine clearance levels may lead to a harmful bacteria persisting system, a bad bacteria-free state and the coexistence of harmful and good bacterial populations in Type I dynamics, while a bi-stable system without coexistence is illustrated in the Type II dynamics. This illustrates a possible phenotypic switch among harmful and good bacterial populations in a microenvironment. We also found that large time delays in the activation of immune responses on the dynamics of those bacterial populations lead to the onset of oscillations in harmful bacteria and immune activities. The mathematical model suggests possible annihilation of time-delay-driven oscillations by therapeutic drugs.


Introduction
The skin is the largest tissue, which is composed of several different layers.The epidermis is located at the outmost part of the skin tissue, which acts as a first barrier for the invasion of physical (water), chemical (proteins) and biological (virus and bacteria) agents.A population of keratinocytes is the major cell type in the epidermis, which constitutes stratum basale, stratum spinosum, stratum granulosum and stratum corneum.Keratinocytes release anti-microbial peptides or pro-inflammatory cytokines to prevent bacterial or viral infection [1].The second outmost layer, the dermis, is situated between the epidermis and subcutaneous tissues, which are composed of fibroblasts, macrophages and adipocytes [2].The dermis contains extracellular matrix components, including collagen and elastin, as well as lymph, blood vessels and many skin-resident immune cell types.The homeostasis of the skin tissue is maintained by appropriate elimination of invading agents and tight regulation of cellular activities.On the other hand, the breakdown of the homeostasis of the skin tissue induces numerous diseases, including cancer, complications after an injury and inflammatory symptoms.Atopic dermatitis (AD) is one of the major skin inflammatory diseases, which are characterized by the elevated level of serum IgE and chronic allergic immune responses [3].Incidence of AD has been increasing in developed countries.Notably, recent genetic studies have revealed that barrier dysfunction of the epidermis due to filaggrin mutation is a major triggering factor of disease progression [4,5].Filaggrin is synthesized in keratinocytes at the stratum granulosum, which is degraded to become a major component of natural moisturizing factor (NMF).Lack of NMFs is associated with water loss skin dryness, leading to the progression of AD and ichthyosis vulgaris [4].Not only dry condition, but also excessive proteolytic activities of the epidermis are implicated as causal factors of AD.Patients with Netherton syndrome who exhibit atopic dermatitis-like chronic inflammation indicate a genetic defect causing excessive serine proteases [6].
The inference of bacteria as an environmental factor has been implicated as a possible factor for the progression of skin inflammatory diseases.Two major pathogenic bacteria species for skin diseases are Staphylococcus aureus (S. aureus) and Streptococcus pyogenes (S. pyogenes).Importantly, S. aureus is a major virulent species, which is implicated to be associated with the progression of atopic dermatitis [7].S. aureus also induces impetigo and another serious symptoms.Methicillin-resistant S. aureus (MRSA) is an antibiotic-resistant strain of S. aureus, the incidence of which is often reported as a nosocomial infection [8].On the other hand, some commensal bacteria can exhibit mutualistic behaviors through the suppression of potentially pathogenic bacterial species via direct and indirect interactions, known as probiotic effects.For instance, Staphylococcus epidermititis(S.epidermititis), a major commensal bacterial species in the skin, can support the host defense by releasing antimicrobial peptides [9,10].The other beneficial microbial species include species belonging to the Lactobacillus genus.Lactobacillus reuteri helps keratinocyte survival from S. aureus-induced cell death by outcompeting S. aureus [11].
The immune system plays a major role in preventing the invasion of numerous agents, including bacteria, fungi, virus and foreign proteins [12].Not only foreign antigens, but also antigens presented by commensal bacteria can be an antigenic stimulation for the host immune system [13].In fact, the number of bacteria in abundance is controlled by the host immune system under normal conditions [14].The pathogenicity of S. aureus can be conferred by numerous immune evasion strategies.In fact, several virulent factors of S. aureus have been reported in [15][16][17].On the contrary, several beneficial roles of S. epidermititis have been reported, although S. epidermititis can be virulent as a nosocomial pathogen for immunocompromised patients [18].S. epidermititis triggers innate immune responses via toll-like receptor (TLR)-2, which mediate the killing of pathogenic bacteria, such as S. aureus [19].Lactobacillus plantarum can utilize the host innate immune system mediated by epithelial cells by modulating the IL-17, IL-23 and TLR-2/4 expressions [20].
Regardless of the fact that many causal and preventive factors for the progression of AD and other skin inflammatory diseases have been identified, each experimental and clinical research only focuses on a specific aspect of the skin biology.The integration of knowledge in each sub-domain is needed in order to achieve a comprehensive understanding of the progression of AD.As described above, the detection of the manifestation of atopic dermatitis requires the integration of weakened barrier function due to a genetic defect or excessive proteolytic activity, the inflammatory response triggered by some of commensal bacteria and abnormal recruitment of immune cells via irregular cellular communication with respect to cytokine signaling.Mathematical modeling and simulation study enable integration and provide some basic insights into the maintenance mechanism underlying the skin homeostasis and disease development as a defect of the homeostatic condition.
We focus on competitive bacterial interactions among S. epidermatitis as good bacteria and S. aureus as bad bacteria that occur at the skin tissue to specifically reflect experimental and empirical observations, such as [19].Although our primary focus is to investigate the effects of bacterial competition on the dynamics of inflammatory responses in the epidermis, the mathematical models presented here can be useful as a general scheme to describe the interactions among bacterial species as an environmental factor with host immune responses on the surface of the body, such as the epidermis and gastro-intestinal (GI) tract.See [21,22] for the diverse roles of proteases in the GI tract in the maintenance of intestinal homeostasis.Hence, we focus on constructing mathematical models to represent a less detailed, but the general manner of interactions among inflammation-related molecules, such as protease, transcription factors and extracellular cytokines with bacterial species in epidermis.
In the present paper, we investigate how chronic inflammation can occur at the skin tissue.Simple mathematical models are employed to describe the invasion of bacteria, the proteolytic activity of keratinocytes, the activation of innate immune response and the release of antimicrobial peptide and cytokines.The organization of the present paper is as follows.In Section 2, two mathematical models (M1 and M2) are formulated.The M1 model is formulated to model the dynamics of inflammation in response to bacterial infection via a transcription factor and extracellular cytokines, as well as active proteases.Time delays may play a central role in the regulation of the bacterial-immune system in this study due to possible delays in the regulation of the intracellular transcription factors, protease induction and secretion of extracellular cytokines from bi-directional communication between a cell in the tissue and the microenvironment.Artificial manipulation (inhibition or enhancement of molecular players) of signaling pathways by therapeutic drugs typically induces time delays [23,24] in generating the final production of immune responses, i.e., extracellular cytokines, such as IL-4, IL-12, TNF-α and IFN-γ [25,26].To explicitly describe the competition between harmful and good bacteria, a formulation via an ordinary differential equation (ODE) and a delay differential equation (DDE) is employed in the M2 model.Mathematical analyses on these models, including the existence and stability of equilibria, are discussed in Section 3. In Section 3.1, numerical simulations are performed to investigate how prominently the chronic inflammatory state is established and maintained.Moreover, we investigate how bacterial competition can lead to a high chronic inflammation state as an imbalanced state (dysbiosis).In Section 3.2, we analyze the competition dynamics of harmful and good bacteria in the absence and presence of time delays and immune boosting drugs.We also perform several in silico experiments, which include the investigation of: (i) the effect of the clearance speed of cytokines on generating three regimes (harmful bacteria persistence, good bacteria persistence and co-existence) in Type I dynamics and a bi-stability system as a possible phenotypic switch in Type II dynamics; (ii) the effect of time delays on generating the oscillatory behaviors of bacterial populations; (iii) the impact of different drug injection regimes on the bacterial populations.In Section 4, we provide a discussion on the fundamental mechanism of the bacterial attacks and immune response, as well as the survival schemes of harmful bacteria in competition with good bacteria and future work in detail.Nondimensionalization and sensitivity analysis of the model are given in the Appendix.

Materials and Methods
In this paper, we present two kinds of mathematical models, the M1 basic model in Section 2.1 and the M2 competition model in Section 2.2.

M1 Model
In this section, we develop a mathematical model based on the schematic diagram in Figure 1.As indicated in Section 1, the key main players of the bacterial infection network in absence of competition with other bacteria are the following variables: B = density of harmful bacteria, P = concentration of protease, A I = concentration of the intracellular transcription factor, A E = concentration of extracellular cytokines, Figure 1A illustrates the dynamical regulation of immune activities in response to bacterial infection.Bacteria grow in the system with a carrying capacity and induce the secretion of proteases for enhanced bacterial invasion.These proteolytic activities are suppressed by protease inhibitors under normal conditions, but the perturbed balances between a protease and its inhibitors induce the activation of transcriptional factors, causing skin troubles, such as atopy.The upregulated transcriptional factors induce immune activities (extracellular cytokines), which in turn try to kill bacteria.In order to incorporate the biological interactions shown in Figure 1A into our model of bacteria-immune dynamics, we began by simplifying this network.Figure 1B shows a representation of Figure 1A.The kinetic interpretation of the arrows and hammerheads in the given network represents induction (arrow) and inhibition (hammerhead).We merged all complex networks between proteases and their inhibitors into one component (blue dotted box in Figure 1A), while we kept the components of bacteria, transcriptional factor and external allergic immune responses (cytokines) in one module (red dotted boxes in Figure 1A), respectively.The scheme includes bacterial growth, the secretion of proteases from bacterial invasion and stimulated transcriptional factors, the activation of the intracellular transcription factor from upregulated proteases and secreted cytokines, the activation of extracellular cytokines from the transcription factor, protein degradation of those key molecules and eradication of those bacteria by cytokines.Activation of transcription factors that are associated with immune responses, such as a member of the interferon regulatory factor (IRF) family, is often mediated by positive feedbacks among these transcription factors [27].It is also known that activation of protease is mediated by molecular interactions among the members of the kallikrein family [28].Production of the cytokines is also facilitated by a positive feedback, which is often referred to as a cytokine storm [29].These activations with positive feedbacks can be modeled by Hill functions.
In this work, we consider the following specific type of functional response known as the Hill function: where • ∈ {IP, BP, PI, BI, EI, IE}.Assume that the activation of protease is mediated by transcription factor A I and bacteria B. Based on these observations, we write the phenomenological equations for the rate change of those key players (B, P, A I , A E ) as follows: where λ is the source of bacterial populations in the tissue from the air, r B is the growth rate of bacteria, K is the carrying capacity of the bacterial population, γ is the killing rate of bacteria by immune cytokines and, finally, δ P , δ I , δ E are the decay/clearance rates of proteases, transcriptional factors and cytokines, respectively.A list of parameters is summarized in Table 1.

M2 Model
In this section, we consider two different bacterial strains.Figure 2A illustrates the dynamical regulation of bacterial infection and immune responses.There exists competition between harmful and good bacteria for bacterial growth.Bacterial infection induces upregulation of the transcriptional factor within the cell for immune activity and the secretion of proteases for enhanced bacterial invasion.Induced extracellular cytokines from the transcription factor then suppress both of those bacteria.In order to incorporate the interaction network shown in Figure 2B into our model of bacterial competition and inflammation, we began by simplifying this network.As indicated in Section 1, five main players of the bacterial infection network are harmful bacteria, good bacteria, protease, intracellular transcription factors and extracellular cytokines.Let the variables B 1 , B 2 , P, A I and A E be the densities or concentrations of harmful bacteria, good bacteria, protease, transcription factor and cytokines, respectively.Figure 2B shows a representation of Figure 2A.The kinetic interpretation of the arrows and hammerheads in the given network represents induction (arrow) and inhibition (hammerhead).The scheme includes bacterial growth of both harmful and good bacteria, mutual inhibition between harmful bacteria and good bacteria, secretion of proteases from those bacterial invasions, activation of the intracellular transcription factor from these bacterial infections, activation of extracellular cytokines via the transcription factor, protein degradation of those key molecules and eradication of those bacteria by cytokines.
It has been known that: (i) the half-life of proteases (P) is short, which indicates the large decay rate of P and that protease reactions occur quickly; and (ii) typical chemical reactions among proteins and genes at a fast time scale lead to the fast internal dynamics.This allows us to use quasi-steady state approximation (QSSA) to simplify the complex models (or interaction networks in Figure 2B).Based on the topological structure and uni-directional activation flows in the immune reaction module (transcription factor activities (A I ), proteolytic activities (P), cytotoxic levels (A E ); gray box in Figure 2B) and the corresponding QSSA, we merged all complex networks between transcriptional factors (A I ), proteases (P) and extracellular cytokines (A E ) into one component (gray box in Figure 2B,C), while we kept the harmful (B 1 ) and bad (B 2 ) bacterial components in one module (yellow dotted boxes in Figure 2B), respectively.Based on these observations, we write the phenomenological equations for the rate change of those key players (B 1 , B 2 , A E ) as follows: where r 1 and r 2 are the growth rate of harmful and good bacteria, respectively, α ii (i = 1, 2) and α ij (i = j, i, j = 1, 2) represent intra-specific and inter-specific competition coefficients between harmful and good bacteria, respectively, γ is the killing rate of those bacteria by the extracellular cytokines, β 1 and β 2 are the activation of immune responses (level of cytokines) from harmful and good bacteria, respectively, and δ is the clearance rate of the immune cytokines.In many cases, measured values of these parameter values (r 1 , r 2 , α ii , γ, β 1 , β 2 , δ) are not available due to technical reasons.In order to determine appropriate ranges of parameter values for correct dynamical behavior reflecting a real biological system and to investigate the sensitivity of the bacterial populations and immune responses to these parameters, we have performed sensitivity analysis for a mathematical model ( 3)-( 5) in Appendix B. Some of the parameter values (α 11 , r 2 , γ, β 1 ) are very sensitive, but others (α 22 ) are not sensitive to these changes.See the Appendix for more details.For computational purposes, we nondimensionalize the variables and parameters of the M1 (Equation ( 2)) and M2 models (Equations ( 3)-( 5)) in Appendix A.
Ecological system in the presence of compe33on between harmful and good bacteria

Results
In this section, we analyze the dynamics of two models (M1 and M2).In the next section, we first investigate the dynamics of M1 model for bacterial infection and its implications on immune responses.

Dynamics of the M1 Model
The M1 model deals with immune responses to bacterial infection via the regulation of proteases (P), internal allergic immune responses (A I ) and the external allergic immune response in terms of cytokines (A E ).We first classify the steady states of the M1 model (2) in the next section.

Classification of Steady States
There exist three types of equilibria for System (2).Let E B := ( B, 0, 0, 0) denote a steady state representing no protease and immune responses under bacterial persistence, where B > K is a positive root of B2 − K B − Kλ/r = 0. Let E * := (B * , P * , A * I , A * E ) denote an inflammatory state that is maintained by additional external stimuli from bacteria.The components of E * are determined by the solution of the following system of equations: By substituting the first and third equations into the second equation of ( 6), we obtain the following equation with respect to A * I : It follows from the first and fourth equations of ( 6) that B * is explicitly written as a positive root of the following quadratic equation with respect to B: Note that B * > 0. A * I must satisfy: Then, there exists a unique positive root of (8), denoted by B = B * > 0. Since σ IE (A I ) is continuous and monotonically increasing with respect to A I , ( 9) is rewritten as: Then, ( 7) is rewritten as: The existence of positive equilibrium E * is determined by the root of (11) with Constraint (10).
Let E * L , E * U and E * H denote three equilibria of (2) ordered by the value of A I s: • represents the strength of inflammation triggered by bacterial antigenic stimuli.
Figure 3 indicates that there are two possible cases: the existence of a unique equilibrium for weak activation of proteases (m IP = 8; Figure 3A) or multiple equilibria for enhanced activation of proteolytic activation (m IP = 50; Figure 3B).Here, m IP is the activation rate of proteases from the Appl.Sci.2016, 6, 252 9 of 34 transcription factor in the cell.In the upper panels of Figure 3, the straight solid line (blue) and the dotted curve (green) show the left-and right-hand sides of (11) as a function of (A * I ), respectively.The intersection of those two curves represents the equilibria (E * L , E * U and E * H ). Stability analysis indicates that: (i) when m IP is small (m IP = 8.0, the upper panel of Figure 3A), E * L (black filled circle) is stable; (ii) when m IP is relatively large (m IP = 50.0, the upper panel of Figure 3B), two steady states (E * L and E * U ; empty circles) are unstable, but one steady state (E * H ; black filled circle) is stable.Figure 3A shows the emergence of the bacterial persistence phenotype in response to the weak activation of proteases (m IP = 8.0).The unique positive stable equilibrium E * L resides in the region (pink box) where the transcription factor activities are suppressed and bacterial growth is active.On the other hand, when protease activation is enhanced more than six-fold (m IP = 50), there exist three positive equilibria (E * L (unstable), E * U (unstable) and E * H (stable)) simultaneously (Figure 3B).The stable steady state E * H resides in the region (blue box in Figure 3B) where bacterial activities are inhibited by persistent internal immune responses.These results predict the dynamical changes for various levels of protease activation and illustrate the importance of protease activation in the regulation of bacterial growth under the surveillance of the immune system in the tissue.In the next section, we analyze the dynamics of the M1 model and discuss the implication of the internal and external immune responses on the regulation of the harmful bacteria population.

Dynamics of the M1 Model System
When the M1 model system ( 2) is in equilibrium, we can solve bacteria density (B) as a function of the activation rate of proteases from the transcription factor (m IP ) for any set of parameters.Figure 4A shows the graph B = B(m IP ) as the S-shaped curves when other parameter values are fixed as in Table 1.While a portion of the upper branch in the lower protease activation range is stable, the remainder of the upper branch corresponding to the intermediate range of the protease activation rate is unstable.On the other hand, the lower branch is stable, and the middle branch is unstable.If m IP is small, then the system (2) is in the upper branch, B is high and the bacterial persistence phenotype emerges.This situation continues to hold as m IP is increased until it reaches criticality.At this point, the system jumps down to the low branch, with a suppressed level of bacterial activities, and the bacterial growth is inhibited (while immune activities are increased).As m IP is decreased from a high level of protease activation, the bacterial growth remains suppressed, until m IP is decreased to the left knee point (red arrow; ∼42), at which time, the bacterial population jumps to the upper branch, and the bacteria return to the growth phase.Figure 4B-D also shows the graphs A I = A I (m IP ), P = P(m IP ) and A E = A E (m IP ) as the hysteresis loops, as well.One notes that the bifurcation curves for those variables in immune responses (intracellular transcription factors (A I ), protease level (P) and extracellular cytokines (A E )) show the flipped images of the B − m IP hysteresis loop in Figure 4A, reflecting the bacteria-immune competition system.In other words, the immune activities (levels of A I , P and A E ) are suppressed compared to bacterial persistence in response to the weak protease activation, while high levels of immune responses are shown compared to inhibited bacterial growth in response to strong protease activation.
Based on the dynamics of the bacterial activities and cytokine levels observed above, we shall define two adaptive types of the bacterial infection system (bacterial persist (T B ) and immune boosting (T I ) systems) as follows: where th B and th AE are the threshold values of bacterial activities and cytokine level, respectively.With this definition (12), the unique stable equilibrium E * L in Figure 3A belongs to the region T B in the case of low m IP , while the stable steady state E * H in Figure 3B resides in the region T I in the case of high m IP .
In Figure 5, we show how the system adapts to the changes in the parameter m IP as predicted in the analysis above (Figure 4). Figure 5A-C illustrates two distinct patterns of the steady state (SS; red filled circles) of the dynamical system in response to a low (m IP = 30, Figure 5A), intermediate (m IP = 50, Figure 5B) and high (m IP = 80, Figure 5C) activation rate of protease (m IP ) in the B-A E phase diagram.A low level of protease activation from the transcription factor (m IP = 30) induces low cytokine levels and high bacterial infection (Figure 5A), while the intermediate or high activation level (m IP = 50, 80) leads to significant immune response and suppressed bacterial activities (Figure 5B,C, respectively) regardless of initial conditions.Figure 5D illustrates two distinct modes in the B-A E plane as described in ( 12): (i) the bacterial persist region (T B ) where bacterial growth is enhanced and cytokine levels are suppressed; (ii) the immune boosting zone (T I where the extracellular cytokine levels are increased and bacterial activities are inhibited.In Figure 5E-G, we show the time courses of bacteria density (B; red) and the concentrations of protease (P, pink), transcription factor (A I ; green) and cytokines (A E , blue) in response to three protease activation rates from the transcription factor (m IP = 30, 50, 80) corresponding to Figure 5A-C, respectively, with the initial condition In the next section, we investigate the dynamics of the competition M2 model.

B low high
Cytokine level ( ) All other parameters are fixed as in Table 1.(E-G) Time courses of the main variables (B, P, A I , A E ) for various activation rates (m IP = 30, 50, 80) with the initial condition: B(0) = 1.7,P(0) = 0.1, A I (0) = 0.1, A E (0) = 0.1.All other parameters are fixed as in Table 1.

Dynamics of Two Bacterial Strains Model (M2 Model)
We first investigate the existence of the equilibria of the M2 model ( 3)-(5).

Existence of Equilibria
By replacing the variable A E with I for notational purposes, the M2 model system is given by: Appl.Sci.2016, 6, 252 13 of 34 where There are four types of equilibria of System ( 13).E 0 = (0, 0, 0) is a trivial equilibrium representing that neither bacteria nor the immune response exist.Let E 1 = ( B1 , 0, Ī1 ) and E 2 = (0, B2 , Ī2 ) denote dominant equilibria in which either B 1 or B 2 exists.The explicit values of each component of E 1 and E 2 are given as follows, respectively.
Let E + := (B * 1 , B 2 2 , I * ) denote a positive equilibrium representing the coexistence of two bacterial species under the pressure of the immune response.It follows from the third equation of ( 13) that I * is given by: By substituting the explicit value of I * into the first and second equations of ( 13), we obtain the following linear system of equations with respect to B * 1 and B * 2 : Hence, the explicit values of each component of E + are given by: where D 0 is given by: For equilibrium E + to be a positive equilibrium requires B * 1 > 0 and B * 2 > 0. Define matrix A = {a ij } and vector b = (b 1 , b 2 ) T (i, j = 1, 2), such that: Then, Note that ( 20) and ( 21) are equivalent to: and: respectively.Define D 1 , D 2 , w 1 and w 2 by: Finally, we consider the conditions for the existence of E + .It follows from ( 22) and ( 23) that + if and only if: or: Note that the stability conditions of E 1 and E 2 are given by ( 34) and (37), respectively (see the next subsection for details).Conditions ( 34) and ( 37) are mutually exclusive with (26), but are identical to (27).In other words, coexistent equilibrium E + exists only if both E 1 and E 2 are unstable or locally stable.In the later case, the system would be expected to exhibit bistability between E 1 and E 2 .
In summary, the existence conditions of internal equilibrium E + are classified in Tables 2 and 3 according to the sign of r 1 − r 2 and w 1 − w 2 .

Table 2. Existence condition of E
In the next section, we check the stability of the equilibria of the M2 model (13).

Stability of Equilibria
Mathematical conditions for local stability of equilibria are derived based on the linearized equations around any of the equilibria ).The Jacobian matrix for E • is given by: The Jacobi matrix for E 0 is given by: Since r 1 > 0 and r 2 > 0, E 0 is always unstable.The Jacobi matrix for E 1 is given by: Characteristic equation P 1 (λ) = 0 defined for J(E 1 ) is given by: Let A 1 (λ) be defined by: Note that the trace and determinant of A 1 satisfy tr(A 1 ) < 0 and det(A 1 ) > 0. Hence, E 1 is locally asymptotically stable if: Note that (33) is rewritten as: In other words, ( 34) is equivalent to one of the following two conditions: Stability conditions of E 2 are derived from the Jacobi matrix for E 2 : In a similar way to E 1 , E 2 is locally asymptotically stable if: Note that (36) is rewritten as: In other words, ( 37) is equivalent to one of the following two conditions: We note that the system of differential equations for bacterial strains without any anti-microbial killing and recruitment and with the same growth rates of bacteria (r 1 = r 2 = r B ): reduces to the classical two-dimensional Lotka-Volterra competition model.If we assume that the magnitude of inter-specific competition is stronger than intra-specific competition, i.e., then a unique positive equilibrium of ( 38) exists, and it is unstable.It can be shown that solutions converge to either (1/α 11 , 0) or (0, 1/α 22 ) depending on the choice of the initial state.
In the next section, we investigate the dynamics of the competition model (M2) and immune system.

Dynamics of the Competition Model M2 in Response to the Immune System
We shall define three adaptive types of the competition system (harmful bacteria-persist (T B ), harmful bacteria-free (T F ) and co-existence (T c )) corresponding to regions in the B 1 -B 2 plane, including equilibria points (E 1 , E 2 , E + ) discussed in the previous section: The basic parameter set for the M2 model is given in Table 4.

Parameter Description Type I Type II
Inter-and intra-competition Intra-specific competition coefficient In Figure 6, we investigate the dynamics of the two-species model (3)-( 5) in the presence of the immune response for a base parameter set (Type I).r 1 = 1.5, r 2 = 1.0,K 1 = 0.5, K 2 = 2.0, α 12 = 1.0, α 21 = 1.0, β 1 = 0.1, β 2 = 1.0, γ = 1.0.Analysis in Section 3.2.2indicates that: (i) E 2 is stable if δ > 0.66666, whereas E 1 is stable if δ < 0.03333; and (ii) E + is expected to be stable if 0.03333 < δ < 0.66666.Figure 6A-C shows the trajectories (B 1 (t), B 2 (t)) of harmful and good bacteria populations for various decay rates of cytokines (δ = 0.033 (Figure 6A), 0.35 (Figure 6B) and 0.7 (Figure 6C)) with four initial conditions: B 1 (0) = 0.07, B 2 (0) = 0.45 (yellow curve); B 1 (0) = 0.38, B 2 (0) = 0.58 (green curve); B 1 (0) = 0.05, B 2 (0) = 0.1 (blue curve); B 1 (0) = 0.35, B 2 (0) = 0.02 (purple curve).The initial condition of the immune response was set to be zero (A E (0) = 0) for all cases.Figure 6E-F shows the trajectories of (B 1 (t), A E (t)) corresponding to Figure 6A-C, respectively.When the decay rate of cytokines is small (δ = 0.033), the system converges to E 1 equilibrium where good bacteria are cleared out and harmful bacteria survive in a battle via the immune system (Figure 6A).Initial strong immune responses (black arrow in Figure 6E) due to the weak clearance rate (δ 1) significantly eliminate both harmful and good bacterial populations (black arrowhead in Figure 6A).While the good bacteria are totally eradicated by this strong immune response due to the relatively low growth rate, the harmful bacteria survive due to the higher growth rate and winning the competition battle with the good ones (blue arrow in Figure 6A).The system adapts a harmful bacteria-free equilibrium when δ is large (δ = 0.7 in Figure 6C).The strong clearance of immune activities in the system leads to relatively weak immune responses (black arrow in Figure 6G).This increases the chances of winning the competition for good bacteria and decreases the harmful bacteria population, pushing the system (B 1 (t), B 2 (t)) in the upper left corner (black arrowhead in Figure 6C).Then, the system converges to the E 2 equilibrium, the attractor (red filled circle in Figure 6C), where harmful bacteria are eradicated by the immune system and the helpful bacteria persist.On the other hand, an intermediate immune response (δ = 0.3 in Figure 6B) leads to the co-existence of good and harmful bacteria populations.The immune system initially successfully attacks and decreases the number of both harmful and good bacteria, but this also reduces the activation of the immune system (cf.Equation ( 5)).The reduced immune activity also increases the chance of the regrowth of both bacterial types (arrowhead in Figure 6B), leading to the coexistence of those two bacterial populations (red filled circle in Figure 6B), corresponding to E + equilibrium.In response to small (δ = 0.0033), intermediate (δ = 0.35) and large (δ = 0.7) decay rates of cytokines, the system transits from harmful bacteria-persist region (T B ) to the co-existence zone (T c ) and then to the harmful bacteria-free (T F ) region (yellow curved arrow).Dysbiosis, or bacterial imbalance, represents a state of reduced species diversity with the emergence of a few extraordinary highly abundant species.Dysbiosis is broadly observed for several microbial ecologies, including in aquatic systems or intestinal systems.It has been shown to be associated with illnesses, such as cancer [30,31], bacterial vaginosis [32], inflammatory bowel disease [33,34], chronic fatigue syndrome [35], obesity [36,37] and colitis [38].Dysbiosis in the gut is known to be associated with major chronic inflammation states [39].Dysbiosis caused by the imbalance of the skin commensal bacterial species composition has been reported [40,41].Importantly, dysbiosis characterized by the increase of S. aureus has been shown to be associated with atopic dermatitis, one of the major skin inflammatory diseases [7].It is suggested that dysbiosis and Staphylococcus aureus colonization drive inflammation in atopic dermatitis [7].In our model, different immune reactions from the weak, intermediate and strong clearance strength (δ) of extracellular cytokines, such as IL-4, IL-12 and IFN-γ, result in the imbalance between harmful and good bacteria, co-existence or healthy tissue homeostasis (Figure 6D).In particular, increased levels of harmful bacteria and reduced levels of the beneficial bacteria, i.e., dysbiosis, may be induced when the clearance of the immune reactions is weak (δ = 0.033; Figure 6A).The tipping point of the balance between beneficial and harmful bacteria is the vigorous and subtle competition between those different kinds of bacteria.
Figure 7A shows the steady state of harmful bacteria (B * 1 ) at three equilibria (E 1 red; E 2 blue; E + black) as a function of δ.Solid and dotted curves illustrate the stable and unstable branches at E 1 , E 2 , E + for the continuous spectrum of δ.Small, intermediate and large values of δ lead to harmful bacteria persisting, co-existence and good bacteria persisting regimes, respectively, as shown in Figure 6A-C.Figure 7B illustrates how the system transits from T B to T c and then to T F in response to an increase in δ by following the stable branches of B * 1 at E 1 , E + and E 2 , respectively, in Figure 7A.This illustrates how phase transitions (T B → T c → T F ) in Figure 6D can be induced by the continuous increase in the clearance rate of the immune system.In Figure 8, we investigate the dynamics of the bi-stable competition system with the parameter set (Type II): In comparison to the previous case in Figure 6, the system does not present the co-existence region (T c ). Figure 8A shows the regions of harmful bacteria-persist (T B ) and bacteria-free (T F ) in the B 1 − B 2 phase-plane.While initial states (B 1 , B 2 ) in the region R B (red) converge to the harmful bacteria-persist equilibrium (T B , red circle), the initial states in the upper-left region (R F ; blue) converge to the harmful bacteria-free equilibrium (T F ; red circle).For example, curves indicate the trajectories (B 1 (t), B 2 (t)) for two very close initial conditions: B 1 (0) = 0.15, B 2 (0) = 0.713 (blue curve) and B 1 (0) = 0.15, B 2 (0) = 0.71 (red curve) near the boundary (green dotted curve) between R B and R F .The initial condition of cytokines was set to be zero (A E (0) = 0).E 1 = (0.9375, 0), E 2 = (0, 0.3158).Figure 8B,C shows the trajectories in the B 1 − A E plane (Figure 8B) and time courses (Figure 8C) of bacterial populations (B 1 , B 2 ) and cytokine level (A E ) for two very close initial conditions in Figure 8A: B 1 (0) = 0.15, B 2 (0) = 0.713 (dotted curves) and B 1 (0) = 0.15, B 2 (0) = 0.71 (solid curves).Figure 8D shows the bi-stable nature of the dynamical system where one kind of the harmful or beneficial bacteria dies out and the other kind survives depending on the initial state B 1 (0), B 2 (0).Therefore, the initial status of exposure to both harmful and beneficial bacteria determines the dysbiosis or healthy tissue.
In the next section, we investigate the effect of time delays in immune responses on the dynamics of the system.
As δ is increased, the system sequentially follows the stable branches in (A), leading to the transition from stable E 1 branch (red solid curve in (A)) to stable E + branch (black solid curve in (A)) and then to stable E 2 branch (blue solid curve in (A)).All other parameters are fixed as in Table 4 (Type I).  4 (Type II).

Effect of Time Delays in Immune Response on the Competition System
In this section, we introduce time delays in the immune response for the reduced competition system (3)-( 5).The governing equations for the simple model with time delays (τ 1 , τ 2 ) are given by: where τ 1 , τ 2 are time delays in the immune response for harmful and beneficial bacteria attacks, respectively.
In Figure 9, we investigate the effect of small time delays (τ 1 = τ 2 = 1.5).In the absence of time delays (τ 1 = τ 2 = 0), the population of bad bacteria (B 1 ) converges to zero (B s,2 1 = 0; blue solid line in Figure 9A), while the population of good bacteria (B 2 ) and the immune system (A E ) converge to positive equilibria B s,2 2 (blue solid line in Figure 9B) and A s,2 E (blue solid line in Figure 9C), respectively.On the other hand, in the presence of time delays (τ 1 = τ 2 = 1.5), the population of good bacteria (B 2 ) converges to zero (B s,1 2 = 0; red dashed line in Figure 9C), while the population of bad bacteria (B 1 ) and the immune system (A E ) persists with positive equilibria B s,1 1 (red dashed line in Figure 9B) and A s,1 E (red dashed line in Figure 9C), respectively.Therefore, an introduction of weak time delays in the system induces a switch from the B2-dominant equilibrium (0, B s,2 2 , A s,2 E ) to the B1-dominant equilibrium (B s, 1  1 , 0, A s,1 E ).See Figure 9D,E A E in (B,C) in the phase plane.The introduction of weak time delays in the system induces a switch from the B2-dominant equilibrium (0, B s,2 2 , A s,2 E ) to the B1-dominant equilibrium (B s,1 1 , 0, A s,1 E ).Initial condition: B 1 (0) = 0.15, B 2 (0) = 0.72, A E (0) = 0.0.All other parameters are fixed as in Table 4 (Type II).
However, for larger time delays, the DDE system completely changes the dynamics.The system induces oscillations in both the population of bad bacteria (B 1 ; Figure 10A) and the levels of immune cytokines (A E ; Figure 10C) in the presence of larger time delays (τ 1 = τ 2 = 3.5).The system maintains the extinction of the good bacterial population under this condition (Figure 10B) as in the small time delays in Figure 9.The dynamics of the ODE case is shown in the blue solid curves in Figure 10 and is the same as in Figure 9.This oscillatory behavior of the in vivo pathogens and specific/non-specific immunity was observed in experiments [42,43], and the time delay may have existed in the specific bacterial kinds.41)- (43).The introduction of time delays in the system induces oscillatory behaviors of both bad bacteria (B 1 ) and immune cytokines (A E ) and the extinction of good bacteria.(A-C) Time courses of the populations of bad bacteria (B 1 in (A)) and good bacteria (B 2 in (B)) and the immune response (A E in (C)) in the absence (blue solid lines) and presence (red dashed lines) of time delays; (D) the corresponding trajectories of B 1 and B 2 in (A,B) in the phase plane; (E) the corresponding trajectories of B 1 and A E in (B,C) in the phase plane; Initial condition: B 1 (0) = 0.15, B 2 (0) = 0.72, A E (0) = 0.0.All other parameters are fixed as in Table 4 (Type II).
Our investigation illustrates that the system undergoes dynamical changes as the time delays (τ = τ 1 + τ 2 ) are increased.In the absence (τ = 0) or small values of time delays, the bi-stable system induces either the imbalance state (B s, 1  1 , 0, A s,1 E ) or disease-free state (0, B s,2 2 , A s,2 E ).As the time delays are further increased, the system induces the oscillatory behaviors of harmful bacteria and immune activities for some initial conditions (B 1 , B 2 , 0).This indicates that the strength of time delays in either the induction of extracellular cytokines or the manipulation of intracellular signaling pathways may be enough to perturb the bistable pathogen-immune dynamics and leads to the recurrence of the harmful bacteria population.
In the next section, we investigate the local stability of the reduced system of DDEs ( 41)-( 43) and present the onset of Hopf bifurcation.

Local Stability Analysis of Delay Differential Equations
For simplicity, we consider the simplified version of the model for analysis by ignoring the good bacterial dynamics: Assume that δ 1. Quasi steady state approximation is applied to (44) to obtain two-dimensional system: where β 11 is explicitly given by: The linearized system of (45), which is defined for E 1 = ( B1 , 0), is given by: Characteristic equation P(z) = 0 defined for ( 47) is explicitly given by: Let us investigate whether P(z) = 0 has a pair of pure imaginary roots z = ±iω, where without loss of generality, we can assume that ω > 0.Then, the real and imaginary parts of P(+iω) = 0 are given by: r B1 α 11 + r B1 β 11 cos ωτ = 0, By adding the square of real and imaginary parts, we obtain that: The equality in (50) holds if and only if β 11 > α 11 .It follows from the first equation of (49) that: By numerical computations, critical time delay τ * at which a system undergoes Hopf bifurcation is determined for the parameter set in Figure 11.More precisely, τ * 1.418.Figure 11 illustrates the dynamic changes of stability at a steady state in Equation ( 45) as the time delay (τ) passes through the critical Hopf bifurcation point τ * 1.418.The stable state of the equilibrium for smaller τ's (τ = 0.0 (Figure 11A), 1.0 (Figure 11B) and 1.417 (Figure 11C)) becomes the unstable state around τ * (τ = 1.4181 in Figure 11C).

Therapeutic Approaches
Results in the previous section indicate that harmful bacteria may not be completely removed by typical immune responses due to recurrence in the presence of time delays, and one has to introduce a drug that enhances the immune system for the eradication of harmful bacteria.We introduce a drug (D) under the following assumptions: (i) drugs enhance the immune activities of extracellular cytokines by inhibiting signaling networks involving the intracellular transcription factors and protease activation; (ii) the drug is administrated with either a constant rate or periodic injection of drug; (iii) drug compounds have a natural decay.The governing equations then are given by: where τ 1 , τ 2 are time delays in the immune response for harmful and beneficial bacteria, respectively, as in the previous section, λ D is the injection rate of drugs and δ D is the decay rate of the drug.In Figure 12, we investigate the effect of constant drug injection on the dynamics of the DDE system ( 52)-( 55) in the presence of large time delays (τ 1 = τ 2 = 3.5).Figure 12A-C shows populations of harmful (B 1 , red curve) and beneficial (B 2 , blue curve) bacteria in the upper panels and the levels of cytokines (A E , solid black curve) and drugs (D, dotted green curve) in the lower panels for low (λ D = 0.01), intermediate (λ D = 0.03) and large (λ D = 0.1) injection rates, respectively.A small injection (λ D = 0.01) of drugs is not so effective to remove the harmful bacteria and still maintain the oscillatory patterns of both harmful bacteria (B 1 ) and immune cytokines (A E ).This still leads to the extinction of beneficial bacteria (B 2 (t) → 0 as t → ∞).See Figure 12A.For an intermediate level of injection (λ D = 0.03; Figure 12B), drugs annihilate the oscillations of the harmful bacteria population and immune responses, leading to the infection-persist state (B s 2 = 0, B s 1 > 0).On the other hand, a large amount of drug injection (λ D = 0.1; Figure 12C) significantly enhances the immune activity and eliminates both harmful and beneficial pathogens (B s 1 = 0, B s 2 = 0).Figure 12D-E shows the trajectories in the B 1 -B 2 (Figure 12D) and B 1 -A E (Figures 12E) planes, respectively, for the corresponding λ D 's (λ D = 0.01 (red), 0.03 (blue) and 0.1 (black)).) and good (B 2 ) bacteria for various drug injection rates (λ D = 0.01, 0.03, 0.1); (E) trajectories of bad bacteria (B 1 ) and immune activity (A E ) for various drug injection rates (λ D = 0.01, 0.03, 0.1).Initial condition: B 1 (0) = 0.2, B 2 (0) = 0.95, A E (0) = 0.0, D(0) = 0.All other parameters are fixed as in Table 4 (Type II).

Oscilla(on
In real intravenous injection, the drug is administrated in a periodic infusion; we investigate the effect of drugs in a more realistic situation in the clinical setting.For this, we replace Equation (55) with the following: where λ D is the injection rate, δ D is the decay rate of drugs, N D is the number of drug injections and is the characteristic equations that give one over the time interval [t j , t i + t d ] with duration t d and zero otherwise.Here, the injection period is fixed: τ D = t j+1 − t j , ∀j = 1, . . ., N D .
In Figure 13, we investigate the dynamics of the system in response to the periodic injection of the high and low doses of drugs that boost patients' immunity.For a low dose of drugs (λ D = 0.2, t d = 1 h, N D = 10), the system still maintains the oscillatory behaviors of bad bacteria, and the immune-boosting effect from drugs is not significant enough to eradicate the bad bacteria (red curve (B 1 ) in Figure 13A).See the relatively low immune responses during IV injection periods in Figure 13D (blue curve (A E ) for t < 240).However, the promoted immune activity in response to a higher dose of drugs (λ D = 1.0, t d = 1 h, N D = 10) results in bacterial extinction (red curve (B 1 ) in Figure 13C) due to elevated levels of initial immune response (blue curve (A E ) in Figure 13D).(Type II).
In Figure 14, we investigate the effect of therapeutic drugs on the regulation of the eradication or recurrence of harmful bacteria in response to various combinations of infusion rates (λ D ) and injection periods (τ D ).For a fixed value of injection period (τ D ), the system switches from the recurrence phase to the eradication phase for harmful bacteria as the injection rate (λ D ) is increased.For a fixed λ D , the larger interval length between drug injections tends to increase the chance of the recurrence of harmful bacteria.The model predicts that the larger infusion rate (λ D ) and shorter injection interval would lead to the elimination of harmful bacteria.(Type II).

Conclusions and Discussion
This paper investigates the progression of skin inflammatory disease by mathematical modeling and simulation.Three mathematical models were built to investigate how bacterial antigenic stimuli initiate and maintain the inflammatory response of keratinocytes.First, the effect of positive feedback regulation among protease and the transcription factor was incorporated, where we consider a single bacterial species as the source of antigenic stimuli (model M1).The existence of multiple positive equilibria indicated by equilibrium analysis implies that feedback switch occurs for model M1.To investigate how high inflammatory response is maintained, parameter m IP representing the activation rate of proteases from the transcription factor was varied.Model M1 exhibits qualitatively different types of behaviors: one is the persistence of bacteria under a low inflammatory state; another is a high inflammatory state.Numerical computations indicate that the transition from low to high inflammatory states can occur when parameter m IP varies (Figures 3-5).From the biological point of view, these computational results suggest that excessive protease activity can lead to a high inflammatory response.A general scheme presented in this paper applies to at least two different types of chronic inflammatory diseases.The first one is atopic dermatitis, which is often caused by primarily defection of the barrier function at the epidermis via excessive protease activity.In the previous study, the switch-like behavior was extensively investigated, which qualitatively explains the progression of atopic dermatitis [44].The second one is an inflammatory bowel disease, which is recognized as a major inflammatory disease in the gut.Several research works imply the association of excessive protease activity with the progression of chronic inflammation [21].
To investigate the effect of species competition among bacterial species, model M1 was extended to include commensal or beneficial bacteria, which compete with a harmful bacterial species.Mathematical model M2 was constructed to investigate the bacterial competition under immune suppression.Model M2 consists of a classical Lotka-Volterra competition model with immune suppression as a negative feedback effect on both species.Numerical simulations were implemented extensively to understand the qualitative behavior of two bacteria under immune suppression (Figures 6 and 8).The outcomes of competition among two species were classified by means of the existence and stability conditions of equilibria (Section 3.2).The mathematical condition for the stability of the dominant equilibrium in which only a single bacteria species exists was derived.Condition (E 1 S 1 ) represents a situation in which harmful bacteria outcompete beneficial bacteria when the growth rate of harmful bacteria is faster than that of beneficial bacteria (r 2 > r 1 ).Moreover, the second condition D 1 = α 21 K 1 − 1 > 0 represents the transversal eigenvalue of equilibrium E 1 in the B 2 -direction.D 1 > 0 implies that the direction of the transversal eigenvalue is negative.Hence, beneficial bacteria cannot invade and grow when the harmful bacteria exist and have reached their carrying capacity.Hence, the condition (E 1 S 1 ) implies the non-invasibility of beneficial bacteria.On the other hand, the condition (E 1 S 2 ) represents an interesting situation.Despite the fact that harmful bacteria have a slower growth rate than beneficial bacteria, this can prevent the invasion of beneficial bacteria (D 1 > 0) and, importantly, suppress the growth rates of both bacteria by utilizing the boosting of immune suppression.In other words, harmful bacteria take advantage of growth inhibition by immune suppression, which leads to the dominance of harmful bacteria.From the biological point of view, there exists a possibility that some of the bacteria favor the inflammatory condition that suppresses potential competitors.In [7], it is shown that S. aureus increases in abundance during the process of dysbiosis, which can drive the inflammatory response, leading to atopic dermatitis progression.Hence, the implication derived from the mathematical computation and analysis results can partly explain how S. aureus grows in abundance while suppressing potential competitors.
Finally, time delays were introduced to represent the time required to activate an immune response.In the skin tissue, this is generally mediated by innate immune cell types, such as neutrophils, which have to be recruited from peripheral blood to an infection site.Hence, a time delay is inevitable to consider the process of immune activation.By the introduction of a time delay, interestingly, beneficial bacteria can outcompete harmful bacteria even though having a disadvantage in competition (Figures 9-11).These numerical computation results indicate the possibility that the time delay in immune cell recruitment can affect the outcome when two species are competitively bistable.A case study on Francisella tularensis infection, which may causes hypercytokinemia, reported at least a one-day delay in neutrophil recruitment post infection [45].The introduction of a time delay in the recruitment of activated innate immune cells to the infection site in the skin exhibits interesting behaviors.Bacterial antigenic stimuli trigger the immune response via TLR4, antigen recognition receptors, which specifically detect lipopolysaccharide (LPS) expressed on the surface of Gram-negative bacteria.If we assume that good and bad bacteria are both Gram-negative bacteria, then it would be reasonable to assume that the time required to activate immune responses in response to bacterial stimuli would be the same.Hence, the same value was utilized for the two delays in this work.We also investigated the effects of the immune-boosting drugs on the selection of harmful or beneficial bacteria and developed strategies to prevent the cycles or recurrence of harmful bacteria populations (Figures 12-14).
The results of this work can serve as a starting point for better comprehensive modeling and experimentation.Some problems and extensions of the model that can be addressed in the future are as follows.

•
In the present model, we concentrated on two major pathogenic and commensal bacterial species to obtain basic insight into how microbial interactions mediated chronic inflammation.However, more than hundreds of bacterial species have been demonstrated to coexist in the skin tissue.Metagenomic analysis targeting the gut-and skin-resident microbiome has revealed that numerous uncultured species exist and potentially affect the maintenance of skin homeostasis, as well as the progression of skin inflammatory disease [46].The existence of spatial compartmentalization by forming heterogeneous clusters of colonies across the epidermis and dermis has been shown [47,48].Although a few numbers of dominant species exist in terms of population abundance, bacterial diversity in the skin is highly maintained [49].Complex interactions among commensal bacteria, the host immune system and different sources of environmental fluctuations should be essential factors for the maintenance of species diversity.Therefore, the incorporation of more than two bacterial species into the model would be more realistic.Colonization of harmful bacteria would be prevented by community-level resistance by a bacterial community.The incorporation of multiple species interactions will provide new insights on how the loss of bacterial diversity would lead to high inflammatory states.

•
We considered the same time delays in the M2 model in this work.There exists the possibility of an immune escape mechanism, which might justify the use of different time delays.For instance, certain types of bacteria downregulate antigenicity when they invade tissue in order to escape from immune surveillance [50].This would lead to a time delay in the activation of the immune system.Major extensions of the current model to include different time delays are warranted.

•
The mathematical models presented here do not distinguish immune cell types, which are crucial to determine the difference between the epidermis and the GI tract.For instance, Langerhans cells are the major resident immune cell type that stays below the second layer of the stratum granulosum (below the tight junction) and captures the antigen.After capturing the antigen, Langerhans cells move to a draining lymph node to present the antigen to lymphocytes, known as homing.In the intestine, invading bacteria that attach to the gut epithelial cells trigger inflammatory responses, and finally, these bacteria are eliminated by immune cells recruited from the Payer's patch or gastric mucosal lymphoid follicles.

•
Explicit incorporation of spatial structure is essential to represent specific and unique information to the epidermis or the GI tract.In the present paper, however, we focused on the role of bacterial species to induce inflammatory responses rather than spatial structure, which forms specific and unique interactions among invading bacteria and immune cells.The ongoing project aims to incorporate spatial structure and heterogeneity in immune cell subtypes, but it is currently under investigation.

•
The major signaling networks that control the intracellular regulation of transcriptional factor, proteases and protease inhibitors need to be addressed.

•
The microenvironment also plays an important role in the regulation of epidermis and stem cell dynamics [51].These include other immune cells, endothelial cells and stromal cells, such as fibroblasts, as well as growth factors secreted by these cells.

•
Cell-mechanical regulations, such as actin and serum response factor, were also shown to transduce bio-physical cues from the microenvironment to control epidermal stem cell fate [52].
Our understanding of the complex biochemical interactions between the epidermal cells and the microenvironment is very limited.Hybrid approaches may be used to take into account these intracellular signaling pathways in addition to the mechanical interactions of cells with microenvironment for detailed proteolytic activities, growth and invasion at the cellular level and viscoelastic response of the whole tissue [23,[53][54][55][56][57].Yet a more comprehensive understanding of the role of the microenvironment in epidermal homeostasis may lead to the development of new therapeutic agents.We will discuss these in future work.immune response is positively correlated with K 1 , β 1 , β 2 , but is negatively correlated with the decay rate (δ).In particular, the immune response is strongly negatively correlated with the killing rate (γ) of both bad and good bacteria at later times (t = 40, 80) due to decreased bacterial populations in response to increased immune response.The parameters α 12 , β 2 , K 2 (= 1/α 22 ) have little correlation with all variables (B 1 , B 2 , A E ).Table B1 summarizes the results of the sensitivity analysis in terms of the populations of bad (B 1 ) and good (B 2 ) bacteria and cytokine concentration (A E ), at t = 1, 40, 80.   3)-( 5)) for ten model parameters (r 1 , K 1 , α 12 , r 2 , α 21 , K 2 , γ, β 1 , β 2 , δ).Red color indicates positive correlations, and blue color indicates negative correlations between the main variable and each parameter at the given time.The analysis was carried out using the method of Marino et al. ( 2008) [66] with a sample size of 10,000.

Figure 1 .
Figure 1.A schematic of the M1 model.(A) A schematic of immune responses to bacterial infection; (B) the final network model that abstracts the key structure of the interaction network in (A).By merging a multi-species compartment (blue dashed box in (A) including proteases and their inhibitors) into a compartment ('P' in (B)), we get a simpler model in (B).Densities of bacteria, proteases, intracellular transcription factor and extracellular cytokines are represented by 'B', 'P', 'A I ' and 'A E ', respectively.

Figure 2 .
Figure 2. A schematic of the M2 model.(A) A schematic of the biological system for the competition between harmful and good bacteria and immune responses.There exists mutual antagonism between harmful and good bacteria.On the other hand, bacterial infection induces upregulation of the transcriptional factor (blue diamond) within the cell for immune activity and secretion of proteases (red quarter pie) for enhanced bacterial invasion.Induced extracellular cytokines (green) from transcription factor (blue) then suppress both harmful (red star) and good (blue star) bacteria.All arrows refer to the induction of gene expression or proteins.The hammerheads from and to bacteria (B 1 , B 2 ) refer to the inhibition or suppression of bacterial growth.(B) Topological networks representing the biological observations in (A).Densities of harmful bacteria, good bacteria, proteases, intracellular transcription factor and extracellular cytokines are represented by 'B 1 ', 'B 2 ', 'P', 'A I ' and 'A E ', respectively.(C) The final network model that abstracts the key structure of the network in (B).By merging a multi-species compartment (gray box in (B) including 'A I ', 'P' and 'A E ') into a compartment 'A E ' (gray box in (C)), we get a simpler model in (C).

Figure 3 .
Figure 3. Characterization of the protease activation and immune response in the M1 model.Circles in the lower panels represent the steady state solutions of (11) for high and low values of a control parameter m IP , the activation rate of proteases from the transcription factor in the cell.The intersection (black circles) of the straight line (left-hand side of (11); blue solid line) and curve (right-hand side of (11); green dotted curve) corresponds to the numerical value of A * I in the upper panel.* Black filled circle = stable, empty circle = unstable.(A) The bacterial persistence phenotype in response to the weak activation of proteases (m IP = 8.0).There exists a unique positive stable equilibrium E * L in the region (pink box) where transcription factor activities are reduced and bacterial persistence is observed.(B) Suppression of bacterial growth by immune activities in response to enhanced activation of proteases (m IP = 50).There exist three positive equilibria E* L , E * U and E * H simultaneously.While two equilibria E * L and E * U are unstable, the equilibrium E * H is stable.The equilibrium E * H resides in the region (blue box) where bacterial activities are reduced and high internal immune responses persist.All other parameters are fixed as in Table 1.

Figure 4 .
Figure 4. Bifurcation curves on the M1 model (2).(A) The B-m IP hysteresis loop: bacterial growth is active when m IP varies in the upper stable branch and suppressed when m IP varies in the lower stable branch.We define the bacterial persistence types by B > th B and the immune boosting region by B < th B and take th B = 0.7.As m IP is increased from a low value (black arrow) in the upper branch, the system loses stability at a Hopf bifurcation point (black arrowhead, marked with 'H').(B-D)The corresponding hysteresis loops for intracellular transcription factors (A I ), protease level (P) and extracellular cytokines (A E ).Other parameters are fixed as in Table1.Black = stable, blue = unstable, red dots = Hopf bifurcation point.

Figure 5 .
Figure 5. Dynamics of the M1 model (2).(A-C) Dynamics of the M1 system in the B-A E phase plane to a low (m IP = 30; (A)), intermediate (m IP = 50; (B)) and high (m IP = 80; (C)) activation rate (m IP ) of protease by the transcription factor.* Filled red circles in (A-C) = stable steady state (S.S.), empty black circle in (B) = unstable S.S. (D) A schematic of two adaptive types of bacterial infection systems (bacteria persist (T B ) and immune boosting (T I ) systems): T B = {(B, A E ) ∈ R 2 : B > th B , 0 ≤ A E < th AE },

Figure 6 .
Figure 6.Co-existence and dynamics of harmful and good bacteria in response to various cytokine clearance levels (δ) in the competition M2 model (3)-(5).(A-C) Trajectories (B 1 (t), B 2 (t)) of bacterial populations for various decay rates of cytokines (δ = 0.033 (A), 0.35 (B), 0.7 (C)) with four initial conditions; (D) characterization of the dynamical system the B 1 -B 2 plane.There exist three phenotypic regions: harmful bacteria-persist (T B , lower pink box near x-axis), harmful bacteria-free (T F , left blue box near y-axis) and co-existence (T c , gray box in the center) regions.As δ is increased (δ = 0.0033 → 0.35 → 0.7), the system undergoes the transition from T B to T c and then to T F .(E-G) Trajectories of (B 1 (t), A E (t)) corresponding to (A-C), respectively.All other parameters are fixed as in Table 4 (Type I).

Figure 7 .
Figure 7. (A) Steady state solutions of bad bacteria (B * 1 ) as a function of δ corresponding to E 1 (red), E 2 (blue) and E + (black) equilibria, respectively.Solid curve = stable, dotted curve = unstable.Green arrows indicate the points δ = 0.0033, 0.35, 0.7 corresponding to Figure 6A-C, respectively.(B) Trajectories of B 1 (t) in a B 1 − δ plane when δ is a monotonic increasing function of time, satisfying dδ dt = 0.00002 with the initial condition B 1 (0) = 0.12, B 2 (0) = 0.01, A E (0) = 0, δ(0) = 0.01 near the stable equilibrium point E 1 = E 1 (δ = 0.01).As δ is increased, the system sequentially follows the stable branches in (A), leading to the transition from stable E 1 branch (red solid curve in (A)) to stable E + branch (black solid curve in (A)) and then to stable E 2 branch (blue solid curve in (A)).All other parameters are fixed as in Table4(Type I).

Figure 8 .
Figure 8. Bi-stability dynamics of harmful and beneficial bacteria in response to the immune system in the two-species M2 model (3)-(5).(A) Dynamics of harmful bacteria-persist (T B ) and bacteria-free (T F ) in the B 1 -B 2 phase-plane.While initial states (B 1 , B 2 ) in the region R B (red) converge to the harmful bacteria-persist equilibrium (T B ), the initial states in the upper-left region (R F ; blue) converge to the harmful bacteria-free equilibrium (T F ). Blue and red curves indicate the trajectories (B 1 (t), B 2 (t)) for two very close initial conditions: B 1 (0) = 0.15, B 2 (0) = 0.713 (blue curve) and B 1 (0) = 0.15, B 2 (0) = 0.71 (red curve) near the asterisk on the boundary (green dotted curve) between R B and R F .A E (0) = 0. * Filled red circle in (A,B) = stable S.S.: E 1 = (0.9375, 0), E 2 = (0, 0.3158).(B,C) Trajectories (B 1 (t), A E (t)) and time courses (C) of bacterial populations (B 1 , B 2 ) and cytokine level (A E ) for two very close initial conditions in (A): B 1 (0) = 0.15, B 2 (0) = 0.713 (dotted curves) and B 1 (0) = 0.15, B 2 (0) = 0.71 (solid curves).(D) Characterization of the system: the dynamics adapts to the bi-stability system where the dynamical system chooses either harmful or good bacteria persisting tissue based on the initial exposure to those bacterial kinds.All other parameters are fixed as in Table 4 (Type II).

Figure 9 .
Figure 9.Effect of small time delays (τ 1 = τ 2 = 1.5) in the competition M2 model (41)-(43).(A-C) Time courses of the populations of bad bacteria (B 1 in (A)) and good bacteria (B 2 in (B)) and immune response (A E in (C)) in the absence (blue solid lines) and presence (red dashed lines) of time delays; (D) the corresponding trajectories of B 1 and B 2 in (A,B) in the phase plane; (E) the corresponding trajectories of B 1 andA E in (B,C) in the phase plane.The introduction of weak time delays in the system induces a switch from the B2-dominant equilibrium (0, B s,2 2 , A s,2 E ) to the B1-dominant equilibrium (B s,1 1 , 0, A s,1 E ).Initial condition: B 1 (0) = 0.15, B 2 (0) = 0.72, A E (0) = 0.0.All other parameters are fixed as in Table4(Type II).

Figure 10 .
Figure 10.Effect of large time delays (τ 1 = τ 2 = 3.5) in the competition M2 model (41)-(43).The introduction of time delays in the system induces oscillatory behaviors of both bad bacteria (B 1 ) and immune cytokines (A E ) and the extinction of good bacteria.(A-C) Time courses of the populations of bad bacteria (B 1 in (A)) and good bacteria (B 2 in (B)) and the immune response (A E in (C)) in the absence (blue solid lines) and presence (red dashed lines) of time delays; (D) the corresponding trajectories of B 1 and B 2 in (A,B) in the phase plane; (E) the corresponding trajectories of B 1 and A E in (B,C) in the phase plane; Initial condition: B 1 (0) = 0.15, B 2 (0) = 0.72, A E (0) = 0.0.All other parameters are fixed as in Table4(Type II).

Figure 13 .
Figure13.Dynamics of the system in response to low and high doses of drugs in the two-species M2 model (52)-(56).Immune activity was enhanced by the injection of drugs in a periodic fashion with an injection period 24 h (τ D = 24 h).(A) Time courses of populations of bad (B 1 ) and good bacteria (B 2 ) in response to periodic injection of drugs with a lower infusion rate λ D = 0.2 (duration t d = 1 h fixed); (B) time courses of immune response (A E ) and drug levels (D) corresponding to (A); (C) time courses of populations of bad (B 1 ) and good bacteria (B 2 ) in response to periodic injection of drugs with higher infusion rate λ D = 1.0 (duration t d = 1 h); (D) time courses of immune response (A E ) and drug levels (D) corresponding to (C).Initial condition:B 1 (0) = 0.2, B 2 (0) = 0.95, A E (0) = 0.0, D(0) = 0.Parameter values: N D = 10, δ D = 0.1, τ 1 = τ 2 = 3.5; all other parameters are fixed as in Table4(Type II).

B 1 (B 1 Figure B1 .
Figure B1.Sensitivity analysis of a mathematical model, Equations (3)-(5).General Latin hypercube sampling (LHS) scheme and partial rank correlation coefficient (PRCC) performed on the model.The reference outputs are the densities of bad bacteria (B 1 ) and good bacteria (B 2 ) and the concentration of inflammatory cytokines (A E ) at time t = 1, 40, 80.The colors indicate PRCC values of all variables (B 1 , B 2 , A E ) in the simple model (Equations (3)-(5)) for ten model parameters (r 1 , K 1 , α 12 , r 2 , α 21 , K 2 , γ, β 1 , β 2 , δ).Red color indicates positive correlations, and blue color indicates negative correlations between the main variable and each parameter at the given time.The analysis was carried out using the method ofMarino et al. (2008) [66] with a sample size of 10,000.

Table 1 .
Dimensionless parameter values in the M1 model.TF = transcription factor.
An internal equilibrium exists if and only if:

Table 4 .
Parameters used in the M2 model.
Table A1 lists the reference values in the model.

Table A1 .
Reference value.tw = estimated in this work.