# A Skin Microbiome Model with AMP interactions and Analysis of Quasi-Stability vs Stability in Population Dynamics

Eléa Thibault Greugny, François Fages, Ovidiu Radulescu, Peter Szmolyan, Georgios Stamatas

► To cite this version:

Eléa Thibault Greugny, François Fages, Ovidiu Radulescu, Peter Szmolyan, Georgios Stamatas. A Skin Microbiome Model with AMP interactions and Analysis of Quasi-Stability vs Stability in Population Dynamics. Theoretical Computer Science, 2023, pp.114294. 10.1016/j.tcs.2023.114294. hal-04253028

HAL Id: hal-04253028

https://inria.hal.science/hal-04253028v1

Submitted on 21 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
A Skin Microbiome Model with AMP interactions and
Analysis of Quasi-Stability vs Stability in Population
Dynamics

Eléa Thibault Greugny<sup>a,b</sup>, François Fages<sup>b</sup>, Ovidiu Radulescu<sup>c</sup>, Peter

Szmolyan<sup>d</sup>, Georgios N. Stamatas<sup>a</sup>

<sup>a</sup> Johnson & Johnson Santé Beauté France, Issy-les-Moulineaux, France

<sup>b</sup>Inria Saclay, Lifeware Team, Palaiseau, France

<sup>c</sup> University of Montpellier, France

<sup>d</sup>Technical University, Wien, Austria

Abstract

The skin microbiome plays an important role in the maintenance of a healthy
skin. It is an ecosystem, composed of several species, competing for resources
and interacting with the skin cells. Imbalance in the cutaneous microbiome, also
called dysbiosis, has been correlated with several skin conditions, including acne
and atopic dermatitis. Generally, dysbiosis is linked to colonization of the skin
by a population of opportunistic pathogenic bacteria. Treatments consisting in
non-specific elimination of cutaneous microflora have shown conflicting results.
In this article, we introduce a mathematical model based on ordinary differ-
ential equations, with 2 types of bacteria populations (skin commensals and
opportunistic pathogens) and including the production of antimicrobial pep-
tides to study the mechanisms driving the dominance of one population over
the other. By using published experimental data, assumed to correspond to the
observation of stable states in our model, we reduce the number of parameters
of the model from 13 to 5. We then use a formal specification in quantitative
temporal logic to calibrate our model by global parameter optimization and
perform sensitivity analyses. On the time scale of 2 days of the experiments,
the model predicts that certain changes of the environment, like the elevation
of skin surface pH, create favorable conditions for the emergence and coloniza-
tion of the skin by the opportunistic pathogen population, while the production
of human AMPs has non-linear effect on the balance between pathogens and
commensals. Surprisingly, simulations on longer time scales reveal that the
equilibrium reached around 2 days can in fact be a quasi-stable state followed
by the reaching of a reversed stable state after 12 days or more. We analyse
the conditions of quasi-stability observed in this model using tropical algebraic
methods, and show their non-generic character in contrast to slow-fast systems.
These conditions are then generalized to a large class of population dynamics
models over any number of species.

Keywords: skin microbiome; atopic dermatitis; antimicrobial peptides;
population dynamics; quasi-stability; meta-stability; ODE models; model reduction; steady-state reasoning; quantitative temporal logics; sensitivity analyses; tropical algebra.

# 1. Introduction

Located at the interface between the organism and the surrounding environment, the skin constitutes the first line of defense against external threats, including irritants and pathogens. In order to control potential colonization of the skin surface by pathogens, the epidermal cells, called keratinocytes, produce antimicrobial peptides (AMPs) [44]. The physiologically acidic skin surface pH also contributes to control the growth of bacterial populations [45, 28]. Another contributor to the defense against pathogen colonization are commensal bacteria in the community of microorganisms living on the skin, commonly referred to as the skin microbiome. Over the past decade, several studies have highlighted the key role played by such commensal bacterial species defending against invading pathogens, as well as their contribution to the regulation of the immune system [31, 9, 30, 26, 3, 5].

Alterations in the composition of the skin microbiome resulting in a dominance by a pathogenic species, also called dysbiosis, have been associated with skin conditions such as acne or atopic dermatitis (AD) [32, 27]. In the case of AD, the patient skin is often colonized by *Staphylococcus aureus* (*S. aureus*), especially on the lesions [27], while in the case of acne it is colonized by another bacteria, *C. acnes*. Treatment strategies targeting non-specific elimination of cutaneous microflora, such as bleach baths, have shown conflicting results regarding their capacity to reduce the disease severity [8]. On the other hand, treatments involving introduction of commensal species, like *Staphylococcus hominis* [40] on the skin surface appear promising. Accordingly, the interactions between the commensal populations, pathogens and skin cells seem at the heart of maintaining microbiome balance. There is therefore a necessity to investigate further those interactions and the drivers of dominance of one population over others. Unfortunately, it is challenging to perform *in vitro* experiments involving more than one or two different species, even more so on skin explants or skin equivalents.

Mathematical models of population dynamics have been developed and used for more than 200 years [34]. Here<fnref n="1" />, we introduce a model based on ordinary differential equations (ODEs), describing the interactions of a population of

<footnote n="1">¹This article is an extended version of a previous communication at CMSB 2022 published in [19]. The equilibrium constraint, formalized in quantitative temporal logic and used here for parameter search and sensitivity analyses, is a stronger property than in the previous version. Results of global parameter optimization and sensitivity analyses have now been added for all parameters. The main addition however is the mathematical analysis of the quasi-stability phenomenon exhibited by the model, the analysis of its non-generic character, in contrast to well-studied slow-fast systems, and its generalization to a large class of population dynamics models over any number of species in Sec. 7.</footnote>
commensal species with one of opportunistic pathogens and the skin cells. To our knowledge, there are only two other mathematical models describing the interactions between different bacteria populations and the skin [38, 36]. Nakaoka et al. introduced a model combining ODE and delayed differential equations to look at the competition dynamic between 2 populations of bacteria exposed to cytokines. Miyano et al. designed a QSP model including the interactions between *S. aureus* and coagulase negative Staphylococcus to evaluate AD treatment strategies targeting specifically *S. aureus*. To the best of our knowledge, the model we present here is the first to consider explicitly the production of antimicrobial peptides by skin commensal bacteria. This allows us to study the factors influencing the dominance of one population over the other, especially the elevation of pH and the production of antimicrobial peptides, in addition to the respective growth rates of the bacterial populations.

More specifically, based on published experimental data of Nakatsuji et al. [39] and Kohda et al. [25], assumed to correspond to steady states of our parametric ODE model, we derive constraints on the parameter values which allow us to reduce the parametric dimension of our model from 13 to 5 parameters. The landscape of those 5 remaining parameter values is then explored using global parameter optimization and sensitivity analyses procedures, implemented in our modeling software Biocham [6], using a specification of the properties of the different bacterial population equilibria formalized in quantitative temporal logic FO-LTL($\mathbb{R}_{lin}$) [50].

On the time scale of the experiments of Kohda et al. [25], our model predicts that certain changes in the environment, like an elevation of skin surface pH, create favorable conditions for the emergence and colonization of the skin by the opportunistic pathogen population. The model also predicts that the concentration of human AMPs affects in a non-linear manner the balance between pathogens and commensals. Such predictions can help identify potential therapeutic strategies for the treatment of skin conditions involving microbiome dysbiosis, or optimize AMP levels on the skin surface in topical treatments.

Interestingly, in the reduced model, we observe an unexpected phenomenon of meta-stability [4, 56, 47, 53, 52, 14] or quasi-stability [37], a terminology we prefer to adopt in the absence of stochasticity, in which the seemingly stable state observed at the time scale of the experiments of 2 days is followed by a switch after 12 days toward a reversed stable state. We develop a complete mathematical analysis of the conditions of quasi-stability observed in our population dynamics model using singular perturbations and tropical algebraic methods. We show the non-generic character of this form of quasi-stability in our model, in sharp contrast to slow-fast systems which are mainly studied. Finally, we generalize the conditions first obtained by tropical analysis of the parameter values of our model, to a large class of population dynamics models over any number of species which will similarly exhibit quasi-stability phenomenon on relatively long timescales.
<figure><img src="image_2.png" /><figcaption>Figure 1: Model overview of the interactions between the commensal and pathogenic bacterial populations at the surface of the skin, including the production of bacterial and human AMPs. Green arrows represent positive production effects and red T-lines represent negative killing effect.</figcaption></figure>

## 2. Initial ODE model with 13 parameters

The ODE model built in this article considers two types of bacterial populations. The first population, $S_c$, regroups commensal bacteria species having an overall beneficial effect for the skin, and the second population, $S_p$, represents opportunistic pathogens. One originality of our model is that it considers their interactions through the production of both human and bacterial AMPs, as depicted in Fig. 1.

The terms of the differential equations governing the growth of both bacterial populations are based on the common logistic growth model [59], considering non-explicitly the limitations in food and space. The limited resources are included in the parameters $K_{sc}$ and $K_{sp}$, representing the optimum concentration of the populations in a given environment, considering the available resources. $r_{sc}$ and $r_{sp}$ represent the growth rate of commensal and pathogenic bacteria respectively.

The bactericidal effect of antimicrobial peptides (AMPs) produced by skin cells, $Amp_h$, on $S_p$ is included with a Hill function of order 1 representing a simple saturation effect. Parameter $d_{sph}$ represents the maximal killing rate of $S_p$ by $Amp_h$ and $C_{ah}$ is the concentration of $Amp_h$ inducing half this maximum killing rate. This type of highly non-linear functions is commonly used to model the saturating effect of antibiotics on bacterial populations [7, 35, 55]. For the sake of simplicity, the AMPs produced by skin cells is introduced as a constant parameter, [$Amp_h$], in the model. It represents the average concentration of these AMPs among surface cells, under given human genetic background and environmental conditions.
Several studies revealed that commensal bacterial populations, like *S. epidermidis* or *S. hominis*, are also able to produce AMPs targeted against opportunistic pathogens, such as *S. aureus* [10, 39]. For these reasons, we chose to also introduce in our model AMPs of bacterial origin, *Ampb*, acting similarly to *Amph* on the pathogenic population *S*<sub>p</sub>. *Ampb* is produced at rate *k*<sub>c</sub> by *S*<sub>c</sub>, and degraded at rate *d*<sub>a</sub>. Furthermore, we include a defense mechanism of *S*<sub>p</sub> against *S*<sub>c</sub> with a direct killing effect. This killing effect is also modelled with a Hill function or order 1, where *d*<sub>sc</sub> is the maximum killing rate of *S*<sub>c</sub> by *S*<sub>p</sub> and *C*<sub>1</sub> the concentration of *S*<sub>p</sub> necessary to reach half this maximum killing rate.

Altogether, this gives us the following ODE system with 3 variables and 13 parameters, all taking non-negative values:

$$
\left\{
\begin{aligned}
\frac{d[S_c]}{dt} &= \left( r_{sc} \left( 1 - \frac{[S_c]}{K_{sc}} \right) - \frac{d_{sc}[S_p]}{C_1 + [S_p]} \right) [S_c] \\
\frac{d[S_p]}{dt} &= \left( r_{sp} \left( 1 - \frac{[S_p]}{K_{sp}} \right) - \frac{d_{spb}[Amp_b]}{C_{ab} + [Amp_b]} - \frac{d_{sph}[Amp_h]}{C_{ah} + [Amp_h]} \right) [S_p] \\
\frac{d[Amp_b]}{dt} &= k_c[S_c] - d_a[Amp_b]
\end{aligned}
\right.
\quad (1)
$$

Table 1 recapitulates the 3 variables and 13 parameters of this model with their units.

### 3. Using published experimental data to eliminate model parameters by steady-state reasoning

The amount of quantitative experimental data available for the model calibration is very limited due to the difficulty of carrying out experiments involving co-cultures of different bacterial species. Most of the published work focuses on single species or on measuring the relative abundances of species living on the skin, which is highly variable between individuals and skin sites [20]. In the case of AD specifically, *S. aureus* is considered pathogenic and *S. epidermidis* commensal. Published data exist however for those species which we can use to constrain the parameter values of the model.

Two series of *in vitro* experiments are considered [39, 25]. While *in vitro* cultures, even on epidermal equivalent, do not entirely capture the native growth of bacteria on human skin, they provide useful quantitative data that would be very difficult to measure *in vivo*.

In the first experiment [25], mono-cultures and co-cultures of *S. epidermidis* and *S. aureus* were allowed to develop on a 3D epidermal equivalent. Table 2 recapitulates the population sizes of the two species measured after 48 hours of incubation. Kohda *et al.* also performed another co-culture experiment where *S.epidermidis* was inoculated 4 hours prior to *S.aureus* in the media. This data
<table><caption>Table 1: List of the parameters and variables of our mathematical model with their units. CFU = Colony forming unit, AU = Arbitrary Unit, ASU = Arbitrary Surface Unit</caption><thead><tr><th>Variable</th><th>Interpretation (unit)</th></tr></thead><tbody><tr><td>[S<sub>c</sub>]</td><td>Surface apparent concentration of S<sub>c</sub> (CFU.ASU<sup>-1</sup>)</td></tr><tr><td>[S<sub>p</sub>]</td><td>Surface apparent concentration of S<sub>p</sub> (CFU.ASU<sup>-1</sup>)</td></tr><tr><td>[Amp<sub>b</sub>]</td><td>Concentration of Amp<sub>b</sub> (AU.ASU<sup>-1</sup>)</td></tr><tr><th>Parameter</th><th>Interpretation (unit)</th></tr><tr><td>r<sub>sc</sub></td><td>Growth rate of S<sub>c</sub> (h<sup>-1</sup>)</td></tr><tr><td>r<sub>sp</sub></td><td>Growth rate of S<sub>p</sub>, (h<sup>-1</sup>)</td></tr><tr><td>K<sub>sc</sub></td><td>Optimum concentration of S<sub>c</sub> (CFU.ASU<sup>-1</sup>)</td></tr><tr><td>K<sub>sp</sub></td><td>Optimum concentration of S<sub>p</sub> (CFU.ASU<sup>-1</sup>)</td></tr><tr><td>d<sub>sc</sub></td><td>Maximal killing rate of S<sub>c</sub> by S<sub>p</sub> (h<sup>-1</sup>)</td></tr><tr><td>C<sub>1</sub></td><td>Concentration of S<sub>p</sub> inducing half the maximum killing rate d<sub>sc</sub> (CFU.ASU<sup>-1</sup>)</td></tr><tr><td>d<sub>spb</sub></td><td>Maximal killing rate of S<sub>p</sub> by Amp<sub>b</sub>, (h<sup>-1</sup>)</td></tr><tr><td>C<sub>ab</sub></td><td>Concentration of Amp<sub>b</sub> inducing half the maximum killing rate d<sub>spb</sub> (AU.ASU<sup>-1</sup>)</td></tr><tr><td>d<sub>sph</sub></td><td>Maximal killing rate of S<sub>p</sub> by Amp<sub>b</sub>, (h<sup>-1</sup>)</td></tr><tr><td>C<sub>ah</sub></td><td>Concentration of Amp<sub>b</sub> inducing half the maximum killing rate d<sub>sph</sub> (AU.ASU<sup>-1</sup>)</td></tr><tr><td>[Amp<sub>b</sub>]</td><td>Concentration of AMPs produced by the skin cells (AU.ASU<sup>-1</sup>)</td></tr><tr><td>k<sub>c</sub></td><td>Production rate of Amp<sub>b</sub> by S<sub>c</sub> (AU.h<sup>-1</sup>.CFU<sup>-1</sup>)</td></tr><tr><td>d<sub>a</sub></td><td>Degradation rate of Amp<sub>b</sub> (AU.h<sup>-1</sup>)</td></tr></tbody></table>
<table><caption>Table 2: Experimental data from Kohda et. al [25] used for identifying parameter values.</caption><thead><tr><th></th><th><i>S. epidermidis</i> (CFU/well)</th><th><i>S. aureus</i> (CFU/well)</th></tr></thead><tbody><tr><td>Mono-cultures</td><td>4.10<sup>8</sup></td><td>3.10<sup>9</sup></td></tr><tr><td>Co-cultures</td><td>1.10<sup>8</sup></td><td>1.10<sup>9</sup></td></tr></tbody></table>

<table><caption>Table 3: Experimental data from Nakatsuji et. al [39] used for identifying parameter relations.</caption><thead><tr><th><i>Sh</i>-lantibiotics (μM)</th><th>LL-37 (μM)</th><th><i>S. aureus</i> (CFU/mL)</th></tr></thead><tbody><tr><td>0</td><td>4</td><td>10<sup>9</sup></td></tr><tr><td>0</td><td>8</td><td>6.10<sup>5</sup></td></tr><tr><td>0.32</td><td>0</td><td>5.10<sup>8</sup></td></tr><tr><td>0.64</td><td>0</td><td>3.10<sup>3</sup></td></tr></tbody></table>

is not used here as it requires additional manipulation to match the situation represented by the model. However, it would be interesting to use it in the future for model validation.

In the second experiment [39] the impact of human (LL-37) and bacterial (*Sh*-lantibiotics) AMPs on *S. aureus* survival was studied. The experiments were performed *in vitro*, and the *S. aureus* population size was measured after 24 hours of incubation. Table 3 summarizes their observations.

### 3.1. Parameter values inferred from mono-culture experiment data

We consider first the monocultures experiments from Kohda et al. [25], representing the simplest experimental conditions. *S. epidermidis* is a representative of the commensal population *S*<sub>c</sub>, and *S. aureus* of the pathogenic one, *S*<sub>p</sub>. Since the two species are not interacting, the set of equations simplifies to:

$$ 
\left\{
\begin{aligned}
\frac{d[S_c]}{dt} &= \left( r_{sc} \left( 1 - \frac{[S_c]}{K_{sc}} \right) \right) [S_c] \\
\frac{d[S_p]}{dt} &= \left( r_{sp} \left( 1 - \frac{[S_p]}{K_{sp}} \right) \right) [S_p]
\end{aligned}
\right.
\qquad (2)
$$

At steady-state, the population concentrations are either zero, or equal to their optimum capacities (*K*<sub>sc</sub> or *K*<sub>sp</sub>) when the initial population concentration is non-zero. Given the rapid growth of bacterial population, the experimental
measurements done after 48 hours of incubation can be considered as corresponding to a steady-state, which gives:

$$K_{sc} = 4.10^8 \text{ CFU.ASU}^{-1} \qquad (3)$$

$$K_{sp} = 3.10^9 \text{ CFU.ASU}^{-1} \qquad (4)$$

### 3.2. Parameter relations inferred from experimental data on AMP

The experimental conditions of Nakatsuji et al. [39] correspond to the special case where there is no commensal bacteria alive in the environment, only the bacterial AMPs, in addition to those produced by the skin cells. Our system of equations then reduces to:

$$\frac{d[S_p]}{dt} = \left( r_{sp} \left( 1 - \frac{[S_p]}{K_{sp}} \right) - \frac{d_{spb}[Amp_b]}{C_{ab} + [Amp_b]} - \frac{d_{sph}[Amp_h]}{C_{ah} + [Amp_h]} \right) [S_p] \quad (5)$$

The concentrations in LL-37 and Sh-lantibiotics, translated in our model into $[Amp_h]$ and $[Amp_b]$ respectively, are part of the experimental settings. Therefore, we consider them as constants over time. At steady state, we get:

$$[S_p]^* = 0 \quad \text{or} \quad [S_p]^* = K_{sp} \left( 1 - \frac{d_{spb}[Amp_b]}{r_{sp}(C_{ab} + [Amp_b])} - \frac{d_{sph}[Amp_h]}{r_{sp}(C_{ah} + [Amp_h])} \right) \quad (6)$$

Let us first focus on the special case where no Sh-lantibiotics were introduced in the media, translating into $[Amp_b] = 0$ in our model. We consider again that the biological observations after 24 hours of incubation correspond to steady-state and substitute the experimental values measured

$[Amp_h] = 4 \mu M$; $[S_p]^* = 10^9 \text{ CFU}$, and $[Amp_h] = 8 \mu M$; $[S_p]^* = 6.10^5 \text{ CFU}$, together with the values of $K_{sc}$ and $K_{sp}$ (from (3) and (4)) in (6), to obtain the following equations:

$$\left\{ \begin{aligned} \frac{d_{sph}}{r_{sp}} &= \frac{4 + C_{ah}}{6} \\ \frac{d_{sph}}{r_{sp}} &= \frac{(10^4 - 2)(C_{ah} + 8)}{8.10^4} \end{aligned} \right. \qquad (7)$$

which reduce to $C_{ah} = 8$ and $\frac{d_{sph}}{r_{sp}} = 2$.

Following the same method with the experimental conditions without any LL-37 (i.e. $[Amp_h] = 0$) and using two data points

([Amp_b] = 0.32 \mu M$; $[S_p]^* = 5.10^8 \text{ CFU}$)

and ([Amp_b] = 0.64 \mu M$; $[S_p]^* = 3.10^3 \text{ CFU}$),

we get $C_{ab} = 0.16$ and $\frac{d_{spb}}{r_{sp}} = \frac{5}{4}$.
It is notable that the maximum killing rates of $S_p$ by $Amp_b$ and $Amp_h$ are both proportional to $S_p$ growth rate. Interestingly, such proportional relation has been observed experimentally between the killing rate of *Escherichia coli* by an antibiotic and the bacterial growth rate [57].

To be consistent with the ranges of Sh-lantibiotics concentrations described in Nakatsuji et al. [39], $[Amp_b]$ should take positive values below 10. Given that $[Amp_b]^* = \frac{k_c[S_c]^*}{d_a}$ at steady-state, and that $K_{sc} = 4.10^8$ CFU is the upper bound for $[S_c]^*$, we obtain the following constraint:

$$ \frac{k_c}{d_a} \le \frac{1}{4.10^7} \qquad (8) $$

### 3.3. Parameter relations inferred from co-culture data

The initial model described earlier is representative of the experimental settings of the co-culture conditions described in Kohda et al. [25]. At steady-state, the system (1) gives:

$$ [S_c]^* = 0 \quad \text{or} \quad [S_c]^* = K_{sc} \left( 1 - \frac{d_{sc}[S_p]^*}{r_{sc}(C_1 + [S_p]^*)} \right) \qquad (9) $$

$$ [S_p]^* = 0 \quad \text{or} \quad [S_p]^* = K_{sp} \left( 1 - \frac{d_{spb}[Amp_b]}{r_{sp}(C_{ab} + [Amp_b])} - \frac{d_{sph}[Amp_h]}{r_{sp}(C_{ah} + [Amp_h])} \right) \qquad (10) $$

$$ [Amp_b]^* = \frac{k_c[S_c]^*}{d_a} \qquad (11) $$

Considering that what is observed experimentally after 48 hours of incubation is at steady-state, one can replace $[S_c]^*$ and $[S_p]^*$ with the experimental data point ($S. epidermidis$ = $10^8$ CFU; $S. aureus$ = $10^9$ CFU) in (9) and (10) to get the following parameter relation:

$$ \frac{d_{sc}}{r_{sc}} = \frac{3}{4.10^9} C_1 + \frac{3}{4} \qquad (12) $$

$$ \frac{2}{3} r_{sp} = \frac{d_{sph}[Amp_h]}{C_{ah} + [Amp_h]} + \frac{10^8 d_{spb} k_c}{d_a C_{ab} + 10^8 k_c} \qquad (13) $$

By integrating the values found for $C_{ah}$ and $C_{ab}$, and the relations involving $d_{sph}$ and $d_{spb}$ into (13), we end up with:

$$ d_a = 10^8 k_c \frac{56 + 31[Amp_h]}{2.56(4 - [Amp_h])} \quad \text{with } [Amp_h] < 4 \qquad (14) $$
<table><caption>Table 4: Summary of the parameter relations embedded in the reduced model.</caption><thead><tr><th>Parameter</th><th>Value or relation to other parameters</th></tr></thead><tbody><tr><td>K<sub>sc</sub></td><td>4.10<sup>8</sup></td></tr><tr><td>K<sub>sp</sub></td><td>3.10<sup>9</sup></td></tr><tr><td>C<sub>ah</sub></td><td>8</td></tr><tr><td>C<sub>ab</sub></td><td>0.16</td></tr><tr><td>d<sub>sph</sub></td><td>2 r<sub>sp</sub></td></tr><tr><td>d<sub>spb</sub></td><td><sup>5</sup>&frasl;<sub>4</sub> r<sub>sp</sub></td></tr><tr><td>d<sub>sc</sub></td><td>r<sub>sc</sub> ( <sup>3</sup>&frasl;<sub>4.10<sup>9</sup></sub> C<sub>1</sub> + <sup>3</sup>&frasl;<sub>4</sub> )</td></tr><tr><td>d<sub>a</sub></td><td>10<sup>8</sup> k<sub>c</sub> <sup>56 + 31 [Amp<sub>h</sub>]</sup>&frasl;<sub>2.56 (4 - [Amp<sub>h</sub>])</sub> with [Amp<sub>h</sub>] &lt; 4</td></tr></tbody></table>

## 4. Reduced model with 5 parameters

Using the previously mentioned experimental data, and assuming they represent steady state conditions of the initial model (1), we have reduced the parametric dimension of the model from 13 to 5. Specifically, out of the original 13 parameters, we could define the values of 4 of them, and derive 4 functional dependencies from the values of the remaining parameters, as summarized in Table 4).

In our reduced skin microbiome model with AMPs (1), the parameters that remain unknown are thus:

* $r_{sc}$, the growth rate of $S_c$ which can reasonably take values between 0 and $2 h^{-1}$ following [11, 7];

* $r_{sp}$, the growth rate of $S_p$, taking similar values in the interval between 0 and $2 h^{-1}$;

* $C_1$, the concentration of $S_p$ that induces half the maximum killing rate $d_{sc}$ (in CFU.ASU$^{-1}$) and is thus bounded by the optimum concentration of $S_p$, i.e. $K_{sp} = 3.10^9$ CFU.ASU$^{-1}$, as calculated in section 3.1 from [25];

* $k_c$, the production rate of $[Amp_b]$ chosen to take values between 0 and $0.1 AU.h^{-1}.CFU^{-1}$, and shown to have a limited impact on the steady-state values in section 4.2;

* $[Amp_h]$, the concentration in AU.ASU$^{-1}$ of AMPs produced by skin cells between 0 and 4 (equation (14)).
This gives us the following reduced system:

$$
\left\{
\begin{aligned}
\frac{d[S_c]}{dt} &= r_{sc}[S_c] \left( 1 - \frac{[S_c]}{4.10^8} - \frac{3(10^{-9}C_1 + 1)[S_p]}{C_1 + [S_p]} \right) \\
\frac{d[S_p]}{dt} &= r_{sp}[S_p] \left( 1 - \frac{[S_p]}{3.10^9} - \frac{5[Amp_b]}{0.64 + 4[Amp_b]} - \frac{2[Amp_h]}{8 + [Amp_h]} \right) \\
\frac{d[Amp_b]}{dt} &= k_c \left( [S_c] - 10^8[Amp_b] \frac{56 + 31[Amp_h]}{2.56(4 - [Amp_h])} \right)
\end{aligned}
\right.
\quad (15)
$$

A continuum of different parameter values exists however to fit either the experimental data of the pathogenic microbiome reported in [25], or what can be considered as a healthy microbiome phenotype. The landscape of those parameter values can nevertheless be explored using some formalization of the desired behaviors and global parameter optimization methods.

4.1. Formal specification of dynamical behaviors in quantitative temporal logic

Quantitative temporal logics, such as Metric Temporal Logic MTL [1], Signal Temporal Logic [33], or First-Order Linear Time Logic with real-valued linear constraints FO-LTL($\mathbb{R}_{\text{lin}}$) [50, 17], provide very expressive formal languages to specify behaviors of continuous time dynamical systems [15].

In particular, the definition of a continuous degree of satisfaction in the in-
terval [0, 1] for FO-LTL($\mathbb{R}_{\text{lin}}$) formulae [49] is key to achieve several purposes:
first, to guide search in continuous optimization techniques to fit parameter
values to satisfy model's temporal specification; second, to evaluate the robust-
ness of the model with respect to its specification and parameter variations; and
third, to compute parameter sensitivity indices to determine the most influential
parameters.

This is the method used here to fit our reduced model to the pathogenic phe-
notype reported in Kohda et al. experiments [25], and to a healthy phenotype.
We use the implementation of temporal logic FO-LTL($\mathbb{R}_{\text{lin}}$) in Biocham [50] for
this purpose and for sensitivity analyses$^2$.

For example, in the pathogenic experiments, the property of interest is the
reaching of an equilibrium around 48 hours of the bacterial populations with a
predominance of the pathogenic bacteria of the following orders: [Sp] = 10⁹ and
[Sc] = 10⁸ [25]. That first property is a simple curve fitting problem with two
time points that can be formalized in FO-LTL(R<sub>lin</sub>) by the following temporal
logic formula with free variables x1, y1, x2, y2 for the population sizes

$$
\mathbf{F}(\text{Time} \sim 40 \land \mathit{Sc} = x1 \land \mathit{Sp} = y1 \land \mathbf{F}(\mathbf{G}(\mathit{Sc} = x2 \land \mathit{Sp} = y2))) \quad (16)
$$

<footnote n="2">²All computation results reported in this article have been obtained with a Biocham notebook available in different formats at https://lifeware.inria.fr/wiki/Main/Software#TCS23.</footnote>
and objective values

$$x1 = 10^8, y1 = 10^9, x2 = 10^8, y2 = 10^9 \qquad (17)$$

(finally **F**) around time 40 and (finally **F** globally **G**) in all future time points (i.e. at the last time point of a finite trace, 50 hours in the case of the experiments and of our model simulations).

The temporal formula for fitting the model to a balanced microbiome expresses the dominance of the commensal bacteria over the pathogenic bacterial populations by at least one order of magnitude, as follows:

$$\exists c \exists p \mathbf{F}(Time \sim 40 \land S_c > k*S_p \land S_c = c \land S_p = p \land \mathbf{F}(G(r*S_p = p \land s*S_c = c))) \qquad (18)$$

with objective values

$$k = 10, r = 1, s = 1. \qquad (19)$$

The abstraction by variables of the objective values makes it possible to compute first, the validity domain of the variables that satisfy the formula, as a finite union of polytopes [18], and second, a violation degree of the formula with respect to the objective values, as the Haussdorf distance between the validity domain of the free variables and that objective point. The landscape of parameter values satisfying the formula can then be explored in this setting by global parameter optimization.

The continuous degree of satisfaction in [0, 1] of the specification is obtained as the inverse of the violation degree, and can be used as optimization criterion for parameter search [50]. In Biocham, the covariance matrix adaptive evolution strategy (CMA-ES) [21] is used to search parameter values maximising the satisfaction degree of a FO-LTL formula given with respect to the objective values given for its free variables.

The satisfaction degree of a FO-LTL formula and objectives also provides as computational measure of robustness of the model and sensitivity indices to the parameter values, by taking the mean satisfaction degree when varying the parameter values according to a normal distribution $N(\mu, \sigma^2)$ around their nominal values $\mu$, with standard deviation $\sigma = COV.\mu$ given by a coefficient of variation COV.

## 4.2. Model parameter fitting to pathogenic experiments and sensitivity analyses

In order to reproduce what was observed by Kohda et. al [25], that is the reaching of an equilibrium with a dominant pathogenic population after 40 hours, which can thus be considered as dysbiosis in our skin microbiome model, we fix the initial concentrations for $[S_c]$ and $[S_p]$ to the doses of S. epidermidis and S. aureus applied at the surface of the 3D epidermal equivalent at the beginning of the experiment (i.e. $10^5$ CFU/mL and $10^3$ CFU/mL, respectively).

The parameterization problem consists in finding values for the remaining unknown parameters, $[Amp_h]$, $r_{sc}$, $r_{sp}$, $C_1$, $k_c$, within the intervals of possible values described above, and where $r_{sp}$ can be fixed to 1 arbitrarily.
<table><caption>Table 5: Optimal parameter values satisfying the reaching of an equilibrium after 40h (temporal formula 16) with a dominance of the pathogenic bacteria by one order of magnitude (objective values 17) when fixing [$Amp_h$] = 1.5, $r_{sp}$ = 1. The parameter values retained are in bold.</caption><tr><td>$C_1$</td><td>1e3</td><td>1e4</td><td>1e5</td><td>1e6</td><td><b>5e6</b></td><td>1e7</td><td>1e8</td><td>1e9</td></tr><tr><td>$r_{sc}$</td><td>1.3</td><td>1</td><td>0.76</td><td>0.59</td><td><b>0.5</b></td><td>0.46</td><td>0.38</td><td>0.34</td></tr></table>

Parameter optimization with respect to the property of reaching an equilibrium reached at time 40h (FO-LTL formula 16) with a dominance of the pathogenic bacterial population by one order of magnitude (objective values 17) shows the existence of a continuum of solutions. Not surprisingly, the solutions found assign small values to parameters [$Amp_h$], $r_{sc}$ or $C_1$, which could be anticipated from the structure of the model. On the other hand, it is worth noting that no parameter values are found to reach the opposite equilibrium by one order of magnitude only, yet solutions are found for reaching a healthy microbiome with a dominance of the commensal bacteria over the pathogenic bacteria populations by 4 orders of magnitude<fnref n="3" />.

In order to reproduce culture data, by opposition to in vivo measures, it is reasonable to fix parameter [$Amp_h$] to a relatively low value with respect to its range of typical values in the interval [0, 4]. When fixing [$Amp_h$] = 1.5 for instance, parameter search shows the irrelevance of parameter $k_c$ and interestingly, a functional dependency of optimal parameter values $r_{sp}$ as a function of $C_1$. Table 5 summarizes the optimal values found for $C_1$ and $r_{sp}$ in this setting.

We choose accordingly the parameter set [[$Amp_h$] = 1.5, $r_{sc}$ = 0.5, $r_{sp}$ = 1, $C_1$ = 5.10⁶, $k_c$ = 0.01 for fitting those experiments. Fig. 2 shows the result of a numerical simulation of our model with this parameter set which is in accordance to the co-culture experiments of Kohda et. al and does reproduce a consistent qualitative behavior [25].

The sensitivity analysis reported in Table 6 shows the expectation of the satisfaction degree of the temporal specification when varying each parameter value one-at-a-time (OAT) using a normal distribution centered around its nominal value and some indicated coefficient of variation. This reveals that the dominance of the pathogenic population is not sensitive to the values of $C_1$ and $k_c$ which can vary by orders of magnitude, and that all model parameters can vary by 20% without affecting the qualitative behaviour. The highest sensitivity indices concern the growth rates ($r_{sp}$ and $r_{sc}$) which should remain with a COV below 20%, and the threshold concentration $C_{ah}$, followed by the optimum concentrations $K_{sc}$, $K_{sp}$ and $C_{ab}$.

<footnote n="3">³The Biocham notebook mentioned above shows the detailed results of parameter search and simulation figures on logscale.</footnote>
<table><caption>Table 6: OAT and global sensitivity indices of the pathogenic model to variations of all kinetic parameters and initial concentrations around their nominal values $r_{sp} = 1, r_{sc} = 0.5, [Amp_h] = 1.5, C_1 = 5e6, k_c = 0.01, K_{sc} = 4e8, K_{sp} = 3e9, C_{ab} = 0.16, C_{ah} = 8$, for the equilibrium property of Kohda's co-culture experiments (i.e. formula 16 with objective values 17). The most influential parameters given in decreasing order of OAT sensitivity indices with $COV = 0.2$ are $r_{sp}, r_{sc}, C_{ah}, K_{sc}, K_{sp}, C_{ab}, [Amp_h], [S_c]^0, [S_p]^0, C_1, k_c$. Global sensitivity analysis provides an indication of the robustness of the model with respect to the equilibrium property for simultaneous variations of all model parameters around their nominal values with a COV of 10%.</caption><thead><tr><th>Parameter</th><th>Nominal value</th><th>COV</th><th>Sensitivity index</th></tr></thead><tbody><tr><td>r<sub>sp</sub></td><td>1</td><td>0.2</td><td>0.271</td></tr><tr><td>r<sub>sc</sub></td><td>0.5</td><td>0.2</td><td>0.253</td></tr><tr><td>C<sub>ah</sub></td><td>8</td><td>0.2</td><td>0.249</td></tr><tr><td>K<sub>sc</sub></td><td>4.10<sup>8</sup></td><td>0.2</td><td>0.194</td></tr><tr><td>K<sub>sp</sub></td><td>3.10<sup>9</sup></td><td>0.2</td><td>0.164</td></tr><tr><td>C<sub>ab</sub></td><td>0.16</td><td>0.2</td><td>0.159</td></tr><tr><td>[Amp<sub>h</sub>]</td><td>1.5</td><td>0.2</td><td>0.065</td></tr><tr><td>[Amp<sub>h</sub>]</td><td>1.5</td><td>0.3</td><td>0.11</td></tr><tr><td>[S<sub>c</sub>]<sup>0</sup></td><td>10<sup>5</sup></td><td>0.2</td><td>0.029</td></tr><tr><td>[S<sub>c</sub>]<sup>0</sup></td><td>10<sup>5</sup></td><td>0.8</td><td>0.131</td></tr><tr><td>[S<sub>p</sub>]<sup>0</sup></td><td>10<sup>3</sup></td><td>0.2</td><td>0.026</td></tr><tr><td>[S<sub>p</sub>]<sup>0</sup></td><td>10<sup>3</sup></td><td>0.5</td><td>0.131</td></tr><tr><td>C<sub>1</sub></td><td>5.10<sup>6</sup></td><td>0.2</td><td>0.021</td></tr><tr><td>C<sub>1</sub></td><td>5.10<sup>6</sup></td><td>2</td><td>0.135</td></tr><tr><td>k<sub>c</sub></td><td>0.01</td><td>0.2</td><td>0.01</td></tr><tr><td>k<sub>c</sub></td><td>0.01</td><td>10</td><td>0.01</td></tr><tr><td>{[S<sub>p</sub>]<sup>0</sup>, [S<sub>c</sub>]<sup>0</sup>}</td><td>{10<sup>3</sup>, 10<sup>5</sup>}</td><td>0.4</td><td>0.107</td></tr><tr><td>{r<sub>sp</sub>, r<sub>sc</sub>, [Amp<sub>h</sub>], C<sub>1</sub>, k<sub>c</sub>}</td><td>”</td><td>0.1</td><td>0.206</td></tr><tr><td>set of all 11 parameters</td><td>”</td><td>0.1</td><td>0.3</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>Figure 2: Numerical simulation of the reduced ODE model over 50 hours, with initial conditions $[S_c] = 10^5$, $[S_p] = 10^3$, $[Amp_b] = 0$ and parameter values $[Amp_h] = 1.5$, $r_{sc} = 0.5$, $r_{sp} = 1$, $C_1 = 5.10^6$, $k_c = 0.01$, to fit Kohda et al. co-culture data [25] (Table 2).</figcaption></figure>

### 4.3. Parameter values fitting a balanced microbiome and sensitivity analyses

Our model can also be used to reproduce what is considered a balanced microbiome, corresponding to the commensal population being significantly more abundant than the pathogenic one [27]. Given the high inter-individual variability, one cannot define precise target values for the population sizes corresponding to a balanced microbiome. Therefore, we will consider in the following analysis that the microbiome is balanced when it stabilizes with a commensal population more abundant than the pathogenic one by at least one order of magnitude (see formula 18 with objectives 19).

Satisfying this property requires modifying some parameter values to represent a less virulent pathogenic population, closer to the physiological context, given that the experiments from Kohda et. al [25] were performed using a virulent methicillin-resistant S. aureus strain.

Parameter search shows again a continuum of solutions for the parameter values, with optimal values for $[Amp_h]$ around 2.5, but no solution leading to a dominance of the commensals by just one order of magnitude, but by 4 orders of magnitude.

Sensitivity analysis reported in Table 7 reveals that the stable dominance of the commensal population is mostly sensitive to variations of $C_{ah}$, $C_{ab}$, $K_{sc}$ and interestingly $[Amp_h]$. On the other hand, the results are not sensitive to the initial bacterial population sizes suggesting a form of perfect adaptation of the balance between both population. In fact, this property follows from the absence of multiple non-degenerate steady states in this model, a property checked in the Biocham notebook of the model using the graphical requirements method for multistationarity of [2].

In order to get closer to physiological context however, we prefer to retain
<table><caption>Table 7: OAT and global sensitivity indices of the safe microbiome model to variations of the parameters and initial concentrations around their nominal values with a COV of 20% for the property of maintaining stable populations of bacteria at times 40h-50h with a pre-dominance of the commensal bacteria by at least one order of magnitude (i.e. formula 18 with objective 19). OAT sensitivity indices show the preponderant influence of parameters $C_{ah}$, $C_{ab}$, $K_{sc}$, $[Amp_h]$, $r_{sc}$ for maintaining a healthy balance, while the other kinetic parameters and the initial conditions appear much less important. Global sensitivity indices are an indication for the robustness of the model.</caption><thead><tr><th>Parameters</th><th>Nominal value</th><th>COV</th><th>Sensitivity index</th></tr></thead><tbody><tr><td>$C_{ah}$</td><td>8</td><td>0.2</td><td>0.242</td></tr><tr><td>$C_{ab}$</td><td>0.16</td><td>0.2</td><td>0.207</td></tr><tr><td>$K_{sc}$</td><td>4.10<sup>8</sup></td><td>0.2</td><td>0.188</td></tr><tr><td>$[Amp_h]$</td><td>2.5</td><td>0.2</td><td>0.187</td></tr><tr><td>$r_{sc}$</td><td>0.5</td><td>0.2</td><td>0.007</td></tr><tr><td>$K_{sp}$</td><td>3.10<sup>9</sup></td><td>0.2</td><td>0.005</td></tr><tr><td>$[S_p]^0$</td><td>10<sup>3</sup></td><td>0.2</td><td>0.005</td></tr><tr><td>$[S_c]^0$</td><td>10<sup>5</sup></td><td>0.2</td><td>0.005</td></tr><tr><td>$r_{sp}$</td><td>0.5</td><td>0.2</td><td>0.005</td></tr><tr><td>$C_1$</td><td>2.10<sup>9</sup></td><td>0.2</td><td>0.005</td></tr><tr><td>$k_c$</td><td>0.01</td><td>0.2</td><td>0.005</td></tr><tr><td>$(r_{sp}, r_{sc}, [Amp_h], C_1, k_c)$</td><td>”</td><td>0.2</td><td>0.184</td></tr><tr><td>tuple of all 11 parameters</td><td>”</td><td>0.2</td><td>0.321</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Figure 3: Numerical simulation of the reduced ODE model over 50 hours, with initial conditions [$S_c$] = $10^5$, [$S_p$] = $10^3$, [$Amp_b$] = 0 and parameter values $r_{sc}$ = $r_{sp}$ = 0.5, $C_1$ = $2.10^8$, $k_c$ = 0.01, [$Amp_h$] = 3 retained here for a healthy microbiome in a physiological context, with a value of $Sp$ = $6.10^4$ at equilibrium, less than $Sc$ by 4 orders of magnitude.</figcaption></figure>

the slightly suboptimal parameter set with a higher production of AMPs by the skin cells, [$Amp_h$] = 3, to compensate for feedback loops or stimuli that might be missing in the 3D epidermal equivalent used, and $r_{sp}$ = $r_{sc}$ = 0.5, $C_1$ = $2.10^8$. Fig. 3 shows a simulation trace obtained under those conditions which clearly satisfies the dominance property after 10h and reaching equilibrium around 25h, as shown also in the notebook with more details.

**5. Conditions favoring the pathogenic population**

Whether the dysbiosis observed in AD is the cause or the result of the disease is unclear [23, 24]. Infants developing AD do not necessarily have more *S. aureus* present on their skin prior to the onset of the disease compared to the healthy group [22]. This suggests that atopic skin has some characteristics enabling the dominance of *S. aureus* over the other species of the microbiome. To test this hypothesis, we investigate two changes of the skin properties observed in AD patients (skin surface pH elevation [16] and reduced production of AMPs [43]), and their impact on the dominant species at steady-state.

More specifically, we study the behavior of the system following the introduction of a pathogen and whether the pathogen will colonize the media depending on the initial concentrations of the bacterial populations and the particular skin properties mentioned before.

**5.1. Skin surface pH elevation**

According to Proksch [45], the physiological range for skin surface pH is 4.1-5.8. However, in certain skin conditions, like AD, an elevation of this pH has
<table><caption>Table 8: Experimental data from Dasgupta et. al [12] showing the influence of pH on growth rates of *S. epidermidis* and *S. aureus*</caption><thead><tr><th rowspan="2">pH</th><th colspan="2">Growth rate (ΔOD/hour)</th></tr><tr><th><i>S. aureus</i></th><th><i>S. epidermidis</i></th></tr></thead><tbody><tr><td>5</td><td>0.03</td><td>0.05</td></tr><tr><td>5.5</td><td>0.04</td><td>0.07</td></tr><tr><td>6</td><td>0.09</td><td>0.08</td></tr><tr><td>6.5</td><td>0.13</td><td>0.09</td></tr><tr><td>7</td><td>0.14</td><td>0.10</td></tr></tbody></table>

been observed. Dasgupta *et al.* studied *in vitro* the influence of pH on the growth rates of *S. aureus* and *S. epidermidis* [12]. Their experimental results show that, when the pH is increased from 5 to 6.5, the growth rate of *S. epidermidis* is multiplied by 1.8, whereas the one of *S. aureus* is multiplied by more than 4 (Table 8).

Their data can be used to select values for the growth rates $r_{sc}$ and $r_{sp}$ in our model, corresponding to healthy skin with a skin surface pH of 5 and compromised skin with a pH of 6.5. Because the experiments from Dasgupta *et al.* were performed *in vitro* and the bacterial population sizes measured with optical density (OD) instead of CFU, the growth rates cannot be directly translated into $r_{sc}$ and $r_{sp}$. Indeed, the experiments done by Dasgupta *et al.* were performed *in vitro* whereas the model is designed to represent the microbial dynamics in vivo. Even though in vitro experiments provide insightful results on the biological mechanisms studied, they also lack important factors such as interactions with the skin providing both nutrients and antimicrobial peptides, and skin renewal. In in vitro experiments the bacteria usually have access to an optimum amount of nutrients that can be renewed during the experiment or not. In vivo, the nutrients are provided by the skin. Depending on the skin state (dry, moisturized, with eczema,...) this quantity might not be optimal for bacterial growth. However, nutrients are continuously renewed due to new surface cells arising.

We use $r_{sc} = 0.5$ as the reference value for the commensal growth rate at pH 5, following on from previous simulation (Fig. 3). Maintaining the ratio between the two population growth rates at pH 5 and the multiplying factors following the pH elevation from Dasgupta *et al.* experimental data, we can define two sets of values for $r_{sc}$ and $r_{sp}$:

$$ \text{skin surface pH of 5} \Rightarrow r_{sc} = 0.5, r_{sp} = 0.3 $$

$$ \text{skin surface pH of 6.5} \Rightarrow r_{sc} = 0.9, r_{sp} = 1.3 $$
<figure><img src="image_5.png" /><figcaption>Figure 4: Landscape of satisfaction degree of the temporal formula corresponding to healthy skin with a skin surface pH of 5 ($r_{sc}$ = 0.5 and $r_{sp}$ = 0.3). The x and y axis represent variations of the initial quantities of $[S_p]$ and $[S_c]$, respectively. The color coding corresponds to the satisfaction degree of the temporal logic formula. Values used for the other parameters: $C_1 = 2.10^8$, $k_c = 0.01$, $[Amp_h] = 3$.</figcaption></figure>

Considering the healthy skin scenario with a skin surface pH of 5, the influence of the bacterial populations initial concentrations on the dominant species after 50 hours is evaluated using the temporal logic formula:

$$F(\text{Time} == 40 \land ([S_c] > u1 [S_p]) \land F(G([S_c] > u2 [S_p]))))$$

where $u1$ and $u2$ are free variables representing the abundance factors between both populations, evaluated at Time= 40 and at the last time point of the trace respectively (F stands for finally and G for globally at all future time points), i.e. at the time horizon of the experiments of 50 hours. It is worth remarking that with that temporal formula, we do not impose that an equilibrium is reached since we want to focus primarily on the dominance criterion in this analysis.

When given with an objective value, e.g. $u1 = 10$, the distance between that value and the validity domain of the formula, i.e. the set of values for $u1$ that satisfy the formula, provides a violation degree which is used to evaluate the satisfaction degree of the property.

Here, we evaluate how much the temporal formula $F(\text{Time} == 40 \land ([S_c] > u1 [S_p]) \land F(G([S_c] > u2 [S_p]))))$, $u1 \to 10$, $u2 \to 10$, is satisfied given variations of the initial concentrations of two populations (Fig. 4). The model predicts that, under the healthy skin condition, the commensal population will always dominate after 50 hours, except when introduced at a relatively low concentration ($< 2.10^4$) while the initial concentration of the pathogenic population is high ($> 5.10^5$).
<figure><img src="image_6.png" /><figcaption>Figure 5: Landscape of satisfaction degree of the temporal formula corresponding to compromised skin with a skin surface pH of 6.5 ($r_{sc}$ = 0.9 and $r_{sp}$ = 1.3). The x and y axis represent variations of the initial quantities of [$S_p$] and [$S_c$], respectively. The color coding corresponds to the satisfaction degree of the temporal logic formula. Values used for the other parameters: $C_1$ = 2.10$^8$, $k_c$ = 0.01, [$Amp_h$] = 3.</figcaption></figure>

The model predicts a higher vulnerability of the skin regarding invading pathogens with an elevated skin surface pH. When evaluating the same temporal formula with growth rates values corresponding to a skin surface pH of 6.5, we observe that even when the initial concentration of commensal is high (> 10$^7$), the pathogenic population is able to colonize the skin when introduced at a concentration as low as 3.10$^4$ (Fig. 5). Such predictions highlight the protective effect of the skin surface acidic pH against the invasion of pathogenic bacteria.

It was expected that an elevated pH would favor the pathogens. It is worth remarking however that the magnitude of the pathogens dominance compared to the small difference in growth rates ($r_{sp}$ 44% higher than $r_{sc}$) was surprising to us. Modeling more explicitly the skin surface pH would require a much more complex model since the acidic skin surface pH depends on the concentration in Natural Moisturizing Factor (NMF). NMF is constituted of small molecules resulting from the degradation of other molecules such as filaggrin deeper within the skin. The exact mechanisms involved in establishing and maintaining the acidic skin surface pH are still not completely understood. Short of modeling the pH impact on bacteria on a molecular scale, for which there is not much data available, changes in pH can only be modeled through a change in bacterial growth rates, death rates or carrying capacities ($K_{sp}$ and $K_{sc}$). Here, we chose to impact the growth rates based on the data that was available.

## 5.2. Reduced production of skin AMPs

As mentioned before, human keratinocytes constitutively produce AMPs as a defense against pathogens. In atopic dermatitis, the expression of AMPs
is dysregulated, leading to lower concentration levels of AMPs in the epidermis [39]. Similarly to the analysis done for skin surface pH, our model can be used to study how the skin microbiome reacts to modulation of the AMPs production by the skin cells. Two situations are considered: an impaired production of AMPs by the skin cells ($[Amp_h] = 0.5$) and a higher concentration with $[Amp_h] = 3$. Using the same methodology as in the case of skin surface pH, the temporal logic formula $F(\text{Time} == 40 \land ([S_c] > u1 [S_p]) \land F(G([S_c] > u2 [S_p]))), u1 \rightarrow 10, u2 \rightarrow 10$, is evaluated for variations of the initial concentrations of both populations for $[Amp_h] = 0.5$ and $[Amp_h] = 3$ (Fig. 6).

The model predicts a slightly protective effect of $Amp_h$ regarding the colonization of the skin by a pathogenic population, for low initial concentrations. However when both populations are introduced in high concentrations, the increase of $[Amp_h]$ appears to have the opposite effect of facilitating the colonization by the pathogenic population.

This mitigated effect is unexpected from a biological standpoint and might be due to the presence of $[Amp_h]$ in the constraint related to the degradation rate of $[Amp_b]$ (Equation (14)). That equation was derived from a steady-state assumption made on our model to relate it to some available experimental observations. It resulted in making $[Amp_b]$ degradation rate $d_a$ proportional to the concentration in $[Amp_h]$. There is no evidence in the literature that an increased concentration in human AMPs might enhance the degradation of bacterial AMPs. The predicted non-linear protective effect of human AMPs on the skin should thus be further investigated, either experimentally or by refinement of our model. In case of discrepancy with new experimental evidence, that model prediction could serve as a guide towards a modification of the model structure or a reconsideration of the previous mathematical steady-state assumption made on the experimental data.

## 6. Quasi-stability revealed by simulation on a long time scale

Interestingly, by extending the simulation time horizon to a longer time scale of 500 hours, one can observe a quasi-stability phenomenon, shown in Fig. 7. The equilibrium observed in Fig. 3 at the relevant time scale of 50 hours of the experiments, is thus not a mathematically stable state, but a quasi-stable state, that slowly evolves, with $\frac{d[S_c]}{dt} \neq 0$ and $\frac{d[S_p]}{dt} \neq 0$, towards a true stable state of the model reached after 300 hours in which the population density are reversed.

The $S_c$ population almost reaches its optimum capacity $K_{sc}$ after approximately 30 hours and stays relatively stable for around 100 hours more, that is over 4 days, which can reasonably be considered stable on the microbiological time scale. Meanwhile, the $S_p$ population is kept at a low concentration compared to $S_c$, even though it is continuously increasing and eventually leading to its overtake of $S_c$.

By varying the parameters values, it appears that this quasi-stability phenomenon emerges above a threshold value of 2.5 for $[Amp_h]$, that is for almost
<figure><img src="image_7.png" /><figcaption>Figure 6: Landscape of satisfaction degree of the healthy condition formula with a low concentration of human AMPs on the upper graph ([$Amp_h$] = 0.5) and a high concentration at the bottom ([$Amp_h$] = 3). The x and y axis represent variations of the initial quantities of [$S_p$] and [$S_c$] respectively. The color coding corresponds to the satisfaction degree of the temporal logic formula. Values used for the other parameters: $r_{sc} = r_{sp} = 0.5$, $C_1 = 2.10^8$, $k_c = 0.01$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Numerical simulation of the reduced ODE model on a longer time scale of 500 hours, with the same initial concentrations and parameter values as in Fig. 3, showing an inversion of the dominant bacterial population after 220 hours.</figcaption></figure>

half of its possible values (see Section 4).

Such phenomena of quasi-stability, also called meta-stability when the switches result from stochasticity [4] appear in different forms in the theory of deterministic dynamical systems [56, 47, 53, 52, 14, 37]. The observation of this phenomenon in our simple population dynamics model is however quite surprising and deserves an original mathematical analysis to understand its sources.

# 7. Mathematical Analysis of Quasi-stability versus Stability in Population Dynamics

## 7.1. Quasi-stability versus stability

Dynamical systems theory is widely applied in biology. Notions such as attractors and steady states are systematically used for the analysis of biological models. Stable steady states (point attractors) are zero velocity states attracting all neighboring orbits (state S in Figure 8 a)). Quasi-stability phenomena are particularly well-studied in the case of oscillatory systems, e.g. in models of brain activity [56]. They are also central in computational systems biology for simplifying biochemical reaction networks. For instance, they provide a model reduction method based on the identification of different regimes corresponding to different preponderant terms of the ODEs, for which simplified dynamics can be defined [51, 48, 29, 13], and chained within a hybrid automaton [47, 14]. The simplified dynamics takes place on attractive, normally hyperbolic invariant manifolds, sometimes defined as quasi-steady states like state Q in Figure 8 a. These states attract neighboring orbits along some fast, stable directions, and evolve very slowly along other directions [29, 13].

Quasi-stable states were discussed for models in ecology [37], but are certainly more general. A non-exhaustive list of nonlinear phenomena leading to
<figure><img src="image_9.png" /><figcaption>Figure 8: Quasi-stability versus stability illustrated for the 2D skin microbiome model obtained by eliminating the fast variable [$Amp_b$] and using the parameters of Fig. 3. a) A typical orbit leaves the neighborhood of the unstable state ($S_c = 0, S_p = 0$), reaches the quasi-stable state Q where it stays for a long time and ends with the stable state S. b) The corresponding time series with different values of $S_p(0)$ shows that for a long time $S_p$ is constant whereas $S_c$ increases exponentially but very slowly starting with very low values. The transition to the stable steady state occurs when $S_c = 2 \times 10^8$ and the time spent in the quasi-stable state is given by $\log(2 \times 10^8/S_c(0))/\alpha$ where $S_c(0)$ and $\alpha$ are the initial values and effective exponential growth rate of the pathogen, respectively.</figcaption></figure>
quasi-stability include: i) slow attractive, normally hyperbolic invariant manifolds, that exist for slow/fast systems [58, 29, 13]; ii) saddle points with slow unstable directions, where slowness can be due to the proximity of a saddle-node bifurcation [37]; and iii) ghost attractors, also due to the proximity of a saddle-node bifurcation [37].

## 7.2. Quasi-stability of the reduced parametric model

We first look for slow invariant manifolds by using tropical geometry scaling methods introduced elsewhere [41, 42, 54, 46]. By these methods we identify fast and slow variables and use quasi-steady state as zero order approximations to invariant manifolds [51, 48, 29, 13].

For the sake of simplicity we rename the model variables and parameters from (1). The analysis of the model is performed using formal parameters. However, we use the nominal parameter values from Table 4 and Figure 3 to compute parameter orders of magnitude that are considered fixed throughout this section. These choices lead to

$$ 
\begin{aligned}
\dot{x}_1 &= k_1 x_1 \left( 1 - \frac{x_1}{k_2} - \frac{k_3 x_2}{k_4 + x_2} \right), \\
\dot{x}_2 &= k_5 x_2 \left( k_9 - \frac{x_2}{k_6} - \frac{k_7 x_3}{k_8 + x_3} \right), \\
\dot{x}_3 &= k_{10} (x_1 - k_{11} x_3),
\end{aligned} 
\quad (20) $$

where $x_1 = [S_c], x_2 = [S_p], x_3 = [Amp_b], k_1 = r_{sc} = 0.5, k_2 = K_{sc} = 4 \times 10^8, k_3 = d_{sc}/r_{sc} = 9/10, k_4 = C_1 = 2 \times 10^8, k_5 = r_{sp} = 0.5, k_6 = K_{sp} = 3 \times 10^9, k_7 = d_{spb}/r_{sp} = 5/4, k_8 = C_{ab} = 0.16, k_9 = 1 - \frac{d_{sph}[Amp_h]}{r_{sp}(C_{ah} + [Amp_h])} = 5/11, k_{10} = K_c = 0.01, k_{11} = d_a/K_c = \frac{149}{256} \times 10^{10}$.

The first step of tropical scaling is to rescale parameters as powers of a small positive parameter $\epsilon$,

$$ k_i = \bar{k}_i \epsilon^{\gamma_i}, $$

where $\gamma_i = \text{round}(\log(k_i)/\log(\epsilon))$ [51, 48, 29, 13].

For instance, for $\epsilon = 1/11$, we have $\bar{k}_1 = 0.5, \gamma_1 = 0, \bar{k}_2 = 1.86, \gamma_2 = -8, \bar{k}_3 = 9/10, \gamma_3 = 0, \bar{k}_4 = 0.933, \gamma_4 = -8, \bar{k}_5 = 0.5, \gamma_5 = 0, \bar{k}_6 = 1.27, \gamma_6 = -9, \bar{k}_7 = 5/4, \gamma_7 = 0, \bar{k}_8 = 1.76, \gamma_8 = 1, \bar{k}_9 = 5/11, \gamma_9 = 0, \bar{k}_{10} = 1.21, \gamma_{10} = 2, \bar{k}_{11} = 2.46, \gamma_{11} = -9$. Here the parameters $k_1, k_5, k_4, k_{10}$ and $k_{11}$ were chosen as in Fig. 3 to illustrate quasi-stability.

By doing so, parameters are represented by their orders of magnitude $\gamma_i$. In the generic situation when the rescaled, zero order parameters $\bar{k}_i$ are independent, $\gamma_i$ determine alone the orders of magnitude of the system's timescales [51, 48, 29, 13]. However, when $\bar{k}_i$ are dependent (for instance when particular expressions of these parameters are small in absolute value) the orders of magnitude of the timescales may also depend on $\bar{k}_i$. We will see that our problem corresponds to the non-generic case.
The second step of the tropical scaling method is to introduce all of order variables $x_i = y_i \epsilon^{a_i}$ and to consider that the orders of magnitude $a_i$ of the variables satisfy tropical equilibration conditions.

Although the tropical scaling method has been most often presented in the framework of polynomial ODEs [51, 48, 29, 13], early tropicalization ideas use rational ODEs [41]. The idea is to decompose the r.h.s. of each ODE into a sum of rational terms that have constant signs. For each ODE, a tropical equilibration condition means that its r.h.s. contains at least two such rational terms of opposite signs and dominant (smallest) order of magnitude. The orders of magnitude of rational terms are computed using the tropical algebra (with addition replaced by minimum and multiplication by plus). Therefore, the tropical equilibration conditions read:

$$ 
\begin{aligned}
0 &= \min(-\min(-a_2 - 8, 0), a_1 + 8) \\
0 &= \min(a_2 + 9, -\min(1 - a_3, 0)) \\
0 &= a_1 - a_3 + 9
\end{aligned} 
\qquad (21) $$

Using the last equation we obtain the equivalent system

$$ 
\begin{aligned}
0 &= \min(-\min(-a_2 - 8, 0), a_1 + 8) \\
0 &= \min(a_2 + 9, -\min(-a_1 - 8, 0)) \\
a_3 &= a_1 + 9
\end{aligned} 
\qquad (22) $$

The system (22) has two one dimensional branches of solutions

$$ \textbf{I.} \quad a_1 = -8, a_2 > -8, a_3 = 1. $$

$$ \textbf{II.} \quad a_2 = -9, a_1 > -8, a_3 = a_1 + 9. $$

The first branch corresponds to large $x_1$ (commensal bacteria) and small $x_2$ (pathogens), whereas the second branch corresponds to small $x_1$ and large $x_2$.

For the initial conditions in this study the branch I is appropriate. The rescaled system is then

$$ 
\begin{aligned}
\dot{y}_1 &= \bar{k}_1 y_1 \left( 1 - \frac{y_1}{\bar{k}_2} - \frac{\bar{k}_3 y_2 \epsilon^{a_2+8}}{\bar{k}_4 + y_2 \epsilon^{a_2+8}} \right), \\
\dot{y}_2 &= \bar{k}_5 y_2 \left( \bar{k}_9 - \epsilon^{a_2+9} \frac{\bar{y}_2}{\bar{k}_6} - \frac{\bar{k}_7 y_3}{\bar{k}_8 + y_3} \right), \\
\dot{y}_3 &= \epsilon^{-7} \bar{k}_{10} (y_1 - \bar{k}_{11} y_3).
\end{aligned} 
\qquad (23) $$

From the rescaled system one can see that $y_3$ is fast and $y_1$ and $y_2$ are slow. Therefore, one can reduce the system of 3 ODEs to a system of 2 ODEs by using the quasi-steady approximation for the variable $y_3$, namely $y_3 = y_1/\bar{k}_{11}$. This reduction is robust as it is valid for low and high values of the sub-populations.
On the branch I, because $a_2 + 8 > 0$ the terms scaling like $\epsilon^{a_2+8}$ and $\epsilon^{a_2+9}$ can be neglected. The 2D reduced order model reads

$$ 
\begin{aligned}
\dot{y}_1 &= \bar{k}_1 y_1 \left(1 - \frac{y_1}{\bar{k}_2}\right), \\
\dot{y}_2 &= \bar{k}_5 y_2 \left(\bar{k}_9 - \frac{\bar{k}_7 y_1}{\bar{k}_8 \bar{k}_{11} + y_1}\right).
\end{aligned} 
\quad (24) $$

In the system (24) the first equation is decoupled from the second. In a short time, $y_1 \to \bar{k}_2$ and then $y_2$ follows the equation

$$ \dot{y}_2 = \bar{k}_5 y_2 \left( \bar{k}_9 - \frac{\bar{k}_7 \bar{k}_2}{\bar{k}_8 \bar{k}_{11} + \bar{k}_2} \right) \quad (25) $$

The condition for this to happen is that $y_1$ is faster than $y_2$ which reads

$$ 0 < \bar{k}_5(\bar{k}_9 - \frac{\bar{k}_7\bar{k}_2}{\bar{k}_8\bar{k}_{11} + \bar{k}_2}) \ll \bar{k}_1 \quad (26) $$

Using the numerical values of the parameters we find

$$ \bar{k}_5(\bar{k}_9 - \frac{\bar{k}_7\bar{k}_2}{\bar{k}_8\bar{k}_{11} + \bar{k}_2}) = 0.0790\bar{k}_1. $$

The variable $y_2$ is slow with a characteristic time

$$ \left[ \bar{k}_5 \left( \bar{k}_9 - \frac{\bar{k}_7 \bar{k}_2}{\bar{k}_8 \bar{k}_{11} + \bar{k}_2} \right) \right]^{-1} = (0.5 \times 0.0790)^{-1} = 25.32h, $$

meaning that $x_2$ (like $y_2$) needs 58.32h to change by a decade. That explains the observed quasi-stability.

The numerical value of the characteristic time was obtained from the nominal parameters used in Figure 3. For different values of the zero order rescaled parameters $\bar{k}_i$ one should verify our condition (26) in order to conclude about quasi-stability. Furthermore, if orders of magnitude of the parameters change, the tropical scaling and thus the entire analysis change. A full formal parametric analysis, discussing all possible tropical scalings depending on parameter orders of magnitude, is possible but left for a different publication.

### 7.3. Non-genericity of this form of quasi-stability

In contrast to slow-fast systems, where quasi-stability is generic (depends only on orders of magnitudes of parameters) here one has non-generic slowing down because of compensation between zero order parameters $\bar{k}_i$ (condition (26)).

The effective dynamics in the quasi-stable state consists in very slow exponential growth of the pathogen population $x_2$ with the effective rate

$$ \alpha = \bar{k}_5(\bar{k}_9 - \frac{\bar{k}_7\bar{k}_2}{\bar{k}_8\bar{k}_{11} + \bar{k}_2}), \quad (27) $$
while the commensal population $x_1$ is constant, see also Figure 8. When the pathogen population becomes too large the conditions needed for the existence of the branch I are no longer satisfied and all the above approximations break down. The system exits quasi-stability and evolves towards the stable steady state S (Figure 8). As a matter of fact, the time needed to exit Q corresponds to the time needed for $x_2$ to reach $k_4$. When this happens, the term $\frac{\bar{k}_3 y_2 \epsilon^{a_2+8}}{k_4+y_2 \epsilon^{a_2+8}}$ can no longer be neglected in the first equation of the system (24). This time depends on the initial value $x_2(0)$, scales like $[\log(k_4) - \log(x_2(0))]\alpha^{-1}$, and becomes infinite for $x_2(0) = 0$.

### 7.4. General conditions in population dynamics models

One would expect that this slowing down mechanism by compensation of zero order parameters occurs more generally in models of antagonistic species. In order to illustrate this, let us consider the following model describing competition of N species:

$$ \dot{x}_i = r_i x_i \left(1 - \frac{x_i}{K_i} - \sum_{j \neq i} \frac{\beta_{ij} x_j}{\theta_{ij} + x_j}\right), \quad 1 \le i \le N, \qquad (28) $$

where $x_i$ are the species concentrations.

We consider that the orders of magnitude of $r_i$, $\alpha_i$, $K_i$, $\beta_{ij}$ and $\theta_{ij}$ are all zero. More generally, for parameters with different orders of magnitude, we should first use tropical scaling to obtain a rough fast/slow decomposition of the ODE system. In the general case, (28) represents the slowest dynamics resulting from tropical rescaling and describing interactions among slow species, direct or mediated by fast species.

Let us further consider two categories of species defined by indices subsets $S_1$ and $S_2$, where $S_1 \cap S_2 = \emptyset$, $S_1 \cup S_2 = \{1, \dots, N\}$. $S_1$ includes the indices of saturated species, whose concentrations will tend to the capacities $K_i$ and $S_2$ includes the indices of species whose initial concentrations are close to zero, initially. For simplicity, we consider that species $S_1$ do not compete one another, i.e. $\beta_{ij} = 0$ for $i, j \in S_1$. Under these conditions, the dominating terms in the ODE model are

$$ \dot{x}_i = r_i x_i \left(1 - \frac{x_i}{K_i}\right), \quad \text{for } i \in S_1, \qquad (29) $$

$$ \dot{x}_i = r_i x_i \left(1 - \sum_{j \in S_1} \frac{\beta_{ij} x_j}{\theta_{ij} + x_j}\right), \quad \text{for } i \in S_2. \qquad (30) $$

The approximate Eq.(29) was obtained by neglecting the terms $\sum_{j \in S_2} \frac{\beta_{ij} x_j}{\theta_{ij} + x_j}$, which is valid as long as $x_j \ll \theta_{ij}$ for all $i \in S_1, j \in S_2$.

According to (29) the species $S_1$ stabilize to saturation levels $K_i$ and according to (30) the species $S_2$ increase exponentially with rates $r_i(1 - \sum_{j \in S_1} \frac{\beta_{ij} K_j}{\theta_{ij} + K_j})$. Therefore, the quasi-stability conditions read:

$$ 0 < r_i \left(1 - \sum_{j \in S_1} \frac{\beta_{ij} K_j}{\theta_{ij} + K_j}\right) \ll 1, \quad \text{for all } i \in S_2 \qquad (31) $$

$$ x_j(0) \ll \theta_{ij}, \quad \text{for all } i \in S_1, j \in S_2 \qquad (32) $$
The non-negativity and the much-less-than conditions in (31) ensure that concentrations of species from $S_2$ increase, and that the increase is slow, respectively. The much-less-than condition in (32) ensures that the approximate equation (29) is initially valid. If (31),(32) are satisfied, then the saturated species reach almost constant levels $x_i = K_i$, $i \in S_1$ after a short transient, and the non-saturated species increase very slowly. The lifetime of this quasi-stable state is the smallest time for some non-saturated species $i$ to reach the threshold level $\theta_{ij}$ in which case (29) is no longer true. This lifetime reads

$$ \min_{i \in S_2, j \in S_1} \frac{\log(\theta_{ij}) - \log(x_i(0))}{r_i(1 - \sum_{l \in S_1} \frac{\beta_{il}K_l}{\theta_{il} + K_l})} $$

## 8. Conclusion

In the perspective of identifying conditions which might favor or inhibit the emergence of pathogenic populations in the skin microbiome, we have developed a simple ODE model of skin microbiome with 3 variables and 13 parameters, reduced to 5 parameters by using, on the one hand, published data from the literature, and on the other hand, steady state reasoning for taking into account available data from biological experiments. Our bacterial population model is generic in the sense that we did not take into account the peculiarities of some specific bacterial populations, but on some general formulas of adversary population dynamics and influence factors. We showed through sensitivity analyses that our model predictions are reasonably robust with respect to parameter variations.

The presented model predicts that a) elevation of skin surface pH, creates favorable conditions for the emergence and colonization of the skin by the opportunistic pathogen populations to the detriment of commensals and b) the concentration of human AMPs is non-linearly affecting the balance between pathogens and commensals. The first prediction is in line for example with current knowledge of the effect of pH on staphylococcus aureus in AD. The model can be used to generalize this observation to other pathogens and can justify the use of optimal pH levels in topical formulations. The second prediction is novel and the specific non-linearity can be used to design strategies that optimize AMP levels on the skin surface using appropriately formulated topical treatments.

Perhaps surprisingly, we also showed that this simple model exhibits a quasi-stability phenomenon over a large range of biologically relevant parameter values, revealed by allowing the simulation to continue for times one order of magnitude longer than the reported time of the biological experiments. This observation questions the existence and importance of quasi-stability phenomena in real biological processes, whereas a natural assumption made in mathematical modeling, and model fitting to data, is that the experimental data are observed in states corresponding to real stable states of the mathematical model.

Using tropical scaling, we have shown that, in contrast to well-studied slow-fast systems, quasi-stability in our model does not depend on the orders of
magnitude of the parameters, but results from a non-generic slowing-down of
zero order parameters, for which our analysis provides analytic estimates of the
timescale dependence on the parameters. Finally, those conditions for quasi-
stability, obtained by tropical algebra analysis of our skin microbiome model,
have been shown to generalize to population dynamics models, with an arbitrary
number of species, for their potential applications in other domains.

Acknowledgments.

We are grateful to Mathieu Hemery, Aurélien Naldi and Sylvain Soliman for
interesting discussions on this work, and to the anonymous reviewers for their
useful comments to improve the presentation of our results.

**References**

[1] Alur, R., Feder, T., Henzinger, T.A.: The benefits of relaxing punctuality. J. ACM 43(1), 116–146 (Jan 1996). https://doi.org/10.1145/227595.227602

[2] Baudier, A., Fages, F., Soliman, S.: Graphical requirements for multistationarity in reaction networks and their verification in biomodels. Journal of Theoretical Biology 459, 79–89 (Dec 2018). https://doi.org/10.1016/j.jtbi.2018.09.024, https://hal.archives-ouvertes.fr/hal-01879735

[3] Belkaid, Y., Segre, J.A.: Dialogue between skin microbiota and immunity. Science 346(6212), 954–959 (Nov 2014). https://doi.org/10.1126/science.1260144, https://www.science.org/doi/10.1126/science.1260144

[4] Bovier, A., den Hollander, F.: Metastability. Methods of contemporary mathematical statistical physics 1970, 177–221 (1970)

[5] Byrd, A.L., Belkaid, Y., Segre, J.A.: The human skin microbiome. Nature Reviews Microbiology 16(3), 143–155 (Feb 2018). https://doi.org/10.1038/nrmicro.2017.157

[6] Calzone, L., Fages, F., Soliman, S.: BIOCHAM: An environment for modeling biological systems and formalizing experimental knowledge. Bioinformatics 22(14), 1805–1807 (2006). https://doi.org/10.1093/bioinformatics/btl172

[7] Campion, J.J., McNamara, P.J., Evans, M.E.: Pharmacodynamic Modeling of Ciprofloxacin Resistance in Staphylococcus aureus. Antimicrobial Agents and Chemotherapy 49(1), 209–219 (Jan 2005). https://doi.org/10.1128/AAC.49.1.209-219.2005, https://journals.asm.org/doi/10.1128/AAC.49.1.209-219.2005, publisher: American Society for Microbiology
[8] Chopra, R., Vakharia, P.P., Sacotte, R., Silverberg, J.I.: Efficacy of bleach baths in reducing severity of atopic dermatitis: A systematic review and meta-analysis. Annals of Allergy, Asthma & Immunology 119(5), 435-440 (Nov 2017). https://doi.org/10.1016/j.anai.2017.08.289, https://www.annallergy.org/article/S1081-1206(17)30958-4/fulltext, publisher: Elsevier

[9] Cogen, A.L., Yamasaki, K., Muto, J., Sanchez, K.M., Alexander, L.C., Tanios, J., Lai, Y., Kim, J.E., Nizet, V., Gallo, R.L.: Staphylococcus epidermidis antimicrobia δ-toxin (phenol-soluble modulin-γ) cooperates with host antimicrobial peptides to kill group a streptococcus. PLOS ONE 5(1), e8557 (Jan 2010). https://doi.org/10.1371/journal.pone.0008557, https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0008557, publisher: Public Library of Science

[10] Cogen, A.L., Yamasaki, K., Sanchez, K.M., Dorschner, R.A., Lai, Y., MacLeod, D.T., Torpey, J.W., Otto, M., Nizet, V., Kim, J.E., Gallo, R.L.: Selective antimicrobial action is provided by phenol-soluble modulins derived from staphylococcus epidermidis, a normal resident of the skin. Journal of Investigative Dermatology 130(1), 192-200 (Jan 2010). https://doi.org/10.1038/jid.2009.243, https://linkinghub.elsevier.com/retrieve/pii/S0022202X15345486

[11] Czock, D., Keller, F.: Mechanism-based pharmacokinetic-pharmacodynamic modeling of antimicrobial drug effects. Journal of Pharmacokinetics and Pharmacodynamics 34(6), 727-751 (Dec 2007). https://doi.org/10.1007/s10928-007-9069-x

[12] Dasgupta, A., Iyer, V., Raut, J., Qualls, A.: 16502 effect of pH on growth of skin commensals and pathogens. Journal of the American Academy of Dermatology 83(6), AB180 (Dec 2020). https://doi.org/10.1016/j.jaad.2020.06.080, https://www.jaad.org/article/S0190-9622(20)31908-3/abstract, publisher: Elsevier

[13] Desoeuvres, A., Iosif, A., Lüders, C., Radulescu, O., Rahkooy, H., Seiß, M., Sturm, T.: Reduction of chemical reaction networks with approximate conservation laws. arXiv preprint arXiv:2212.13474 (2022)

[14] Desoeuvres, A., Szmolyan, P., Radulescu, O.: Qualitative dynamics of chemical reaction networks: An investigation using partial tropical equilibrations. In: Petre, I., Păun, A. (eds.) Computational Methods in Systems Biology. pp. 61-85. Springer International Publishing, Cham (2022)

[15] Donzé, A., Maler, O.: Robust satisfaction of temporal logic over real-valued signals. In: FORMATS 2010. Lecture Notes in Computer Science, vol. 6246, pp. 92-106. Springer-Verlag (2010)
[16] Eberlein-König, B., Schäfer, T., Huss-Marp, J., Darsow, U., Möhrenschlager, M., Herbert, O., Abeck, D., Krämer, U., Behrendt, H., Ring, J.: Skin surface ph, stratum corneum hydration, trans-epidermal water loss and skin roughness related to atopic eczema and skin dryness in a population of primary school children. Acta Dermato-Venereologica 80(3), 188–191 (May 2000). https://doi.org/10.1080/000155500750042943

[17] Fages, F., Rizk, A.: On temporal logic constraint solving for the analysis of numerical data time series. Theoretical Computer Science 408(1), 55–65 (Nov 2008). https://doi.org/doi:10.1016/j.tcs.2008.07.004

[18] Fages, F., Soliman, S.: Abstract interpretation and types for systems biology. Theoretical Computer Science 403(1), 52–70 (2008). https://doi.org/10.1016/j.tcs.2008.04.024

[19] Greugny, E.T., Stamatas, G.N., Fages, F.: Stability versus meta-stability in a skin microbiome model. In: CMSB'22: Proceedings of the twentieth international conference on Computational Methods in Systems Biology. Lecture Notes in BioInformatics, vol. 13447. Springer-Verlag (Sep 2022), https://hal.inria.fr/hal-03698344

[20] Grice, E.A., Kong, H.H., Conlan, S., Deming, C.B., Davis, J., Young, A.C., Bouffard, G.G., Blakesley, R.W., Murray, P.R., Green, E.D., Turner, M.L., Segre, J.A.: Topographical and temporal diversity of the human skin microbiome. Science (New York, N.Y.) 324(5931), 1190–1192 (May 2009). https://doi.org/10.1126/science.1171700, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2805064/

[21] Hansen, N., Ostermeier, A.: Completely derandomized self-adaptation in evolution strategies. Evolutionary Computation 9(2), 159–195 (2001)

[22] Kennedy, E.A., Connolly, J., Hourihane, J.O., Fallon, P.G., McLean, W.H.I., Murray, D., Jo, J.H., Segre, J.A., Kong, H.H., Irvine, A.D.: Skin microbiome before development of atopic dermatitis: Early colonization with commensal staphylococci at 2 months is associated with a lower risk of atopic dermatitis at 1 year. J Allergy Clin. Immunol. 139(1), 7 (2017)

[23] Kobayashi, T., Glatz, M., Horiuchi, K., Kawasaki, H., Akiyama, H., Kaplan, D.H., Kong, H.H., Amagai, M., Naga, K.: Dysbiosis and staphylococcus aureus colonization drives inflammation in atopic dermatitis. Immunity 42(4), 756–766 (Apr 2015). https://doi.org/10.1016/j.immuni.2015.03.014, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4407815/

[24] Koh, L.F., Ong, R.Y., Common, J.E.: Skin microbiome of atopic dermatitis. Allergology International p. S1323893021001404 (Nov 2021). https://doi.org/10.1016/j.alit.2021.11.001, https://linkinghub.elsevier.com/retrieve/pii/S1323893021001404
[25] Kohda, K., Li, X., Soga, N., Nagura, R., Duerna, T., Nakajima, S., Nakagawa, I., Ito, M., Ikeuchi, A.: An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology. Frontiers in Cellular and Infection Microbiology 11, 712360 (Sep 2021). https://doi.org/10.3389/fcimb.2021.712360, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8481888/

[26] Kong, H.H.: Skin microbiome: genomics-based insights into the diversity and role of skin microbes. Trends in Molecular Medicine 17(6), 320-328 (Jun 2011). https://doi.org/10.1016/j.molmed.2011.01.013, https://www.sciencedirect.com/science/article/pii/S1471491411000232

[27] Kong, H.H., Oh, J., Deming, C., Conlan, S., Grice, E.A., Beatson, M.A., Nomicos, E., Polley, E.C., Komarow, H.D., NISC Comparative Sequence Program, Murray, P.R., Turner, M.L., Segre, J.A.: Temporal shifts in the skin microbiome associated with disease flares and treatment in children with atopic dermatitis. Genome Research 22(5), 850-859 (May 2012). https://doi.org/10.1101/gr.131029.111, http://genome.cshlp.org/lookup/doi/10.1101/gr.131029.111

[28] Korting, H.C., Hübner, K., Greiner, K., Hamm, G., Braun-Falco, O.: Differences in the skin surface ph and bacterial microflora due to the long-term application of synthetic detergent preparations of ph 5.5 and ph 7.0. results of a crossover trial in healthy volunteers. Acta Dermato-Venereologica 70(5), 429-431 (1990)

[29] Kruff, N., Lüders, C., Radulescu, O., Sturm, T., Walcher, S.: Algorithmic reduction of biological networks with multiple time scales. Mathematics in Computer Science 15(3), 499-534 (2021)

[30] Lai, Y., Cogen, A.L., Radek, K.A., Park, H.J., MacLeod, D.T., Leichtle, A., Ryan, A.F., Di Nardo, A., Gallo, R.L.: Activation of tlr2 by a small molecule produced by staphylococcus epidermidis increases antimicrobial defense against bacterial skin infections. The Journal of investigative dermatology 130(9), 2211-2221 (Sep 2010). https://doi.org/10.1038/jid.2010.123, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2922455/

[31] Lai, Y., Di Nardo, A., Nakatsuji, T., Leichtle, A., Yang, Y., Cogen, A.L., Wu, Z.R., Hooper, L.V., Schmidt, R.R., von Aulock, S., Radek, K.A., Huang, C.M., Ryan, A.F., Gallo, R.L.: Commensal bacteria regulate toll-like receptor 3-dependent inflammation after skin injury. Nature Medicine 15(12), 1377-1382 (Dec 2009). https://doi.org/10.1038/nm.2062

[32] Leyden, J.J., McGinley, K.J., Mills, O.H., Kligman, A.M.: Propionibacterium levels in patients with and without acne
vulgaris. Journal of Investigative Dermatology **65**(4), 382-384 (Oct 1975). https://doi.org/10.1111/1523-1747.ep12607634, https://www.sciencedirect.com/science/article/pii/S0022202X1544610X

[33] Maler, O., Nickovic, D.: Monitoring Temporal Properties of Continuous Signals, pp. 152-166. Springer Berlin Heidelberg (2004). https://doi.org/10.1007/978-3-540-30206-3_12

[34] Malthus, T.R.T.R.: An essay on the principle of population, as it affects the future improvement of society. With remarks on the speculations of Mr. Godwin, M. Condorcet and other writers. London, J. Johnson (1798), http://archive.org/details/essayonprinciples00malt

[35] Meredith, H.R., Lopatkin, A.J., Anderson, D.J., You, L.: Bacterial temporal dynamics enable optimal design of antibiotic treatment. PLOS Computational Biology **11**(4), e1004201 (Apr 2015). https://doi.org/10.1371/journal.pcbi.1004201, https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1004201, publisher: Public Library of Science

[36] Miyano, T., Irvine, A.D., Tanaka, R.J.: Model-based meta-analysis to optimize staphylococcus aureus-targeted therapies for atopic dermatitis. JID Innov (2002)

[37] Morozov, A., Abbott, K., Cuddington, K., Francis, T., Gellner, G., Hastings, A., Lai, Y.C., Petrovskii, S., Scranton, K., Zeeman, M.L.: Long transients in ecology: Theory and applications. Physics of Life Reviews **32**, 1-40 (2020)

[38] Nakaoka, S., Kuwahara, S., Lee, C.H., Jeon, H., Lee, J., Takeuchi, Y., Kim, Y.: Chronic inflammation in the epidermis: A mathematical model. Applied Sciences (2016)

[39] Nakatsuji, T. Chen, T.H. Narala, S., Chun, K.A., Two. A.M., Yun, T., Shafiq, F., Kotol, P.F., Bouslimani, A., Melnik, A.V., Latif, H., Kim, J.N., Lockhart, A., Artis, K., David, G., Taylor, P., Streib, J., Dorrestein, P.C., Grier, A., Gill, S.R., Zengler, K., Hata, T.R., Leung, D.Y.M., Gallo, R.L.: Antimicrobials from human skin commensal bacteria protect against staphylococcus aureus and are deficient in atopic dermatitis. Science translational medicine **9**(378), eaah4680 (Feb 2017). https://doi.org/10.1126/scitranslmed.aah4680, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5600545/

[40] Nakatsuji, T., Hata, T.R., Tong, Y., Cheng, J.Y., Shafiq, F., Butcher, A.M., Salem, S.S., Brinton, S.L., Rudman Spergel, A.K., Johnson, K., Jepson, B., Calatroni, A., David, G., Ramirez-Gama, M., Taylor, P., Leung, D.Y.M., Gallo, R.L.: Development of a human skin commensal microbe for bacteriotherapy of atopic dermatitis and use in a phase 1 randomized clinical trial. Nature Medicine
27(4), 700-709 (Apr 2021). https://doi.org/10.1038/s41591-021-01256-2
https://www.nature.com/articles/s41591-021-01256-2

[41] Noel, V., Grigoriev, D., Vakulenko, S., Radulescu, O.: Tropical geometries and dynamics of biochemical networks application to hybrid cell cycle models. In: Feret, J., Levchenko, A. (eds.) Proc. SASB 2011, ENTCS, vol. 284, pp. 75-91. Elsevier (2012). https://doi.org/10.1016/j.entcs.2012.05.016

[42] Noel, V., Grigoriev, D., Vakulenko, S., Radulescu, O.: Tropicalization and tropical equilibration of chemical reactions. In: Litvinov, G.L., Sergeev, S.N. (eds.) Tropical and Idempotent Mathematics and Applications, Contemporary Mathematics, vol. 616, pp. 261-277. AMS (2014). https://doi.org/10.1090/conm/616/12316

[43] Ong, P.Y., Ohtake, T., Brandt, C., Strickland, I., Boguniewicz, M., Ganz, T., Gallo, R.L., Leung, D.Y.M.: Endogenous antimicrobial peptides and skin infections in atopic dermatitis. The New England Journal of Medicine 347(15), 1151-1160 (Oct 2002). https://doi.org/10.1056/NEJMoa021481

[44] Pazgier, M., Hoover, D.M., Yang, D., Lu, W., Lubkowski, J.: Human β-defensins. Cellular and Molecular Life Sciences CMLS 63(11), 1294-1313 (Jun 2006). https://doi.org/10.1007/s00018-005-5540-2, http://link.springer.com/10.1007/s00018-005-5540-2

[45] Proksch, E.: pH in nature, humans and skin. The Journal of Dermatology 45(9), 1044-1052 (2018). https://doi.org/10.1111/1346-8138.14489, https://onlinelibrary.wiley.com/doi/abs/10.1111/1346-8138.14489

[46] Radulescu, O.: Tropical geometry of biological systems (invited talk). In: Computer Algebra in Scientific Computing, LNCS 12291. pp. 1-13. Springer (2020)

[47] Radulescu, O., Samal, S.S., Naldi, A., Grigoriev, D., Weber, A.: Symbolic dynamics of biochemical pathways as finite states machines. In: Roux, O.F., Bourdon, J. (eds.) Computational Methods in Systems Biology - 13th International Conference, CMSB 2015, Nantes, France, September 16-18, 2015, Proceedings. Lecture Notes in Computer Science, vol. 9308, pp. 104-120. Springer (2015). https://doi.org/10.1007/978-3-319-23401-4_10

[48] Radulescu, O., Vakulenko, S., Grigoriev, D.: Model reduction of biochemical reactions networks by tropical analysis methods. Math. Model. Nat. Pheno. 10(3), 124-138 (Jun 2015). https://doi.org/10.1051/mmnp/201510310

[49] Rizk, A., Batt, G., Fages, F., Soliman, S.: A general computational method for robustness analysis with applications to synthetic gene networks. Bioinformatics 12(25), il69-il78 (Jun 2009). https://doi.org/10.1093/bioinformatics/btp200
[50] Rizk, A., Batt, G., Fages, F., Soliman, S.: Continuous valuations of temporal logic specifications with applications to parameter optimization and robustness measures. Theoretical Computer Science 412(26), 2827-2839 (2011). https://doi.org/10.1016/j.tcs.2010.05.008

[51] Samal, S.S., Grigoriev, D., Fröhlich, H., Weber, A., Radulescu, O.: A geometric method for model reduction of biochemical networks with polynomial rate functions. Bulletin of Mathematical Biology 77(12), 2180-2211 (2015)

[52] Samal, S.S., Krishnan, J., Lueders, C., Weber, A., Radulescu, O., et al.: Metastable regimes and tipping points of biochemical networks with potential applications in precision medicine. bioRxiv preprint (2018)

[53] Samal, S.S., Naldi, A., Grigoriev, D., Weber, A., Théret, N., Radulescu, O.: Geometric analysis of pathways dynamics: Application to versatility of tgf-β receptors. Biosystems 149, 3-14 (2016)

[54] Soliman, S., Fages, F., Radulescu, O.: A constraint solving approach to model reduction by tropical equilibration. Algorithms for Molecular Biology 9(1), 24 (2014)

[55] Spalding, C., Keen, E., Smith, D.J., Krachler, A.M., Jabbari, S.: Mathematical modelling of the antibiotic-induced morphological transition of pseudomonas aeruginosa. PLoS computational biology (2018)

[56] Tognoli, E., Kelso, J.A.S.: The metastable brain. Neuron 81(1), 35-48 (01 2014). https://doi.org/10.1016/j.neuron.2013.12.022, https://pubmed.ncbi.nlm.nih.gov/24411730

[57] Tuomanen, E., Cozens, R., Tosch, W., Zak, O., Tomasz, A.: The rate of killing of Escherichia coli by beta-lactam antibiotics is strictly proportional to the rate of bacterial growth. Journal of General Microbiology 132(5), 1297-1304 (May 1986). https://doi.org/10.1099/00221287-132-5-1297

[58] Wiggins, S.: Normally hyperbolic invariant manifolds in dynamical systems, vol. 105. Springer Science & Business Media (1994)

[59] Zwietering, M.H., Jongenburger, I., Rombouts, F.M., van 't Riet, K.: Modeling of the bacterial growth curve. Applied and Environmental Microbiology 56(6), 1875-1881 (Jun 1990), https://www.ncbi.nlm.nih.gov/pmc/articles/PMC184525/