Eléa Thibault Greugny

François Fages

Ovidiu Radulescu

Peter Szmolyan

Georgios N Stamatas

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

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

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 differential equations, with 2 types of bacteria populations (skin commensals and opportunistic pathogens) and including the production of antimicrobial peptides 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 colonization 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.

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) [START_REF] Pazgier | Human β-defensins[END_REF]. The physiologically acidic skin surface pH also contributes to control the growth of bacterial populations [START_REF] Proksch | pH in nature, humans and skin[END_REF][START_REF] Korting | 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[END_REF]. 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 [START_REF] Lai | Commensal bacteria regulate tolllike receptor 3-dependent inflammation after skin injury[END_REF][START_REF] Cogen | Staphylococcus epidermidis antimicrobia δ-toxin (phenol-soluble modulin-γ) cooperates with host antimicrobial peptides to kill group a streptococcus[END_REF][START_REF] Lai | Activation of tlr2 by a small molecule produced by staphylococcus epidermidis increases antimicrobial defense against bacterial skin infections[END_REF][START_REF] Kong | Skin microbiome: genomics-based insights into the diversity and role of skin microbes[END_REF][START_REF] Belkaid | Dialogue between skin microbiota and immunity[END_REF][START_REF] Byrd | The human skin microbiome[END_REF].

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) [START_REF] Leyden | Propionibacterium levels in patients with and without acne vulgaris[END_REF][START_REF] Kong | Temporal shifts in the skin microbiome associated with disease flares and treatment in children with atopic dermatitis[END_REF]. In the case of AD, the patient skin is often colonized by Staphylococcus aureus (S. aureus), especially on the lesions [START_REF] Kong | Temporal shifts in the skin microbiome associated with disease flares and treatment in children with atopic dermatitis[END_REF], 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 [START_REF] Chopra | Efficacy of bleach baths in reducing severity of atopic dermatitis: A systematic review and meta-analysis[END_REF]. On the other hand, treatments involving introduction of commensal species, like Staphylococcus hominis [START_REF] Nakatsuji | Development of a human skin commensal microbe for bacteriotherapy of atopic dermatitis and use in a phase 1 randomized clinical trial[END_REF] 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 that 200 years [START_REF] Malthus | An essay on the principle of population, as it affects the future improvement of society[END_REF]. Here 1 , we introduce a model based on ordinary differential equations (ODEs), describing the interactions of a population of 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 [START_REF] Nakaoka | Chronic inflammation in the epidermis: A mathematical model[END_REF][START_REF] Miyano | Model-based meta-analysis to optimize staphylococcus aureus-targeted therapies for atopic dermatitis[END_REF]. 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. [START_REF] Nakatsuji | Antimicrobials from human skin commensal bacteria protect against staphylococcus aureus and are deficient in atopic dermatitis[END_REF] and Kohda et al. [START_REF] Kohda | An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology[END_REF], 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 [START_REF] Calzone | BIOCHAM: An environment for modeling biological systems and formalizing experimental knowledge[END_REF], using a specification of the properties of the different bacterial population equilibria formalized in quantitative temporal logic FO-LTL(R lin ) [START_REF] Rizk | Continuous valuations of temporal logic specifications with applications to parameter optimization and robustness measures[END_REF].

On the time scale of the experiments of Kohda et al. [START_REF] Kohda | An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology[END_REF], 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 [START_REF] Bovier | Metastability. Methods of contemporary mathematical statistical physics[END_REF][START_REF] Tognoli | The metastable brain[END_REF][START_REF] Radulescu | Symbolic dynamics of biochemical pathways as finite states machines[END_REF][START_REF] Samal | Geometric analysis of pathways dynamics: Application to versatility of tgf-β receptors[END_REF][START_REF] Samal | Metastable regimes and tipping points of biochemical networks with potential applications in precision medicine[END_REF][START_REF] Desoeuvres | Qualitative dynamics of chemical reaction networks: An investigation using partial tropical equilibrations[END_REF] or quasi-stability [START_REF] Morozov | Long transients in ecology: Theory and applications[END_REF], 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 phenomena on relatively long timescales.

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 [START_REF] Zwietering | Modeling of the bacterial growth curve[END_REF], 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 [START_REF] Campion | Pharmacodynamic Modeling of Ciprofloxacin Resistance in Staphylococcus aureus[END_REF][START_REF] Meredith | Bacterial temporal dynamics enable optimal design of antibiotic treatment[END_REF][START_REF] Spalding | Mathematical modelling of the antibiotic-induced morphological transition of pseudomonas aeruginosa[END_REF]. 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 [START_REF] Cogen | Selective antimicrobial action is provide by phenolsoluble modulins derived from staphylococcus epidermidis, a normal resident of the skin[END_REF][START_REF] Nakatsuji | Antimicrobials from human skin commensal bacteria protect against staphylococcus aureus and are deficient in atopic dermatitis[END_REF]. For these reasons, we chose to also introduce in our model AMPs of bacterial origin, Amp b , acting similarly to Amp h on the pathogenic population S p . Amp b is produced at rate k c by S c , and degraded at rate d a . Furthermore, we include a defense mechanism of S p against S c with a direct killing effect. This killing effect is also modelled with a Hill function or order 1, where d sc is the maximum killing rate of S c by S p and C 1 the concentration of S p 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:

                         d[S c ] dt = r sc 1 - [S c ] K sc - d sc [S p ] C 1 + [S p ] [S c ] d[S p ] dt = r sp 1 - [S p ] K sp - d spb [Amp b ] C ab + [Amp b ] - d sph [Amp h ] C ah + [Amp h ] [S p ] d[Amp b ] dt = k c [S c ] -d a [Amp b ] (1)
Table 1 recapitulates the 3 variables and 13 parameters of this model with their units.

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 [START_REF] Grice | Topographical and temporal diversity of the human skin microbiome[END_REF]. 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 [START_REF] Nakatsuji | Antimicrobials from human skin commensal bacteria protect against staphylococcus aureus and are deficient in atopic dermatitis[END_REF][START_REF] Kohda | An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology[END_REF]. 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 [START_REF] Kohda | An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology[END_REF], 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 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 [START_REF] Nakatsuji | Antimicrobials from human skin commensal bacteria protect against staphylococcus aureus and are deficient in atopic dermatitis[END_REF] 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.

Parameter values inferred from mono-culture experiment data

We consider first the monocultures experiments from Kohda et al. [START_REF] Kohda | An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology[END_REF], representing the simplest experimental conditions. S. epidermis is a representative of the commensal population S c , and S. aureus of the pathogenic one, S p . Since the two species are not interacting, the set of equations simplifies to:

             d[S c ] dt = r sc 1 - [S c ] K sc [S c ] d[S p ] dt = r sp 1 - [S p ] K sp [S p ] (2)
At steady-state, the population concentrations are either zero, or equal to their optimum capacities (K sc or K sp ) 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 CF U.ASU -1 (3)
K sp = 3.10 9 CF U.ASU -1 (4)

Parameter relations inferred from experimental data on AMP

The experimental conditions of Nakatsuji et al. [START_REF] Nakatsuji | Antimicrobials from human skin commensal bacteria protect against staphylococcus aureus and are deficient in atopic dermatitis[END_REF] 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:

d[S p ] dt = r sp 1 - [S p ] K sp - d spb [Amp b ] C ab + [Amp b ] - d sph [Amp h ] C ah + [Amp h ] [S p ] (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 or [S p ] * = K sp 1 - d spb [Amp b ] r sp (C ab + [Amp b ]) - d sph [Amp h ] r sp (C ah + [Amp h ]) (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 steadystate and substitute the experimental values measured [Amp h ] = 4 µM ; [S p ] * = 10 9 CFU, and [Amp h ] = 8 µM ; [S p ] * = 6.10 5 CFU, together with the values of K sc and K sp (from (3) and ( 4)) in [START_REF] Calzone | BIOCHAM: An environment for modeling biological systems and formalizing experimental knowledge[END_REF], to obtain the following equations:

             d sph r sp = 4 + C ah 6
d sph r sp = (10 4 -2)(C ah + 8) 8.10 4 (7)
which reduce to C ah = 8 and

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 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 [START_REF] Tuomanen | The rate of killing of Escherichia coli by beta-lactam antibiotics is strictly proportional to the rate of bacterial growth[END_REF].

To be consistent with the ranges of Sh-lantibiotics concentrations described in Nakatsuji et al. [START_REF] Nakatsuji | Antimicrobials from human skin commensal bacteria protect against staphylococcus aureus and are deficient in atopic dermatitis[END_REF], [Amp b ] should take positive values below 10. Given that

[Amp b ] * = 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:

k c d a ≤ 1 4.10 7 (8)

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. [START_REF] Kohda | An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology[END_REF]. At steady-state, the system (1) gives:

[S c ] * = 0 or [S c ] * = K sc 1 - d sc [S p ] * r sc (C 1 + [S p ] * ) (9) [S p ] * = 0 or [S p ] * = K sp 1 - d spb [Amp b ] r sp (C ab + [Amp b ]) - d sph [Amp h ] r sp (C ah + [Amp h ]) (10)
[Amp b ] * = k c [S c ] * d a (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 [START_REF] Cogen | Selective antimicrobial action is provide by phenolsoluble modulins derived from staphylococcus epidermidis, a normal resident of the skin[END_REF] to get the following parameter relation:

d sc r sc = 3 4.10 9 C 1 + 3 4 (12) 2 3 r sp = d sph [Amp h ] C ah + [Amp h ] + 10 8 d spb k c d a C ab + 10 8 k c (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 56 + 31[Amp h ] 2.56 (4 -[Amp h ]) with [Amp h ] < 4 ( 14
)

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 [START_REF] Czock | Mechanism-based pharmacokineticpharmacodynamic modeling of antimicrobial drug effects[END_REF][START_REF] Campion | Pharmacodynamic Modeling of Ciprofloxacin Resistance in Staphylococcus aureus[END_REF];

• 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 CF U.ASU -1 ) and is thus bounded by the optimum concentration of S p , i.e. K sp = 3.10 9 CF U.ASU -1 , as calculated in section 3.1 from [START_REF] Kohda | An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology[END_REF];

• k c , the production rate of [Amp b ] chosen to take values between 0 and 0.1 AU.h -1 .CF U -1 , and shown to have a limited impact on the steadystate 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:

                           d[S c ] dt = r sc [S c ] 1 - [S c ] 4.10 8 - 3 4
(10

-9 C 1 + 1)[S p ] C 1 + [S p ] d[S p ] dt = r sp [S p ] 1 - [S p ] 3.10 9 - 5[Amp b ] 0.64 + 4[Amp b ] - 2[Amp h ] 8 + [Amp h ] d[Amp b ] dt = k c [S c ] -10 8 [Amp b ] 56 + 31[Amp h ] 2.56(4 -[Amp h ]) (15)
A continuum of different parameter values exists however to fit either the experimental data of the pathogenic microbiome reported in [START_REF] Kohda | An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology[END_REF], 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.

Formal specification of dynamical behaviors in quantitative temporal logic

Quantitative temporal logics, such as Metric Temporal Logic MTL [START_REF] Alur | The benefits of relaxing punctuality[END_REF], Signal Temporal Logic [START_REF] Maler | Monitoring Temporal Properties of Continuous Signals[END_REF], or First-Order Linear Time Logic with real-valued linear constraints FO-LTL(R lin ) [START_REF] Rizk | Continuous valuations of temporal logic specifications with applications to parameter optimization and robustness measures[END_REF][START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF], provide very expressive formal languages to specify behaviors of continuous time dynamical systems [START_REF] Donzé | Robust satisfaction of temporal logic over real-valued signals[END_REF].

In particular, the definition of a continuous degree of satisfaction in the interval [0, 1] for FO-LTL(R lin ) formulae [START_REF] Rizk | A general computational method for robustness analysis with applications to synthetic gene networks[END_REF] 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 robustness 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 phenotype reported in Kohda et al. experiments [START_REF] Kohda | An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology[END_REF], and to a healthy phenotype. We use the implementation of temporal logic FO-LTL(R lin ) in Biocham [START_REF] Rizk | Continuous valuations of temporal logic specifications with applications to parameter optimization and robustness measures[END_REF] 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 9 and [Sc] = 10 8 [START_REF] Kohda | An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology[END_REF]. That first property is a simple curve fitting problem with two time points that can be formalized in FO-LTL(R lin ) by the following temporal logic formula with free variables x1, y1, x2, y2 for the population sizes 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:

F(T ime ∼ 40 ∧ Sc = x1 ∧ Sp = y1 ∧ F(G(Sc = x2 ∧ Sp = y2))) (16
∃c ∃p F(T ime ∼ 40∧S c > k * S p ∧S c = c∧S p = p∧F(G(r * S p = p∧s * S c = c))) (18) with objective values k = 10, r = 1, s = 1. ( 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 [START_REF] Fages | Abstract interpretation and types for systems biology[END_REF], 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 [START_REF] Rizk | Continuous valuations of temporal logic specifications with applications to parameter optimization and robustness measures[END_REF]. In Biocham, the covariance matrix adaptive evolution strategy (CMA-ES) [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] 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 (µ, σ 2 ) around their nominal values µ, with standard deviation σ = COV.µ given by a coefficient of variation COV .

Model parameter fitting to pathogenic experiments and sensitivity analyses

In order to reproduce what was observed by Kohda et. al [START_REF] Kohda | An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology[END_REF], 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 CF U/mL and 10 3 CF U/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. 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 [START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF] 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 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 6 , 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 [START_REF] Kohda | An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology[END_REF].

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 .

Parameter

Nominal value COV Sensitivity index 2).

r

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 [START_REF] Kong | Temporal shifts in the skin microbiome associated with disease flares and treatment in children with atopic dermatitis[END_REF]. 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 [START_REF] Kohda | An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology[END_REF] 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 [START_REF] Baudier | Graphical requirements for multistationarity in reaction networks and their verification in biomodels[END_REF].

In order to get closer to physiological context however, we prefer to retain 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.

Conditions favoring the pathogenic population

Whether the dysbiosis observed in AD is the cause or the result of the disease is unclear [START_REF] Kobayashi | Dysbiosis and staphylococcus aureus colonization drives inflammation in atopic dermatitis[END_REF][START_REF] Koh | Skin microbiome of atopic dermatitis[END_REF]. 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 [START_REF] Kennedy | 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[END_REF]. 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 [START_REF] Eberlein-König | 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[END_REF] and reduced production of AMPs [START_REF] Ong | Endogenous antimicrobial peptides and skin infections in atopic dermatitis[END_REF]), 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.

Skin surface pH elevation

According to Proksch [START_REF] Proksch | pH in nature, humans and skin[END_REF], 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 [START_REF] Dasgupta | 16502 effect of ph on growth of skin commensals and pathogens[END_REF]. 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 : skin surface pH of 5 ⇒ r sc = 0.5, r sp = 0.3 skin surface pH of 6.5 ⇒ r sc = 0.9, r sp = 1.3 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 evaluate using the temporal logic formula:

F (Time == 40 ∧ ([S c ] > u1 [S p ]) ∧ F (G([S c ] > u2 [S p ])))
where u1 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 (Time == 40 ∧ ([S c ] > u1 [S p ]) ∧ F (G([S c ] > u2 [S p ]))
), u1 → 10, u2 → 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 ). 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.

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 [START_REF] Nakatsuji | Antimicrobials from human skin commensal bacteria protect against staphylococcus aureus and are deficient in atopic dermatitis[END_REF]. 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 (Time == 40 ∧ ([S c ] > u1 [S p ]) ∧ F (G([S c ] > u2 [S p ]))
), u1 → 10, u2 → 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 steadystate 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.

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

d[S c ] dt ̸ = 0 and d[S p ] dt ̸ = 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 half of its possible values (see Section 4).

Such phenomena of quasi-stability, also called meta-stability when the switches result from stochasticity [START_REF] Bovier | Metastability. Methods of contemporary mathematical statistical physics[END_REF] appear in different forms in the theory of deterministic dynamical systems [START_REF] Tognoli | The metastable brain[END_REF][START_REF] Radulescu | Symbolic dynamics of biochemical pathways as finite states machines[END_REF][START_REF] Samal | Geometric analysis of pathways dynamics: Application to versatility of tgf-β receptors[END_REF][START_REF] Samal | Metastable regimes and tipping points of biochemical networks with potential applications in precision medicine[END_REF][START_REF] Desoeuvres | Qualitative dynamics of chemical reaction networks: An investigation using partial tropical equilibrations[END_REF][START_REF] Morozov | Long transients in ecology: Theory and applications[END_REF]. 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.

Mathematical Analysis of Quasi-stability versus Stability in Population Dynamics

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 [START_REF] Tognoli | The metastable brain[END_REF]. 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 [START_REF] Samal | A geometric method for model reduction of biochemical networks with polynomial rate functions[END_REF][START_REF] Radulescu | Model reduction of biochemical reactions networks by tropical analysis methods[END_REF][START_REF] Kruff | Algorithmic reduction of biological networks with multiple time scales[END_REF][START_REF] Desoeuvres | Reduction of chemical reaction networks with approximate conservation laws[END_REF], and chained within a hybrid automaton [START_REF] Radulescu | Symbolic dynamics of biochemical pathways as finite states machines[END_REF][START_REF] Desoeuvres | Qualitative dynamics of chemical reaction networks: An investigation using partial tropical equilibrations[END_REF]. 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 [START_REF] Kruff | Algorithmic reduction of biological networks with multiple time scales[END_REF][START_REF] Desoeuvres | Reduction of chemical reaction networks with approximate conservation laws[END_REF].

Quasi-stable states were discussed for models in ecology [START_REF] Morozov | Long transients in ecology: Theory and applications[END_REF], but are certainly more general. A non-exhaustive list of nonlinear phenomena leading to quasi-stability include: i) slow attractive, normally hyperbolic invariant manifolds, that exist for slow/fast systems [START_REF] Wiggins | Normally hyperbolic invariant manifolds in dynamical systems[END_REF][START_REF] Kruff | Algorithmic reduction of biological networks with multiple time scales[END_REF][START_REF] Desoeuvres | Reduction of chemical reaction networks with approximate conservation laws[END_REF]; ii) saddle points with slow unstable directions, where slowness can be due to the proximity of a saddlenode bifurcation [START_REF] Morozov | Long transients in ecology: Theory and applications[END_REF]; and iii) ghost attractors, also due to the proximity of a saddle-node bifurcation [START_REF] Morozov | Long transients in ecology: Theory and applications[END_REF].

Quasi-stability of the reduced parametric model

We first look for slow invariant manifolds by using tropical geometry scaling methods introduced elsewhere [START_REF] Noel | Tropical geometries and dynamics of biochemical networks application to hybrid cell cycle models[END_REF][START_REF] Noel | Tropicalization and tropical equilibration of chemical reactions[END_REF][START_REF] Soliman | A constraint solving approach to model reduction by tropical equilibration[END_REF][START_REF] Radulescu | Tropical geometry of biological systems (invited talk)[END_REF]. By these methods we identify fast and slow variables and use quasi-steady state as zero order approximations to invariant manifolds [START_REF] Samal | A geometric method for model reduction of biochemical networks with polynomial rate functions[END_REF][START_REF] Radulescu | Model reduction of biochemical reactions networks by tropical analysis methods[END_REF][START_REF] Kruff | Algorithmic reduction of biological networks with multiple time scales[END_REF][START_REF] Desoeuvres | Reduction of chemical reaction networks with approximate conservation laws[END_REF].

For the sake of simplicity we rename the model variables and parameters from [START_REF] Alur | The benefits of relaxing punctuality[END_REF]. 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 where γ i round(log(k i )/ log(ϵ)) [START_REF] Samal | A geometric method for model reduction of biochemical networks with polynomial rate functions[END_REF][START_REF] Radulescu | Model reduction of biochemical reactions networks by tropical analysis methods[END_REF][START_REF] Kruff | Algorithmic reduction of biological networks with multiple time scales[END_REF][START_REF] Desoeuvres | Reduction of chemical reaction networks with approximate conservation laws[END_REF].

ẋ1 = k 1 x 1 1 - x 1 k 2 - k 3 x 2 k 4 + x 2 , ẋ2 = k 5 x 2 k 9 - x 2 k 6 - k 7 x 3 k 8 + x 3 , ẋ3 = k 10 (x 1 -k 11 x 3 ), (20)
For instance, for ϵ = 1/11, we have k1 = 0.5, γ 1 = 0, k2 = 1.86, γ 2 = -8, k3 = 9/10, γ 3 = 0, k4 = 0.933, γ 4 = -8, k5 = 0.5, γ 5 = 0, k6 = 1.27, γ 6 = -9, k7 = 5/4, γ 7 = 0, k8 = 1.76, γ 8 = 1, k9 = 5/11, γ 9 = 0, k10 = 1.21, γ 10 = 2, k11 = 2.46, γ 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 γ i . In the generic situation when the rescaled, zero order parameters ki are independent, γ i determine alone the orders of magnitude of the system's timescales [START_REF] Samal | A geometric method for model reduction of biochemical networks with polynomial rate functions[END_REF][START_REF] Radulescu | Model reduction of biochemical reactions networks by tropical analysis methods[END_REF][START_REF] Kruff | Algorithmic reduction of biological networks with multiple time scales[END_REF][START_REF] Desoeuvres | Reduction of chemical reaction networks with approximate conservation laws[END_REF]. However, when ki 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 ki . 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 ϵ ai 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 [START_REF] Samal | A geometric method for model reduction of biochemical networks with polynomial rate functions[END_REF][START_REF] Radulescu | Model reduction of biochemical reactions networks by tropical analysis methods[END_REF][START_REF] Kruff | Algorithmic reduction of biological networks with multiple time scales[END_REF][START_REF] Desoeuvres | Reduction of chemical reaction networks with approximate conservation laws[END_REF], early tropicalization ideas use rational ODEs [START_REF] Noel | Tropical geometries and dynamics of biochemical networks application to hybrid cell cycle models[END_REF]. 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:

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 (21)
Using the last equation we obtain the equivalent system 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 (22)
The system [START_REF] Kennedy | 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[END_REF] has two one dimensional branches of solutions

I. a 1 = -8, a 2 > -8, a 3 = 1.
II. 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

ẏ1 = k1 y 1 1 - y 1 k2 - k3 y 2 ϵ a2+8 k4 + y 2 ϵ a2+8 , ẏ2 = k5 y 2 k9 -ϵ a2+9 ȳ2 k6 - k7 y 3 k8 + y 3 , ẏ3 = ϵ -7 k10 (y 1 -k11 y 3 ). ( 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 / k11 . 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 ϵ a2+8 and ϵ a2+9 can be neglected. The 2D reduced order model reads

ẏ1 = k1 y 1 1 - y 1 k2 , ẏ2 = k5 y 2 k9 - k7 y 1 k8 k11 + y 1 . (24)
In the system (24) the first equation is decoupled from the second. In a short time, y 1 → k2 and then y 2 follows the equation

ẏ2 = k5 y 2 k9 - k7 k2 k8 k11 + k2 (25)
The condition for this to happen is that y 1 is faster than y 2 which reads 0 < k5 ( k9 -

k7 k2 k8 k11 + k2 ) << k1 (26)
Using the numerical values of the parameters we find k5 ( k9 -k7 k2 k8 k11 + k2 ) = 0.0790 k1 .

The variable y 2 is slow with a characteristic time k5 k9 -k7 k2 k8 k11 + k2 -1 = (0.5 × 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 ki one should verify our condition [START_REF] Kong | Skin microbiome: genomics-based insights into the diversity and role of skin microbes[END_REF] 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.

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 ki (condition [START_REF] Kong | Skin microbiome: genomics-based insights into the diversity and role of skin microbes[END_REF]).

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

α = k5 ( k9 - k7 k2 k8 k11 + k2 ), (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 k3y2ϵ a 2 +8 k4+y2ϵ a 2 +8 can no longer be neglected in the first equation of the system [START_REF] Koh | Skin microbiome of atopic dermatitis[END_REF]. This time depends on the initial value x 2 (0), scales like [log(k 4 ) -log(x 2 (0))]α -1 , and becomes infinite for x 2 (0) = 0.

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:

ẋi = r i x i (1 - x i K i - j̸ =i β ij x j θ ij + x j ), 1 ≤ i ≤ N, (28)
where x i are the species concentrations. We consider that the orders of magnitude of r i , α i , K i , β ij and θ 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, [START_REF] Korting | 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[END_REF] 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 ∩ S 2 = ∅, S 1 ∪ S 2 = {1, . . . , 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. β ij = 0 for i, j ∈ S 1 . Under these conditions, the dominating terms in the ODE model are ẋi = r i x i (1 -

x i K i ), for i ∈ S 1 , (29) ẋi
= r i x i (1 - j∈S1 β ij x j θ ij + x j ), for i ∈ S 2 . (30)
The approximate Eq.( 29) was obtained by neglecting the terms j∈S2 βij xj θij +xj , which is valid as long as x j << θ ij for all i ∈ S 1 , j ∈ 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-j∈S1 βij Kj θij +Kj ). Therefore, the quasi-stability conditions read: 0 < r i (1 -j∈S1 βij Kj θij +Kj ) << 1, for all i ∈ S 2 [START_REF] Lai | Commensal bacteria regulate tolllike receptor 3-dependent inflammation after skin injury[END_REF] x j (0) << θ ij , for all i ∈ S 1 , j ∈ S 2 [START_REF] Leyden | Propionibacterium levels in patients with and without acne vulgaris[END_REF] The non-negativity and the much-less-than conditions in [START_REF] Lai | Commensal bacteria regulate tolllike receptor 3-dependent inflammation after skin injury[END_REF] ensure that concentrations of species from S 2 increase, and that the increase is slow, respectively. The much-less-than condition in [START_REF] Leyden | Propionibacterium levels in patients with and without acne vulgaris[END_REF] ensures that the approximate equation ( 29) is initially valid. If [START_REF] Lai | Commensal bacteria regulate tolllike receptor 3-dependent inflammation after skin injury[END_REF], [START_REF] Leyden | Propionibacterium levels in patients with and without acne vulgaris[END_REF] are satisfied, then the saturated species reach almost constant levels x i = K i , i ∈ 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 θ ij in which case [START_REF] Kruff | Algorithmic reduction of biological networks with multiple time scales[END_REF] is no longer true. This lifetime reads min i∈S2,j∈S1 log(θ ij ) -log(x i (0)) r i (1 -l∈S1 β il K l θ il +K l ) .

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 quasistability 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 slowfast 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 quasistability, 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.

Figure 1 :

 1 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.

  ([Amp b ] = 0.32 µM ; [S p ] * = 5.10 8 CFU) and ([Amp b ] = 0.64 µM ; [S p ] * = 3.10 3 CFU), we get C ab = 0.16 and d spb r sp = 5 4 .

  and objective values x1 = 10 8 , y1 = 10 9 , x2 = 10 8 , y2 = 10 9 (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).

Table 6 :

 6 OAT and global sensitivity indices of the pathogenic model to variations of all kinetic parameters and initial concentrations around their nominal values rsp = 1, rsc = 0.5, [Amp h ] = 1.5, C 1 = 5e6, kc = 0.01, Ksc = 4e8, Ksp = 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 rsp, rsc, C ah , Ksc, Ksp, C ab , [Amp h ], [Sc] 0 , [Sp] 0 , C 1 , kc. 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%.

Figure 2 :

 2 Figure 2: Numerical simulation of the reduced ODE model over 50 hours, with initial conditions [Sc] = 10 5 , [Sp] = 10 3 , [Amp b ] = 0 and parameter values [Amp h ] = 1.5, rsc = 0.5, rsp = 1, C 1 = 5.10 6 , kc = 0.01, to fit Kohda et al. co-culture data [25] (Table2).

Figure 3 :

 3 Figure 3: Numerical simulation of the reduced ODE model over 50 hours, with initial conditions [Sc] = 10 5 , [Sp] = 10 3 , [Amp b ] = 0 and parameter values rsc = rsp = 0.5, C 1 = 2.10 8 , kc = 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.

Figure 4 :

 4 Figure 4: Landscape of satisfaction degree of the temporal formula corresponding to healthy skin with a skin surface pH of 5 (rsc = 0.5 and rsp = 0.3). The x and y axis represent variations of the initial quantities of [Sp] and [Sc], 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 , kc = 0.01, [Amp h ] = 3.

Figure 5 :

 5 Figure 5: Landscape of satisfaction degree of the temporal formula corresponding to compromised skin with a skin surface pH of 6.5 (rsc = 0.9 and rsp = 1.3). The x and y axis represent variations of the initial quantities of [Sp] and [Sc], 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 , kc = 0.01, [Amp h ] = 3.

Figure 6 :

 6 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 [Sp] and [Sc] respectively. The color coding corresponds to the satisfaction degree of the temporal logic formula. Values used for the other parameters: rsc = rsp = 0.5, C 1 = 2.10 8 , kc = 0.01.

Figure 7 :

 7 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.

Figure 8 :

 8 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. a) A typical orbit leaves the neighborhood of the unstable state (Sc = 0,Sp = 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 Sp(0) shows that for a long time Sp is constant whereas Sc increases exponentially but very slowly starting with very low values. The transition to the stable steady state occurs when Sc = 2 × 10 8 and the time spent in the quasi-stable state is given by log(2 × 10 8 /Sc(0))/α where Sc(0) and α are the initial values and effective exponential growth rate of the pathogen, respectively.

where x 1

 1 = [S c ], x 2 = [S p ], x 3 = [Amp b ], k 1 = r sc = 0.5, k 2 = K sc = 4 × 10 8 , k 3 = d sc /r sc = 9/10, k 4 = C 1 = 2 × 10 8 , k 5 = r sp = 0.5,k 6 = K sp = 3 × 10 9 , k 7 = d spb /r sp = 5/4, k 8 = C ab = 0.16, k 9 = 1 -d sph [Amp h ] r sp (C ah + [Amp h ]) = 5/11, k 10 = K c = 0.01, k 11 = d a /K c = 149 256 × 10 10 . The first step of tropical scaling is to rescale parameters as powers of a small positive parameter ϵ, k i = ki ϵ γi ,

Table 1 :

 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

	Variable	Interpretation (unit)
	[S c ]	Surface apparent concentration of S c (CF U.ASU -1 )
	[S p ]	Surface apparent concentration of S p (CF U.ASU -1 )
	[Amp b ]	Concentration of Amp b (AU.ASU -1 )
	Parameter Interpretation (unit)
	r sc	Growth rate of S c (h -1 )
	r sp	Growth rate of S p , (h -1 )
	K sc	Optimum concentration of S c (CF U.ASU -1 )
	K sp	Optimum concentration of S p (CF U.ASU -1 )
	d sc	Maximal killing rate of S c by S p (h -1 )
	C 1	Concentration of S p inducing half the maximum killing rate d sc
		(CF U.ASU -1 )
	d spb	Maximal killing rate of S p by Amp b , (h -1 )
	C ab	Concentration of Amp b inducing half the maximum killing rate
		d spb (AU.ASU -1 )
	d sph	Maximal killing rate of S p by Amp h , (h -1 )
	C ah	Concentration of Amp h inducing half the maximum killing rate
		d sph (AU.ASU -1 )
	[Amp h ]	Concentration of AMPs produced by the skin cells (AU.ASU -1 )
	k c	Production rate of Amp b by S c (AU.h -1 .CF U -1 )
	d a	Degradation rate of Amp b (AU.h -1 )

Table 2 :

 2 Experimental data from Kohda et. al[START_REF] Kohda | An in vitro mixed infection model with commensal and pathogenic staphylococci for the exploration of interspecific interactions and their impacts on skin physiology[END_REF] used for identifying parameter values.

		S. epidermidis (CFU/well) S. aureus (CFU/well)
	Mono-cultures	4.10 8	3.10 9
	Co-cultures	1.10 8	1.10 9

Table 3 :

 3 Experimental data from Nakatsuji et. al[START_REF] Nakatsuji | Antimicrobials from human skin commensal bacteria protect against staphylococcus aureus and are deficient in atopic dermatitis[END_REF] used for identifying parameter relations.

	Sh-lantibiotics (µM ) LL-37 (µM ) S. aureus (CFU/mL)
	0	4	10 9
	0	8	6.10 5
	0.32	0	5.10 8
	0.64	0	3.10 3

Table 4 :

 4 Summary of the parameter relations embedded in the reduced model.

	Parameter Value or relation to other parameters
	K sc		4.10 8
	K sp		3.10 9
	C ah		8
	C ab		0.16
	d sph		2 r sp
	d spb		5 4 r sp
	d sc	r sc	3 4.10 9 C 1 +	3 4
	d a	10 8 k

c 56 + 31[Amp h ] 2.56 (4 -[Amp h ]) with [Amp h ] < 4

  )

	2 All	computation	results	reported	in	this	article	have	been	ob-
	tained	with	a	Biocham	notebook	available	in	different	formats	at
	https://lifeware.inria.fr/wiki/Main/Software#TCS23.

Table 5 :

 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, rsp = 1. The parameter values retained are in bold.

	C 1 1e3 1e4 1e5	1e6 5e6 1e7	1e8	1e9
	r sc 1.3	1	0.76 0.59 0.5 0.46 0.38 0.34

Table 7 :

 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 predominance 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 , Ksc, [Amp h ], rsc 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.

	Parameters	Nominal value COV Sensitivity index
	C ah	8	0.2	0.242
	C ab	0.16	0.2	0.207
	K sc	4.10 8	0.2	0.188
	[Amp h ]	2.5	0.2	0.187
	r sc	0.5	0.2	0.007
	K sp	3.10 9	0.2	0.005
	[S p ] 0	10 3	0.2	0.005
	[S c ] 0	10 5	0.2	0.005
	r sp	0.5	0.2	0.005
	C 1	2.10 9	0.2	0.005
	k c	0.01	0.2	0.005
	(r sp , r sc , [Amp h ], C 1 , k c )	"	0.2	0.184
	tuple of all 11 parameters	"	0.2	0.321

Table 8 :

 8 Experimental data from Dasgupta et. al[START_REF] Dasgupta | 16502 effect of ph on growth of skin commensals and pathogens[END_REF] showing the influence of pH on growth rates of S. epidermidis and S. aureus

	pH	Growth rate (∆OD/hour)
		S. aureus S. epidermidis
	5	0.03	0.05
	5.5	0.04	0.07
	6	0.09	0.08
	6.5	0.13	0.09
	7	0.14	0.10

been observed. Dasgupta et al. studied in vitro the influence of pH on the growth rates of S. aureus and S. epidermidis

The Biocham notebook mentioned above shows the detailed results of parameter search and simulation figures on logscale.

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.