Bifurcation analysis of a mathematical model of atopic dermatitis to determine patient-specific effects of treatments on dynamic phenotypes

Atopic dermatitis (AD) is a common inflammatory skin disease, whose incidence is currently increasing worldwide. AD has a complex etiology, involving genetic, environmental, immunological, and epidermal factors, and its pathogenic mechanisms have not yet been fully elucidated. Identification of AD risk factors and systematic understanding of their interactions are required for exploring effective prevention and treatment strategies for AD. We recently developed a mathematical model for AD pathogenesis to clarify mechanisms underlying AD onset and progression. This model describes a dynamic interplay between skin barrier, immune regulation, and environmental stress, and reproduced four types of dynamic behaviour typically observed in AD patients in response to environmental triggers. Here, we analyse bifurcations of the model to identify mathematical conditions for the system to demonstrate transitions between different types of dynamic behaviour that reflect respective severity of AD symptoms. By mathematically modelling effects of topical application of antibiotics, emollients, corticosteroids, and their combinations with different application schedules and doses, bifurcation analysis allows us to mathematically evaluate effects of the treatments on improving AD symptoms in terms of the patients’ dynamic behaviour. The mathematical method developed in this study can be used to explore and improve patient-specific personalised treatment strategies to control AD symptoms. © 2018 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Atopic dermatitis (AD), also known as atopic eczema, is one of the most common chronic skin diseases, whose prevalence is around 10-20% in developed countries and is also rapidly increasing in developing countries ( Flohr and Mann, 2014;Guo et al., 2016;Weidinger and Novak, 2016 ). AD patients typically suffer from chronic or recurrent skin inflammation and intense itch in affected skin areas. Socioeconomical burdens of AD, together with degraded quality of life for patients and caregivers, is a grave healthcare concern ( Carroll et al., 2005 ). However, pathological mechanisms of AD still remain to be fully clarified. While the cur-rent mainstay of AD treatment aims at controlling AD symptoms by topical therapies ( Eichenfield et al., 2017;, it is crucial to refine and develop personalised therapeutic strategies for AD control ( Bieber et al., 2017 ), as well as to develop new treatment approaches ( Flohr and Mann, 2014 ), based on the systematic understanding of the highly complex interplay between genetic, environmental, immunologic, and epidermal factors that govern AD pathogenesis.
Mathematical modelling and analysis are potent tools to systematically investigate pathogenic processes, thereby exploring better prevention and treatment measures against AD. We recently developed a mathematical model of AD pathogenesis, based on empirical evidence from clinical and experimental studies . This multi-scale hybrid model describes interacting dynamics of major regulatory factors of AD, including skin barrier function, immune responses, and environ- mental stressors. The model behaviour reproduced a variety of clinical and biological data from AD patients and AD animal models, confirming the validity of the model. Model parameters were obtained by fitting the model to in vivo experimental data where possible.
Our model analysis revealed four qualitatively different dynamic phenotypes (recovery, chronic damage, bistability, and oscillation dynamic phenotypes) that appear in response to environmental stressors : The recovery dynamic phenotype corresponds to the healthy case where the system recovers to the healthy steady state quickly and spontaneously without any treatments; the chronic damage dynamic phenotype corresponds to the severe disease case where the system converges to the unhealthy steady state; the bistability and the oscillation dynamic phenotypes are the asymptomatic-but-susceptible cases where appropriate treatments are required to avoid the transition to a more severe phenotype . The dynamic phenotype thus reflects the patient's AD severity, including the vulnerability to environmental stressors and the need for treatments to prevent or revert aggravation of AD symptoms, and varies depending on the model parameters that characterise each patient's genetic, environmental and therapeutic conditions.
In order to systematically explore the patient-specific effects of different therapeutic interventions, here we derive mathematical conditions that determine whether and how the system can be perturbed by treatments to transition to a more physiological dynamic phenotype, for example, from the asymptomatic-butsusceptible oscillation to the healthy recovery dynamic phenotypes. By modelling the effects of three major treatments for AD, namely topical application of antibiotics ( Lin et al., 2007 ), emollients ( Simpson et al., 2014 ), and corticosteroids ( Leyden et al., 1974;Nilsson et al., 1992 ), we evaluate the effects of different treatments, at different doses and combinations, on dynamic phenotypes, via bifurcation analysis of the model of AD pathogenesis. Specifically, we obtain analytical parameter conditions for transitions between different dynamic phenotypes. Analytical parameter conditions, described in terms of the changes in the parameter values, or of the strength, timing, and duration of treatments to be applied, allow us to assess the patient-specific impacts of different treatment strategies systematically.
The method developed in this study, based on bifurcation analysis, enables us to systematically stratify the patients with different dynamic phenotypes and to evaluate patient-specific optimal treatment strategies (the types and combinations of drugs and application schedules) that improve prevention and control of AD by achieving preferable transition of dynamic phenotypes.

Model description
We first introduce the mathematical model used in this study ( Fig. 1 (a)). The model was originally developed by Domínguez-Hüttinger et al. (2017) to describe an interplay between skin barrier function, immune responses, and environmental stressors, that governs AD pathogenesis. The model was then expanded by Christodoulides et al. (2017) to include the effects of treatment by topical application of corticosteroids and emollients. In this study, we also consider the effects of antibiotics ( Lin et al., 2007 ). Here we show the mathematical description of the model and provide a brief explanation of the biological background, while the detailed biological explanation of the model and its relevance to the experimental and clinical data is found in the original papers Domínguez-Hüttinger et al., 2017 ).
The model ( Fig. 1 ) is given by where the three continuous state variables, namely, the level of pathogen load (or the amount of infiltrated environmental stressors in general), P ( t ), the strength of skin barrier integrity, B ( t ), and the concentration of activated dendritic cells (or the level of inflammation markers in general), D ( t ), characterise the disease status at time t (day). The levels of activated innate immune receptors, R ( t ), activated kallikreins, K ( t ), and Gata-3, G ( t ), are described as discontinuous state variables. A ( t ), E ( t ), and C ( t ) correspond to the amounts of antibiotics, emollients, and corticosteroids applied, respectively. The state variables and the model parameter values are listed in Tables 1 and 2 .
The main mechanisms described in the model are as follows ( Fig. 1 (a)). Stimulations with environmental stressors ( P env ) increase the pathogen load ( P ( t )), which activates the innate immune receptors ( R ( t )) and kallikreins ( K ( t )). Activation of innate immune receptors and kallikreins weakens the skin barrier by inhibiting its formation and inducing its desquamation/degradation. This makes the skin barrier more vulnerable, allowing an increase in the pathogen level. The activation of the innate immune receptors also triggers immune responses, eradicating the infiltrated pathogens but causing inflammation and thus AD flares. The activated immune responses further reduce the barrier integrity by activating dendritic cells, which then move to the lymph nodes, where they induce the Gata-3 mediated differentiation and subsequent recruitment of pro-inflammatory T helper 2 (Th2) cells from the lymph nodes to the epidermis.
We assume that the levels of activated innate immune receptors, R ( t ), and of active kallikreins, K ( t ), are determined by a reversible switch under the variation of P ( t ) ( Fig. 1 (b)). They are described by where R on , R off , and K off are constant and K on = m on P (t) − β on depends on P ( t ). If an increasing P ( t ) exceeds an activation threshold value, P + , then the switch is activated and AD flares occur. If a decreasing P ( t ) falls below an inactivation threshold value, P − , then the switch is inactivated and the AD flares cease. This hysteretic switching is an approximate representation of the jump between equilibrium states in a cellular-level model of AD ( Tanaka et al., 2011 ).
We further assume that Gata-3 is activated once the concentration of the dendritic cells in the lymph nodes, D ( t ), increases beyond a threshold value, D + , and that the level of Gata-3 activation, G ( t ), is maintained even if D ( t ) decreases ( Fig. 1 (c)). G ( t ) is therefore described by an irreversible switch, where G on and G off are constant. This irreversible switch is an approximate representation of the jump between the equilibrium

Fig. 1.
Schematic illustration of the mathematical model for AD pathogenesis with drug treatments Domínguez-Hüttinger et al., 2017 ). (a) The model overview. (b) Hysteretic reversible switch for the level of activated immune receptors ( R ( t )) and kallikreins ( K ( t )), with the activation and inactivation threshold values of the pathogen load, P + and P − , respectively. (c) Irreversible switch for the level of Gata-3 activation ( G ( t )), with the activation threshold value of the concentration of dendritic cells, D + . states in the cellular-level model of the T cell differentiation developed in Ref. Höfer et al. (2002) . The switching-on of the Gata-3 activation corresponds to the progression from a mild to more severe AD phenotype. In this study, we assume G ( t ) to be constant, either at G off or G on , as we focus on the bifurcation diagram for each of these cases, rather than the progression of AD. For a constant G ( t ), the dynamics of P ( t ) and B ( t ) in Eqs. (1) and (2) is independent of that of D ( t ) in Eq. (3) , and the original system is reduced to a twodimensional system ( Eqs. (1 ) and (2) ) combined with a reversible switch ( Eq. (4 )).
Our model also considers effects of three widely-used treatments for AD, namely antibiotics ( A ( t )), emollients ( E ( t )), and corticosteroids ( C ( t )) ( Fig. 1 (a)). Antibiotics inhibit the growth of the pathogen load, emollients have a moisturizing and occluding effect on the skin barrier, and corticosteroids block the activation of immune receptors and Gata-3. All these treatments enhance the barrier function either directly (emollients) or indirectly by damping the barrier-damaging effects of inflammation (antibiotics and corticosteroids). The efficacy of corticosteroids is described by the coefficients, β 1 , β 2 , β 3 , and β 4 , which correspond to the four different effects of corticosteroids on the epidermis, namely pathogen eradication, innate immunity-mediated barrier production, adaptive immunity-mediated barrier production, and activation of dendritic cells, respectively.

Four dynamic phenotypes
Depending on the choice of the model parameters, the AD model described by Eqs. (1) and (2) without treatment ( A (t) = E(t ) = C(t ) = 0 ), coupled with Eq. (4) , demonstrates four dynamic phenotypes typically observed in AD patients ( Fig. 2 )  . The dynamic phenotype characterises the dynamical behaviour of the system in response to an increase in the level of environmental stressors ( P env ) or a transient increase in the pathogen load level ( P (0)), and is classified by the equilibrium state to which the system settles down. The dynamic phenotypes can be used to stratify patients in terms of the treatments that they require to avoid worsening of the AD symptoms.
The four dynamic phenotypes are as follows: (i) Recovery dynamic phenotype ( Fig. 2 (a)): Even if the skin barrier level is degraded temporarily due to an increase in the level of environmental stressors, it eventually recovers to a homeostatic level after the pathogen load decreases below P − , and then converges to a stable value that is kept below P + , and thus R (t) = R off . Treatment is not necessary to keep the healthy phenotype, as the trajectory converges to the stable equilibrium, H , corresponding to the healthy state without AD flares. (ii) Chronic damage dynamic phenotype ( Fig. 2 (b)): The skin barrier integrity is almost lost and the pathogen load is kept at a high level above P − , and thus R (t) = R on . A constant or periodic treatment is necessary to mitigate severe AD symptoms. Without treatment, the trajectory converges to the stable equilibrium, U , corresponding to the unhealthy state with persistent AD flares. (iii) Bistability dynamic phenotype ( Fig. 2 (c)): Either the healthy state or the unhealthy state appears depending on the initial condition. If the unhealthy state appears, a transient treatment can induce convergence to the healthy state by perturbing, or deviating, the unhealthy trajectory into the basin of attraction of H , as the trajectory converges to one of the coexisting stable equilibria, H (healthy) or U (unhealthy). This dynamic phenotype corresponds to an "asymptomaticbut-susceptible" phenotype that can easily move to a more severe phenotype . (iv) Oscillation dynamic phenotype ( Fig. 2 (d)): Remission and relapse of AD flares are alternately repeated. The trajectory goes back and forth between the two switching boundaries,  phenotype by the dynamics of the skin barrier integrity, which demonstrates a longer recovery time (for the healthy branch in the bistability dynamical phenotype) or a higher variability over time (for the oscillation dynamical phenotype). The chronic damage dynamical phenotype, corresponding to severe symptoms, can be identified by conventional clinical measurements for AD severity, such as SCORing Atopic Dermatitis (SCORAD).

Stability and bifurcation analysis
To clarify the mathematical conditions required for the transition of dynamic phenotypes in the AD model ( Eqs. (1) , (2) and (4) ), we first analytically derive the equilibrium points that characterise the recovery, chronic damage and bistability dynamic phenotypes, and present a method to computationally calculate the periodic solution characterising the oscillation dynamic phenotype. The results described in this section for the general case will then be used to investigate the transition of dynamic phenotypes for different scenarios in the next sections.

Equilibrium points
We analytically derive the two stable equilibria, H and U , corresponding to the healthy and unhealthy states, respectively, by assuming a constant treatment over time, i.e. A (t) = Ā , E(t ) = Ē , and C(t ) = C , for simplicity. The system is described by the two- ˙ where c j ≡ 1 / (1 + β j C ) for j = 1 , 2 , and 3 (introduced for convenience of notation) are inversely related to the efficacy of corticosteroids on pathogen eradication ( β 1 ), innate immunity-medicated barrier production ( β 2 ), and adaptive immunity-mediated barrier production ( β 3 ), respectively. The value of c j ( j = 1 , 2 , and 3 ) is equal to one when corticosteroids are absent ( C = 0 ) or not effective at all ( β j = 0 ), and decreases with an increase in their efficacy ( β j ). The equilibrium point (P, where R * , K * and G * are either the on or off state of R ( t ), K ( t ) and G ( t ) at the equilibrium point, respectively. The Jacobian matrix of Eqs. (6) and (7) at the equilibrium point ( P * , B * ) is given by where K = 0 if the reversible switch is off and K = m on otherwise.
Noting that B * ∼ 0 if the reversible switch is on, the eigenvalues of J are approximated by the diagonal elements, both of which are negative. Hence, the equilibrium point is stable, if it exists, according to the linear stability theory for dynamical systems ( Hirsch et al., 2012 ). The features of the equilibrium point vary depending on the state of the reversible switch, either off with ( R off , K off ) or on with ( R on , K on ), but are independent of the state of the irreversible switch, either on with G on or off with G off , as described below.
with G * being either G off or G on . This equilibrium point exists if 0 < P H < P + , i.e.
For instance, as κ P increases from a sufficiently small value, the stable equilibrium H collides with the switching boundary + = { (P, B ) | P = P + } and disappears through a boundary equilibrium (BE) bifurcation ( Bernardo et al., 2008 ). The effect of switching between the two states of G * on B H is equivalent to a scaling of κ B .
Therefore, it is not necessary to explicitly analyse the effects of switching on G on the behaviour of the system.
with G * being either G off or G on . This equilibrium point exists only when P U > P − , i.e.
As κ P decreases from a sufficiently large value, the equilibrium U collides with the switching boundary − = { (P, B ) | P = P − } and disappears at κ P = κ − P via a BE bifurcation. Classification into the recovery, chronic damage, and bistability dynamic phenotypes is determined by the equilibrium-point analysis that can evaluate global long-term behaviour of the system by the changes in the stability and existence of focal points, ( P H , B H ) and ( P U , B U ) ( Fig. 3 (a)). The recovery dynamic phenotype is observed when P H < P + and P U < P −, the chronic damage dynamic phenotype when P H ≥ P + and P U ≥ P −, and the bistability dynamic phenotype when P H < P + and P U ≥ P −. The remaining region with P H ≥ P + and P U < P − roughly corresponds to the oscillation dynamic phenotype. In general, however, the precise condition for the oscillation dynamic phenotype, corresponding to a periodic solution, is not specified solely by the equilibrium-point analysis.
In the next subsection, we describe the numerical method to reveal the necessary condition for observing the oscillation dynamic phenotype.

Periodic solution
Here we analyze the stability and bifurcations of a periodic solution for our model, to reveal the necessary condition for observing the oscillation dynamic phenotype, which corresponds to a stable periodic solution.
Let us focus on a trajectory that travels back and forth between the switching boundaries, − and + . For the k -th ( k = 1 , 2 , . . . ) segment of the trajectory that goes from the point (P, Fig. 3 (b)), we define discrete-time maps, T f and T b , as follows ( Kousaka et al., 1999;Tanaka et al., 2008 ): Using the Poincaré map T ≡ T b •T f , given by the composite of T f and T b , the time evolution of B − k can be represented by The fixed point of the Poincaré map, corresponding to a periodic solution in the original model, is obtained by numerically solving where μ is the eigenvalue characterising the type of local bifurcation ( Kuznetsov, 2013 ), i.e. μ = 1 for a saddle-node (SN) bifurcation point, μ = −1 for a period-doubling bifurcation point, and μ = e iθ with θ = 0, π for a Neimark-Sacker bifurcation point.
Solving Eqs. (21) and (22) simultaneously using Newton's method allows us to specify the critical parameter values for the local bifurcation of a periodic solution ( Kawakami, 1984 ) and reveal the conditions of the parameters for which the periodic solution is stable. When a stable periodic solution coexists with other stable equilibrium solutions, we assume that the model behaviour corresponds to the oscillation dynamic phenotype for simplicity.

Baseline bifurcation diagrams
In this section, we obtain the baseline bifurcation diagrams for our model with the nominal parameter values ( Table 2 ) without any treatment ( A (t) = E(t ) = C(t ) = 0 ). We show both 1-parameter and 2-parameter bifurcation diagrams ( Fig. 4 ) to achieve a better understanding of the system. The 2-parameter bifurcation diagram will be further used in the next Section to evaluate the conditions for transition of the dynamic phenotypes, induced by changes in the values of the model parameters or the treatments.
The bifurcation diagrams are obtained in this and next Sections by deriving the bifurcation sets for the dynamic phenotypes as described in the previous section. We confirm our results by a brute force method, i.e. by numerical integration of Eqs. (6) and (7) for each bifurcation parameter(s) and different initial conditions. Specifically, we identify the dynamic phenotype corresponding to each bifurcation parameter(s) depending on the steady states obtained from three initial conditions ( P (0), B (0), D (0)): (10, B H , 0) close to the equilibrium point H , (50,0,0) close to the equi-librium point U , and (33, B H /2, 0) close to the periodic solution O . For each initial condition, we numerically integrate the model using a fourth-order Runge-Kutta method with a time step of 0.005 until t = 10 0 0 and observe the time series after t = 500 to exclude transient behaviour. Each of these three steady states are classified as a periodic solution if P ( t ) oscillates between P − and P + ; as an equilibrium point H if the solution is not oscillatory and B (t) ≥ B H − where = 0 . 1 ; and as an equilibrium point U if the solution is neither a periodic solution nor the equilibrium H . The dynamic phenotype is then determined as recovery if the steady state is H for all the three initial conditions, oscillation if a periodic solution is obtained, bistability if the steady state can be H or U depending on the initial condition, and chronic damage if the steady state is U for all the three initial conditions.

1-parameter bifurcation diagram
The skin permeability is an important factor whose increase may increase the risk of sensitization to allergens. With an increase in the nominal skin permeability, κ P , we observe a transition from the recovery to the oscillation and subsequently to the chronic damage dynamic phenotype ( Fig. 4 (a)).
Within the parameter range for the recovery phenotype, the trajectory of the system converges to the healthy state H with B H ∼ 1, given by Eqs. (11) and (12) . At κ P = κ + P ∼ 0 . 842 , H collides with the switching boundary + and disappears through the BE bifurcation leading to a transition to an oscillatory regime ( Bernardo et al., 2008 ). When κ P is within the range of the oscillation dynamic phenotype, the trajectory approaches the stable periodic solution corresponding to the fixed point of the Poincaré map T defined by Eq. (21) . The bifurcation analysis in Section 3.2 reveals that the stable periodic solution disappears via saddle-node (SN) bifurcation at both end points ( Fig. 4 (a) inset for the right one). At κ P = κ − P ∼ 1 . 265 , the unhealthy state U with B U ∼ 0 given by Eqs. (14) and (15) appears through a BE bifurcation when it collides with the switching boundary − . We assume that a parameter range in which the stable periodic solution and the unhealthy state coexist corresponds to the oscillation dynamic phenotype. After the disappearance of the stable periodic solution at the SN bifurcation, the trajectory globally converges to U , resulting in the chronic damage dynamic phenotype.

2-parameter bifurcation diagram
A reduced ability to eradicate pathogen by weak immune responses is another important risk factor governing the severity of AD. The dynamic phenotype for each patient is determined in our model by the balance between the level of skin permeability and the strength of the immune responses that eradicate the pathogens that have infiltrated the epidermis ( Fig. 4 (b)). The parameter space for ( κ P , α I ), where α I is the pathogen eradication rate, is separated into the regions corresponding to the four dynamic phenotypes, with their boundaries obtained by our bifurcation analysis in the previous section and corroborated based on the brute-force method.
The BE bifurcation set for the equilibrium point H (the vertical thin solid line ranging for 0 < α I < 0.5) is determined by Eq. (13) , which is independent of α I if R (t) = R off = 0 . The BE bifurcation set for the equilibrium point U (the dashed line) is given by Eq. (16) . The left-bottom region surrounded by these BE bifurcation sets corresponds to the bistability dynamic phenotype (B). When ( κ P , α I ) are in the parameter region corresponding to B, the stable equilibrium points H and U coexist and the state space is divided into two basins of attraction with respect to the equilibria, H and U ( Fig. 4 (c)), suggesting that a temporal treatment that moves the physiological state into the basin of attraction of H would achieve recovery. The SN bifurcation sets for the periodic solution O (the bold solid curves with a cusp) are numerically obtained by solving Eqs. (21) and (22) . The two SN bifurcation curves collide at the cusp point. In the vicinity of the cusp point, three attracting states, H, U , and O , coexist, and the state space is divided into three basins of attraction ( Fig. 4 (d)). Given that the basin of H is relatively small, compared to other basins, it appears to be difficult to move the trajectory into the small basin of H and achieve a recovery by a transient treatment, which could nonetheless help avoid the unhealthy state and settle the system into the oscillation state. Such treatments can be rationally designed based on our mathematical methods.

Transition of dynamic phenotypes
In this section, we obtain the bifurcation diagrams for different choices of the parameter set (both for the model parameters and the treatments) of our model, and evaluate when the transition of dynamic phenotypes occurs. We denote by XY the transition from a dynamic phenotype X to Y, where X and Y are the symbols representing the dynamic phenotype (R, O, B, or D) ( Fig. 5 (a)). Among 16 transitions between dynamic phenotypes five transitions (BR, OR, DR, DB and DO) correspond to improvement of the AD severity and are favourable for patients ( Fig. 5 (b)).
As shown in the previous sections, the boundaries of the parameter regions for different dynamic phenotypes are determined by the BE bifurcation sets for the equilibrium points ( Eqs. (13) and (16) ) and the SN bifurcation sets for the periodic solutions that are numerically specified. By assessing how these boundaries change from those for the baseline bifurcation diagram, we evaluate the effects of changes in model parameters (reflecting the patient's specificity) or of application of treatments on the dynamic phenotypes. Specifically, each region separated by the bifurcation sets for the baseline and new bifurcation diagrams can be associated with a transition XY ( Fig. 5 (c)), demonstrating that the patients with the dynamic phenotype X will have the dynamic phenotype Y after changes in model parameters or application of treatments. We also approximately evaluate the likelihood of the transition XY, by the ratio of the area for the XY region in the new bifurcation diagram to that for the X region in the baseline bifurcation diagram, and represent it by the thickness of the arrow for the transition XY in the corresponding transition diagram ( Fig. 5 (b)).

Transition of dynamic phenotypes by changes in model parameters
We first assess the effects of changes in some model parameters whose variations can be caused by differences in the microbiota composition of the epidermis ( Byrd et al., 2017;Meisel et al., 2018 ), which is one of the recent therapeutic targets for AD  ).
An increase in the basal pathogen death rate, δ P , reduces the net growth rate of P ( t ) ( Eq. (6) ), resulting in the right shift of the boundaries between the regions for H and U , enlarging the parameter regions for the recovery and bistability dynamic phenotypes ( Fig. 6 (a)). A decrease in P env results in a similar effect on the bifurcation diagram ( Fig. 6 (b)), but the slope of the boundary line for U changes as seen from Eq. (16 ). Therefore, the transition from the chronic damage dynamic phenotype is more likely for a larger value of κ P . An increase in the rate of barrier-mediated inhibition of pathogen infiltration, γ B , shifts the boundary for H to the right but barely shifts that for U ( Fig. 6 (c)), because B H > 0 in Eq. (13) and B * U ∼ 0 in Eq. (16) .  when (a) δ P is increased from 1 (nominal value) to 1.6, (b) P env is decreased from 95 to 55, and (c) γ B is increased from 1 to 1.5. Other parameters are set at the nominal values ( Table 2 ). In the bifurcation diagrams, the arrows indicate the shift of the bifurcation sets (thick lines) from those for the baseline bifurcation diagram (thin lines) induced by the relevant parameter change: for the stable healthy equilibrium (solid lines) and the stable unhealthy equilibrium (dashed lines). XY (X,Y = R, O B, D) represents the region corresponding to the transition of dynamic phenotypes from X to Y. These plots are for the cases with G (t) = G on . Similar plots can be obtained for G (t) = G off .
These parameter variations show improvement of the AD severity (BR, OR, DR, DB and DO) and are favourable for patients.

Transition of dynamic phenotypes by treatments
By a similar bifurcation analysis, we evaluate whether drug treatments with antibiotics, emollients, and corticosteroids can induce favourable transitions of dynamic phenotypes (BR, OR, DR, DB and DO) when used separately or in combination at different schedules.

Constant treatment with a single drug
We first consider the treatment with a single drug at a constant amount, i.e. (A (t) , E(t ) , C(t )) = ( Ā , 0 , 0) , (0 , Ē , 0) , or (0 , 0 , C ) , and derived the bifurcation and transition diagrams ( Fig. 7 ). When antibiotics or emollients are applied (i.e. Ā or Ē is non-zero), the parameter regions for the recovery and bista-  Tanaka et al. (2011) between P ( t ) and K on βon Y-intercept of the linear relation 6.71 Tanaka et al. (2011) between P ( t ) and K on bility dynamic phenotypes are enlarged, while the regions for the chronic damage and oscillation dynamic phenotypes become smaller ( Fig. 7 (a), (b), (d) and (e)). These numerical results are concordant with the analytical results, demonstrating a right shift in the boundaries of the regions for H and U as Ā or Ē increases ( Eqs. (13 ) and (16) ), as a result of inhibition of the pathogen growth by antibiotics ( Eq. (6) ) or strengthening of the skin barrier integrity by emollients ( Eq. (7) ). These results suggest that a sufficient amount of antibiotics or emollients can, in principle, induce the recovery dynamic phenotype for the patients by BR, OR and DR, although side effects and drug resistance of antibiotics should be considered when a therapeutic intervention with antibiotics is to be designed. Our theoretical analysis thus guarantees the benefits of the use of antibiotics and emollients in our model. However, corticosteroids are not always effective in improving the disease state, as they have both positive and negative effects on the AD pathogenesis in our model. The positive effects of corticosteroids are induced by the increased skin barrier production through innate immunity ( β 2 ) and adaptive immunity ( β 3 ). The negative effects of corticosteroids are represented as the weakened eradication of infectious load ( β 1 ), both directly via inhibition of T cell cytokine transcription and indirectly via skin atrophy that allows more pathogens to penetrate through the skin. The degree of these negative effects are known to vary depending on the potency of the corticosteroids and the body sites to which the treatments are applied ( Eichenfield et al., 2014 ). For example, low potency corticosteroids applied to less sensitive areas (e.g. arms or legs compared to face) are considered to have much smaller negative effects. Here we consider two cases, Case 1 with both the positive and negative effects of corticosteroid included ( β 1 > 0) and Case 2 with only the positive effects ( β 1 = 0 ) as an extreme case without negative effects. Com parison between these two extreme cases allows us to evaluate the two ends at the spectrum of corticosteroid efficacy, instead of investigating the different degrees of the balance between the positive and negative effects of corticosteroids on the transition of dynamic phenotypes.
In Case 1 ( β 1 > 0), the parameter region for the recovery dynamic phenotype is reduced and that for the chronic damage dynamic phenotype is enlarged as C increases ( Fig. 7 (c)). This is analytically shown by the left shift of the boundaries of the regions for H and U , which is mainly caused by the decrease in c 1 = 1 / (1 + β 1 C ) ( Eqs. (13 ) and (16) ). In Case 2 ( β 1 = 0 ), the beneficial role of corticosteroids appears as an enlargement of the parameter region for the recovery dynamic phenotype ( Fig. 7 (f)). It is due to the effects of corticosteroid treatments on increasing the skin barrier level ( Eq. (12) ) and shifting the boundaries of the regions for H and U to the right ( Eqs. (13 ) and (16) ). These analyses demonstrate that corticosteroids can have positive or negative treatment effects depending on the intrinsic property of patients (that are reflected in the model parameter values) and the efficacy parameters for the actions of corticosteroids ( β j for j = 1 , 2 , 3 ), suggesting the importance of personalised selection of the treatment method.

Intermittent treatment with a single drug
We next consider intermittent application of a single drug, as is often recommended in therapeutic interventions for AD ( Williams, 2005 ). The evaluation of appropriate frequency of treatment has been a focus of research to achieve both clinical and cost effectiveness ( Green et al., 2005 ).
We assume an intermittent treatment with T X (0 ≤ T X < 7) days of on-treatment followed by (7 − T X ) days of off-treatment every week, defined by where X = A, E , or C corresponding to the type of treatment applied, i.e. antibiotics, emollients, or corticosteroids, respectively. The operation "a (mod b )" indicates the remainder after division of a by b . We examined how the duration of the on-treatment period ( T C ) and the amount of treatment ( C ) for corticosteroids affect the dynamic phenotypes, for the case where corticosteroids are beneficial Fig. 7. Effects of continuous treatment on the bifurcation diagram, with antibiotics (left panels), emollients (middle panels), and corticosteroids (right panels). Bifurcation diagrams (top) and observed transitions of dynamic phenotypes (bottom) are shown for each case. (a)-(c) Case 1 with G (t) = G on , β 1 = β 2 = β 3 = β 4 = 1 , δ P = 1 , R off = 0 , and K off = 0 , and (d)-(f) Case 2 with β 1 = 0 , β 2 = β 3 = β 4 = 10 , δ P = 1 . 6 , R off = 1 , and K off = 0 . 1 . In the bifurcation diagrams, the arrows indicate the shift of the bifurcation sets (thick lines) from those for the baseline bifurcation diagram (thin lines) induced by the relevant treatment: for the stable healthy equilibrium (solid lines) and the stable unhealthy equilibrium (dashed lines). XY (X,Y = R, O B, D) represents the region corresponding to the transition of dynamic phenotypes from X to Y.
(Case 2 with β 1 = 0 ). Our simulation results demonstrate that transition to the recovery dynamic phenotype is effectively achieved by a certain amount of total corticosteroids applied ( C × T C ) with a tradeoff between C and T C ( Fig. 8 ). Specifically, a low C for a long T C and a high C for a short T C cannot achieve the recovery dynamic phenotype, or remission with a high B . Similar results were obtained for antibiotics and emollients (Fig. S1 ) and for Case 1 (Fig. S2 ).
These analyses suggest that intermittent treatment can induce and maintain the healthy state if the amount of drug and the length of the on-treatment period are appropriately adjusted. The fact that many options of treatment types and schedules can achieve an equivalent effect in terms of the transition of dynamic phenotypes suggests a possibility to tailor the treatment schedule to the needs and convenience of each patient without compromising the treatment efficacy.

Continuous treatment with a combined use of drugs
AD treatments can become more effective by a concomitant use of multiple drugs targeting different factors responsible for the disease. The bifurcation diagrams for different amounts of two drugs under continuous treatment ( Fig. 9 ) provide the information on how to balance the amounts of the two drugs to induce a transition to a more physiological dynamic phenotype. Similarly for the baseline bifurcation diagram, the boundaries separating the parameter regions for different dynamic phenotypes were computed from Eqs. (13) and (16) , and the coloured parameter regions were identified by the brute-force method. The origin of each panel corresponds to the non-treatment case, and the horizontal and vertical axes correspond to the single-drug treatment cases. Fig. 9 (a) shows the combined effect of antibiotics and emollients. The bifurcation curves dividing the regions for different dynamic phenotypes indicate the minimal amounts of antibiotics ( Ā ) and emollients ( Ē ) required to induce a transition to a more physiological dynamic phenotype. Emollients alone ( Ā = 0 on the left  panel) are not sufficient to induce a transition from bistability to recovery dynamic phenotypes, and antibiotics alone ( Ē = 0 on the middle panel) are also not sufficient for a transition from oscillation to recovery dynamic phenotypes, in the ranges of drug amounts considered. However, a combined use of these drugs can result in the recovery dynamic phenotype. These results suggest that a combined use of antibiotics and emollients is effective, even if a transition to the recovery dynamic phenotypes is difficult to be achieved by single drugs. Within the specified ranges of Ā and Ē , the chronic damage dynamic phenotype (the right panel) is moved only to the bistability dynamic phenotype when (κ P , α I ) = (1 . 6 , 0 . 3) . An induction of the recovery phenotype can be achieved if more potent treatment is given. Alternatively, a healthy state can be achieved by applying transient treatment once the constant treatment has led to the bistability dynamic phenotype (where the two equilibria coexist as in Fig. 4 (c)).
Similar synergistic effects are also observed for the combined use of antibiotics and corticosteroids, as shown by the bifurcation curves given by nonlinear functions of Ā and C ( Fig 9 (b)). The treatment effects are not simply determined by the sum of the amounts of the two drugs but are enhanced by the combined use of two drugs. Particularly, transition from the chronic damage to the recovery dynamic phenotypes is not achieved by the use of single drug but by the combined use of antibiotics and corticosteroids (the right panel). Fig. 9 (c) demonstrates the bifurcation diagrams for the combinatorial use of emollients and corticosteroids, both of which have positive effects on the transitions to the recovery dynamic phenotype from bistability and oscillation phenotypes (the left and middle panels). However, corticosteroids are not effective in moving the patients with the chronic damage dynamic phenotype to a more physiological phenotype (the right panel). The minimal amount of emollients required to induce a transition from the chronic damage to the bistability dynamic phenotypes increases with the amount of corticosteroids. These results suggest that the bifurcation analysis can help finding optimal treatment regimens.
The corresponding results for Case 1 (Fig. S3 ) demonstrate nonbeneficial effects of corticosteroids in terms of transitions of dynamic phenotypes, due to the negative effects of corticosteroids ( β 1 > 0).

A realistic treatment strategy
Finally, we numerically investigate the effects of a more realistic treatment strategy on transition of the dynamic phenotypes. In practice, emollients are recommended to be applied on a daily basis ( Berke et al., 2012 ). When emollient use alone is not sufficient to control AD symptoms, the first choice of additional treatment is the topical application of corticosteroids. To reflect this strategy, we assume continuous emollient treatment, E(t ) = Ē , and intermittent corticosteroids treatment described by Eq. (23) with X = C and T C = 2 . We do not consider antibiotics here because they are mainly used when a secondary infection is evident ( Williams, 2005 ). Fig. 10 (a) shows the averaged level of skin barrier integrity, B ave , after the transient period for different combinations of C and Ē .
The parameter plane of ( C , Ē ) is clearly separated into two regions: one with small values of B ave (the blue region) corresponding to the chronic damage dynamic phenotype ( Fig. 10 (b)), and the other with large values of B ave (the yellow region) corresponding to the recovery dynamic phenotype ( Fig. 10 (c)). The different dynamics for these two parameter regions is characterised by the state of the reversible switch, R ( t ), whose activation corresponds to AD flares. The intermittent treatment with corticosteroids alone ( Ē = 0 ) corresponds to the first region and is not sufficient to induce a transition to the recovery dynamic phenotype, but becomes effective if combined with daily application of emollients.

Discussion
In this paper, we investigated patient-specific effects of AD treatments by conducting bifurcation analysis of the switched nonlinear dynamical model of AD pathogenesis that we previously published ( Fig. 1 ) Domínguez-Hüttinger et al., 2017 ). We used four dynamic phenotypes identified by our previous model analysis ( Fig. 2 )   to describe the severity of AD symptoms, and derived the mathematical conditions that determine when and how the transitions of dynamic phenotypes, or AD severity, are induced by different treatment regimens. The four dynamic phenotypes appear depending on the model parameters that characterise individual patient's genetic, environmental and therapeutic conditions.
Specifically, we derived bifurcation diagrams of the model to identify the parameter regions corresponding to each of the four dynamic phenotypes ( Fig. 4 ). The bifurcation diagram of our model was previously sketched in Domínguez-Hüttinger et al. (2017) , by approximating the switched nonlinear dynamical system by a piecewise affine system and specifying the parameter regions for the four different dynamic phenotypes based on the states of the switches. Here we conducted the bifurcation analysis for the switched nonlinear dynamical systems model, based on the analysis of changes in the stability of the equilibrium and periodic solutions to derive the exact bifurcation diagrams ( Fig. 3 ). We derived the analytical parameter conditions for the bifurcation points and curves, allowing us to systematically evaluate the effects of variations in patient-specific model parameters and treatment options on transitions of the dynamic phenotypes ( Fig. 5 ).
Evaluation of the effects of variations in patient-specific model parameters helps stratification of patients with different genetic or environmental backgrounds. In this study, we evaluated the effects of changes in the pathogen death rate ( δ P ), the environmental pathogen load ( P env ), and the barrier-mediated inhibition of pathogen infiltration ( γ B ) ( Fig. 6 ). Variations in these parameters can be caused by differences in the microbiota composition of the epidermis ( Chieosilapatham et al., 2017;Meisel et al., 2018 ), which is one of the recent therapeutic targets for AD . Providing mathematical tools to assess the benefits of such treatments is of potential clinical relevance. Our method can be easily applied to assess variations in other parameters that may become clinically relevant in the future, as it provides us analytical parameter conditions to evaluate the effects of variations in the model parameters on the bifurcation diagrams.
We further applied our bifurcation analysis approach to systematically evaluate the effects of different AD treatments, for different drug types ( Fig. 7 ), schedules ( Fig. 8 ), and combinatorial use of different drug types with different schedules ( Figs. 9 and 10 ), on transitions of dynamic phenotypes, while considering the individual patient's conditions described by model parameters. Our analysis demonstrated that treatments are not always effective for all patients, but the effectiveness depends on the disease stage described by the original dynamic phenotype of the individual patient, the combination of the drugs applied, and their dosing schedules. On the other hand, similar transition of dynamic phenotypes can be achieved by different treatments. These results are consistent with clinical observations that different treatments are required for different disease stages ( Taïeb et al., 2012 ) and effectiveness of treatments depends on dosing schedule ( Tang et al., 2014 ), suggesting that our mathematical method proposed in this paper can be useful to evaluate multiple treatment options and to design patient-specific optimal treatment strategies, by considering patient-specific information such as genetic background, environmental conditions, the disease stage, and available drug combinations.
A systems biology approach and mathematical methods as proposed here will be especially relevant for evaluation of the compound effects of drug combinations. Prediction of the compound effects is often difficult due to nonlinear interactions leading to non-monotonicity, synergism and antagonism, and is challenging and experimentally expensive due to the combinatorial explosion in the number of the conditions to be explored ( Zimmer et al., 2016 ). We assessed the effects of different treatment strategies with combinations of drugs on patient outcomes, while taking both variations in the dosing and the patient-specificity into account. Finding optimal combination is a big technical challenge that has recently be tackled from a systems biology perspective ( Zimmer et al., 2016 ).
Our method also allowed us to characterise the disease trajectory from the damage (D) to the recovery (R) dynamic phenotypes induced by treatments. The transition may occur directly (DR transition), or indirectly via intermediate asymptomatic-butsusceptible states (DO followed by OR, or DB followed by BR as shown in Fig. 7 (a), (b), (d), (f)). In case of the indirect transitions, it is beneficial for patients to continue treatments to achieve a stable recovery phenotype, even if the disease symptoms disappear to achieve a stable recovery phenotype, rather than staying in the asymptotic-but-susceptible phenotype that is prone to disease aggravation in response to natural and small environmental fluctuations . These results are concordant with recent clinical recommendations to continue treatment for sub-clinical symptoms even after apparent symptoms disappear ( Tang et al., 2014 ).
Our bifurcation-based approach to evaluate the effects of AD treatments is complementary to the optimisation-based approach for treatment scheduling we reported previously . The optimisation-based approach could tailor treatment schedules using heuristic optimisation algorithms under a quantitative objective function. However, it could face with a practical difficulty in collecting full information about patient-specific parameter values. The bifurcation-based approach proposed in this paper could assess the model behaviour systematically for a broad range of the parameter space under simple (e.g. constant) treatment schedules; however, it cannot address complex treatment schedules, because the bifurcation analysis is for the attracting states, and not for transient states. A combination of these two approaches could achieve a comprehensive understanding of the pathogenesis, and help planning effective treatment strategies. In this study, we have applied a brute-force method to investigate the necessary conditions to achieve transitions of dynamic phenotypes by intermittent treatment. As a future work, we can derive the bifurcation point that separates different dynamic phenotypes by applying the technique of the numerical bifurcation analysis for switched dynamical systems ( Tanaka and Aihara, 2013 ). It could handle the time-dependent switch in the amount of drugs applied, as well as the state-dependent switches inherent in the AD model.
In this paper, we formally demonstrated, for the first time, that transitions between AD disease states, described by dynamic phenotypes in our AD model [9], are caused by saddle-node and boundary equilibrium bifurcations of the underlying dynamical system ( Fig. 4 ). This finding suggests that imminent transitions can be detected by the early warning signals associated with these bifurcations (especially in the case of saddle-node bifurcations) ( Boettiger et al., 2013 ). It is known that the system's state would undergo a critical slowing down as the stability of the equilibrium state becomes weakened immediately before it reaches a bifurcation point ( Scheffer et al., 2009 ). The usefulness of early warning signals in predicting disease transitions in advance has been already shown ( Chen et al., 2012;Frey and Suki, 2008 ). Since the saddle-node and boundary equilibrium bifurcation points are relatively close with each other ( Fig. 4 ), the prediction of the saddlenode bifurcation curves by early warning signals enables us to estimate the boundary equilibrium bifurcation curves and thus a rough sketch of a two-parameter bifurcation diagram in a modelfree manner. Clinically, this means that worsening of AD symptoms can be detected early and thus treated preventively, thereby dramatically reducing the cost and side effects of treatments for late and severe disease stages Taïeb et al. (2012) . Therefore, a method to detect early warning signals for the critical transitions (bifurcations) ( Bargaje et al., 2017;Chen et al., 2012;Scheffer et al., 2009 ) could identify a minimal dosage of treatment that can maintain a more physiological dynamic phenotype while avoiding excessive application of the treatments, by deducing the treatment required to cross the bifurcation points from a "spike" in early warning signals, even in the absence of precise knowledge of patient-specific parameters. This could help solving a practically important but challenging issue in our multi-scale modelling approach, i.e. the lack of precise knowledge about the patient-specific parameters, for designing patient-specific optimal treatment strategies from the clinical data measured during transient behaviour.