Computational design of treatment strategies for proactive therapy on atopic dermatitis using optimal control theory

Atopic dermatitis (AD) is a common chronic skin disease characterized by recurrent skin inflammation and a weak skin barrier, and is known to be a precursor to other allergic diseases such as asthma. AD affects up to 25% of children worldwide and the incidence continues to rise. There is still uncertainty about the optimal treatment strategy in terms of choice of treatment, potency, duration and frequency. This study aims to develop a computational method to design optimal treatment strategies for the clinically recommended ‘proactive therapy’ for AD. Proactive therapy aims to prevent recurrent flares once the disease has been brought under initial control. Typically, this is done by using an anti-inflammatory treatment such as a potent topical corticosteroid intensively for a few weeks to ‘get control’, followed by intermittent weekly treatment to suppress subclinical inflammation to ‘keep control’. Using a hybrid mathematical model of AD pathogenesis that we recently proposed, we computationally derived the optimal treatment strategies for individual virtual patient cohorts, by recursively solving optimal control problems using a differential evolution algorithm. Our simulation results suggest that such an approach can inform the design of optimal individualized treatment schedules that include application of topical corticosteroids and emollients, based on the disease status of patients observed on their weekly hospital visits. We demonstrate the potential and the gaps of our approach to be applied to clinical settings. This article is part of the themed issue ‘Mathematical methods in medicine: neuroscience, cardiology and pathology’.


Introduction
Atopic dermatitis (AD) is a chronic inflammatory disease, characterized by recurrent skin inflammation and a defective, permeable skin barrier, that is considered to be caused by complex interactions of genetic and environmental risk factors [1,2]. AD affects up to 25% of children worldwide and has associated socioeconomic burdens [3]. AD is associated with constant itching that may result in chronic sleep loss, and the resultant scratching can cause bleeding and skin infections. The current mainstay of AD treatment is to control the AD symptoms by topical application of corticosteroids or calcineurin inhibitors to reduce the inflammation, in addition to application of emollients to improve the barrier integrity [1,3,4]. However, only 24% of AD patients and carers believe that they adequately manage the symptoms by the current main treatments [5]. This is partly because clear guidance and consensus for effective treatment strategies in terms of frequency, duration and potency are yet to be fully established [2,3,6], while patients are often advised to minimize the use of corticosteroids because of a fear of skin thinning [7].
As a result of the complex underlying mechanisms of AD pathogenesis, AD patients demonstrate a wide spectrum of clinical phenotypes thereby greatly benefitting from personalized treatment [8]. Recently, long-term management of AD focuses on prevention of flares by the so-called 'proactive therapy' [9,10]. Following initial induction of remission to 'get control', proactive therapy aims to 'keep control' by preventing AD flares (inflammation) and achieving skin barrier stabilization. This is achieved by intermittent and scheduled use of lowdose topical corticosteroids or calcineurin inhibitors to the areas of the body that frequently have recurrent flares, even in the absence of the flares. Ideally, effective treatment schedules for proactive therapy can be tailored to each patient, based on the patient's information, such as genetic risk factors, how the symptoms have evolved and the responses to the treatments. This paper investigates the potential of a computational method to inform the design of such personalized effective treatment schedules for proactive therapy, in terms of frequency, duration and potency.
We recently proposed a mathematical model of AD pathogenesis, as an in silico and quantitative framework that coherently explains underlying mechanisms of common AD phenotypes [11]. The model is a system-level representation of the complex and dynamic interplays between immune responses, skin barrier function and environmental triggers that determine the AD pathogenesis; specifically how AD flares start and how AD symptoms exacerbate. Our model simulations reproduced several sets of experimental and clinical results, providing plausible mechanistic and quantitative explanations for dynamic mechanisms behind onset, progression and prevention of AD.
In this study, we extend this experimentally validated mathematical model of AD pathogenesis and propose a new model of treatment effects on AD pathogenesis. We then use the mathematical model to computationally design the personalized optimal treatment schedules for proactive therapy by solving the optimal control problems recursively using a differential evolution (DE) algorithm [12,13], which is an efficient global optimization technique to solve our non-convex optimization problem.

Mathematical model of treatment effects on atopic dermatitis pathogenesis
We consider a mathematical model of treatment effects on AD pathogenesis (figure 1a) obtained by incorporating the dynamic effects of treatments in the previous mathematical model of AD pathogenesis [11]. The proposed model is described by a set of three differential equations: where P(t) is the amount of infiltrated environmental stressors, such as pathogens, that trigger the skin inflammation (AD flare) through activation of innate immune receptors (R(t)), B(t) denotes the strength of skin barrier integrity and D(t) denotes the level of inflammation markers, including pro-inflammatory cytokines (such as TSLP or IL33) and activated dendritic cells. A triple of the state variables, (P(t), B(t), D(t)), represents the patient's disease status described by the levels of infection, barrier defects and inflammation. The variables E(t) and C(t), respectively, represent the potency of emollients and corticosteroids that are applied to achieve skin barrier stabilization and to prevent infection and inflammation. The parameters β 1 , β 2 , β 3 and β 4 represent the relative effects of corticosteroids on the relevant processes (table 1). Other model parameters and their nominal parameter values are taken from [11] (table 2). The proposed model describes the effects of the treatments on AD pathogenesis in a simple form, rather than explicitly incorporating the fine details of the complex molecular and cellular processes. The infiltrated stressors or pathogens, P(t), increase by the penetration of environmental stressors, P env , through the barrier, B(t). P is eradicated by innate immune responses triggered by inflammation (R(t)) and is also naturally degraded. Topical application of corticosteroids (C(t)) provides an anti-inflammatory action [4], resulting in a reduced antimicrobial protein (AMP) expression for eradication of pathogens [14]. The production of B(t) is described by a phenomenological representation of its capacity to self-restore the nominal barrier integrity (B = 1) following its disruption, and is compromised by innate immune responses triggered by inflammation and by cytokines produced by differentiated Th2 cells (where the differentiation is controlled by the master transcription regulator Gata 3, G(t)). Topical application of corticosteroids reduces both inflammation and release of cytokines [3,15], resulting in improved skin barrier function [16]. Application of emollients enhances skin barrier integrity [3,4], helping reduce AD symptoms [17]. The degradation of the barrier occurs as a result of desquamation mediated by active kallikreins, K(t). The level of inflammation markers, D(t), increases while inflammation (R(t)) persists and degrades naturally.
The main aim of the proactive therapy is to prevent AD flares that are assumed to occur as a result of activation of the innate immune receptors. We model the AD flare by a perfect reversible switch between the off-and on-states (R = R off and R on ) with P + and P − denoting, respectively, the activation and inactivation thresholds for the infiltrated stressors (pathogens) that are recognized by the innate immune receptors (figure 1b) [11]. An AD flare occurs when the level of the stressors increases above P + , and stops when it decreases below P − . A similar switching behaviour is assumed for K [11].
In our study, we consider moderate to severe AD patients, who could benefit from proactive therapy for flare prevention and skin barrier stabilization. The severity of the AD symptoms is characterized in our model by two model parameters, κ P and α I , that correspond to the level of the genetic risk factors: mutations in the FLG gene (κ P ) and a dysregulated expression of innate immune system components (α I ) including regulators of antimicrobial peptide expression (e.g. TLRs and NF-κB) [11,18]. A lower α I indicates dysfunctional immune responses and a diminished capacity to eradicate the infiltrating pathogens, and a higher κ P results in a higher skin barrier permeability leading to increased infiltration of stressors that can cause the recurrence of AD flares. To simulate these moderate to severe AD patients, we pick the values of (α I , κ P ) from the parameter region for which a unique pathological steady state exists that corresponds to the high level of stressors leading to a persistent AD flare and completely damaged skin barrier. We also assume R off = 1 and K off = 1 to model the effects of subclinical inflammation in non-lesional skin of severe AD patients [10], and G(t) = 1 corresponding to systemic Th2 sensitization due to the persistent AD flare [11].

Optimal control problem formulation
Using the proposed model, we computationally design optimal treatment strategies for proactive therapy with a combination of emollients and corticosteroids. The proactive therapy consists of two phases: an 'induction of remission' phase where we aim to suppress the clinical inflammation, followed by a 'maintenance of remission' phase where we apply intermittent but scheduled treatment to prevent the recurrence of the AD flare [9,10]. We refer to the two phases as 'induction phase' and 'maintenance phase' hereafter. The induction phase aims to 'get control', and the maintenance phase aims to 'keep control' [10]. To comply with the current clinical recommendations [3,19], we assume that emollients are applied constantly throughout both phases at a low level,Ē, which is insufficient to achieve the remission by itself for moderate to severe AD patients [20]. We, therefore, design the optimal schedules for application of corticosteroids that can induce and maintain the remission. We consider the on-off treatment at discrete times, reflecting the daily application or non-application of corticosteroids with different potencies (figure 2a).
We formulate the problem as finding optimal treatment strategies that minimize the duration and potency of the treatments that effectively move the state variables from the initial steady state to the specified target state, while the dynamics of the state variables is determined by our model (figure 2b). The target state is defined by a target level of P(t) that does not lead to recurrence of an AD flare (figure 1b), as the proactive therapy aims to prevent the AD flare.  Figure 2. Optimal control problem formulation to design treatment strategies for proactive therapy for AD. The whole period consists of the induction phase (red) and maintenance phase (blue). (a) An example dynamics of a state variable (P(t)) with onoff treatment of corticosteroids (C(t)) of different potencies at discrete times, in addition to constant application of emollients (E(t)) of a fixed potencyĒ. (b) An example trajectory of a projection of a patient's states, (P, D), when the optimal treatment strategy is applied. The states move from the initial state ((P 0 , D 0 ), red circle) towards the target level,P r , of the induction phase (red vertical line), and then to the target level,P m , of the maintenance phase (blue vertical line). The state B is omitted in the figure. (c) The recursive optimal control problem to be solved. The optimal treatment strategy for each cycle (either for the duration of T r or T m ) is determined by predicting the optimal evolution of the state variables (dashed lines) based on the measurement of the states at the beginning of the period (circles). The actual evolution of the state (solid lines) can be different from the prediction, for example if the calculated optimal treatment strategy is not applied.
The objective function, J, is described by J = k 1 J 1 + k 2 J 2 + k 3 J 3 + k 4 J 4 , where J 1 is the penalty for the treatment duration, corresponding to the patients' burden to apply the treatments, J 2 is the penalty for the total amount of the treatment applied (duration × strength), representing the financial cost as well as the risk of side effects due to the excessive use of corticosteroids, J 3 is the penalty for the final state to be deviated from the target state and J 4 is the penalty for the trajectory to be deviated from the target state. The functions J 1 , J 2 , J 3 and J 4 are defined for the induction and maintenance phases as shown below. We assume k 3 = k 4 (the coefficients for the penalty on the deviation from the target state) for simplicity.

(a) Induction of remission phase
In the induction phase, we assume that corticosteroids are constantly applied during the calculated optimal duration, T r . The optimal problem to be solved is formulated so as to find a pair of values, (T r ,C), where T r is the duration of the induction phase andC is the potency of corticosteroid to be constantly applied during this phase to minimize the objective function under the constraints 0 ≤ T r ≤ T max r and 0 ≤C ≤ C max . We set the target level,P r , to be smaller than P − where the AD flare ceases. The objective function for the induction phase consists of J 1 = (T r /T max r ) 2 , J 2 = (T r /T max r )(C/C max ) 2 , J 3 = Φ r (P(T r )) + (1 − P(T r )/P r ) 2 and ) is a non-convex function that penalizes the failure to achieve remission and takes the values of 100

(b) Maintenance of remission phase
Once the remission is induced, the proactive therapy proceeds to the maintenance phase, where corticosteroids are applied intermittently to prevent the recurrence of AD flares. We design the optimal treatment schedule for a week, based on the measurement of the state variables at the beginning of the week, and repeat the weekly cycle (figure 2c). This scenario corresponds to weekly hospital visits of the patients where clinicians evaluate the disease state and plan the treatment strategy until the next visit.
For the ith cycle of the maintenance phase, we calculate the optimal treatment strategy (T i C ,C i ) that minimizes the objective function, under the constraints on the duration of corticosteroid application, 0 ≤ T i C ≤ T m = 7 days, and its potency 0 ≤C i ≤ C i max . For the whole duration of the maintenance phase, we set the target level,P m , that is smaller than P + to avoid the recurrence of the AD flare (figure 1b). In the next section, we investigate the effects of the choice ofP m on the calculated optimal treatment strategies. The objective function for the ith maintenance cycle consists of , represents the penalty on the re-occurrence of an AD flare during the maintenance phase, and takes the values of 100 + 0.1(P(T r + iT m ) − P − ) if both P(T r + iT m ) > P − and R(t) = R on are satisfied, and 0 otherwise.

Computational identification of optimal treatment schedules for moderate to severe atopic dermatitis patients
We used DE to solve the optimal control problem formulated above, for different scenarios, to test the applicability of our approach. The nominal values used for the simulations for moderate to severe AD patient cohorts, such as the target levels and the constraints on the strength and duration of treatments, are summarized in table 1. We first confirmed that the optimal treatment strategy calculated using the nominal parameter values suggests a length of the induction phase, T r , that is clinically relevant (less than or equal to four weeks [21]). Indeed, the optimal induction period was calculated to be T r = 19 days, while the AD flare stopped within the first 3 days (figure 3a). We also confirmed by global sensitivity analysis that this calculated optimal strategy is robust to changes of the model parameters and choice of the weighting coefficients for the objective function (electronic supplementary material, figures S1 and S2).

(a) Effects of the choice of the maintenance target level
We investigated the effects of the choice of the target level during the maintenance phase on the calculated optimal treatment strategies ( figure 3).
In the nominal case with the nominal parameter values, we set the maintenance target level to beP m = 26, which is lower than the deactivation threshold (P − = 26.6) but higher than the induction target level (P r = 24). The calculated optimal maintenance treatment (figure 3a) demonstrates that intermittent application of corticosteroids by 3 days per week could achieve the maintenance without recurrence of the AD flare for the whole duration of the maintenance phase investigated (eight weeks).
When the maintenance target is chosen to be closer to the activation threshold (P + = 40), for exampleP m = 30, the calculated optimal treatment strategy suggests the application of corticosteroids by 1 day per week to maintain the remission (figure 3b). This higher target level is much easier to be maintained with a smaller amount of corticosteroids. However, the resulting state corresponds to a lower barrier integrity with a higher level of infiltrated stressors, compared with the nominal scenario (figure 3a), due to a smaller amount of corticosteroid applied. This worsening of the state may make patients more vulnerable to an increased level of environmental stressors due to random or natural fluctuations that can retrigger an AD flare.  is set to be the same (P r = 24) for all the scenarios. The dynamics during the induction phase and that during the maintenance phase are shown in red and blue, respectively. Remission is achieved when the level of the stressors is decreased below P − , and the AD flare reoccurs when it increases above P + . The induction phase continues even after the AD flare ceases, and tries to bring the state towards the target level.
On the contrary, when the maintenance target is chosen to be further lower than P − , for exampleP m = 20, the optimal treatment strategy requires an increased amount of corticosteroids, with application of more potent corticosteroids for 6 days per week, to achieve a very low target value (figure 3c). While the calculated optimal treatment strategy can successfully prevent recurrence of inflammation, this strategy is not desired due to the high amount of corticosteroids applied in total (2.7-fold increase from the nominal case).
These results suggest the importance of the choice of the maintenance target level,P m , as a design criterion for optimal treatment strategies. We decided to useP m = 26 for further simulations, as it ensures the successful maintenance of remission without the need of excessive treatment amount during the maintenance phase and it may correspond to the clinically suggested so-called weekend therapy.
We also investigated the effects of choice of the induction target,P r (figure not shown). Decrease ofP r from our nominal value ofP r = 24 resulted in an optimal strategy that requires an increased amount of corticosteroid (in both the potency and application time) during the induction phase, but did not affect the strategy during the maintenance phase.

(b) Stratification of patient cohorts
To demonstrate that our approach is applicable to different patient cohorts, we computationally obtained the optimal treatment strategies for different values of two model parameters, (κ P , α I ), that can specify patient cohorts by their strengths of the two main genetic risk factors, mutations     in the FLG gene (κ P ) and a dysregulated expression of innate immune system components (α I ).
For the case with (κ P , α I ) = (0.85, 0.04) (figure 4a), the calculated optimal treatment strategy was successful in achieving remission and preventing AD flares for eight consecutive weeks, by an induction treatment of 20 days followed by 3 days per week maintenance treatment. However, it requires a higher amount of corticosteroids (48% increase for the induction phase and 3.3% increase for the maintenance phase) of a higher potency, compared to the nominal case (figure 3a), to combat the higher initial stressor load. When the risk factor of dysfunctional immune responses becomes even stronger, as for the case with (κ P , α I ) = (0.85, 0.03) (figure 4b), the calculated optimal treatment schedule failed to achieve remission, leading to sustained or unresolved AD symptoms with a low barrier integrity. Indeed, the 'maintenance' therapy we computationally calculated after 14 days of the 'failed' induction of remission suggests a continuous use of corticosteroids, meaning that these virtual patients will require constant, rather than intermittent, application of corticosteroids (figure 4b).
For the cohorts with an increased barrier permeability (κ P = 0.9), an increased amount of corticosteroid (by 54.7%) is required during the maintenance phase (figure 4c), compared to the nominal cohort with the same α I = 0.05. Further increase in the barrier permeability (κ P = 0.95) led to failure in inducing remission and the optimal strategy suggests to continue the constant application of corticosteroids (figure 4d). Synergistic effects of the two risk factors, by increasing κ P and decreasing α I simultaneously ((κ P , α I ) = (0.9, 0.04), figure not shown), resulted in an optimal strategy with a 5.8-fold increase in the total amount of corticosteroid (2.35-fold in the induction   The dynamics during the induction phase starting from the initial state (red circles) is shown in red, and that for the maintenance phase is in blue. phase and 11.3-fold in the maintenance phase), compared to the nominal cohort. However, the optimal strategy still could not achieve remission.
These results suggest that our approach is applicable to different virtual patient cohorts, and that it could help stratify the virtual patients into those who would benefit from the calculated optimal treatment strategies, and those who would require additional or stronger treatments, such as systemic treatment, to achieve remission. Our results also demonstrate how the treatment efforts scale with the level of the two common AD risk factors. It will be interesting to compare the computational predictions with the patients' data that relate the severity of the patients' symptoms to the required treatments and the actual treatments prescribed. Evaluation of how the treatment efforts scale with the initial disease severity before the start of treatments in the clinic is an interesting future research topic.
(c) Effects of poor adherence to suggested optimal treatment schedule So far, we assumed that the patients follow the calculated optimal treatment schedule. However, AD patients do not necessarily always follow the treatment guidelines. This problem of poor adherence to treatment could have negative effects on long-term treatment of AD [22,23]. Using   our proposed approach, we investigated the effects of poor adherence to the calculated optimal treatment schedule, particularly how the future optimal treatment schedule and the evolution of the disease state are affected. Consider the case where the virtual patients do not complete the calculated optimal induction phase and stop using corticosteroids after 5 days when the AD flare disappears. This scenario can occur, for example, due to corticosteroid phobia. If they continue their daily emollient treatment (figure 5a), the optimal strategy in the subsequent weeks suggests the use of an increased amount of corticosteroid, with a 1.4-fold more potent corticosteroid for a longer period (by 1 day) during the first maintenance cycle, compared with the nominal case. The effect of non-adherence is not predicted to be dramatic, possibly because the AD flare ceased already within the first 5 days of corticosteroid use. However, if the virtual patients also stop the daily emollient treatment (figure 5b), we observe the recurrence of an AD flare with a severe worsening of the symptoms (a sharp decrease in B and a dramatic increase in P and D). As a result, the calculated optimal strategy in the subsequent cycle suggests the use of a more potent corticosteroid (by 1.4-fold) than that applied during the induction phase, for 5 days. This corresponds to a prolonged induction phase with a more potent corticosteroid. If we prolong the maintenance cycle to 20 days from 7 days, the calculated optimal strategy suggests the application of a corticosteroid of 10% more potent than that used for the induction phase, for an extra duration of 10 days (figure not shown). These computational predictions suggest the benefits of continual use of emollients in reducing AD symptoms, as shown in [3].
Another scenario to be considered is the non-adherence during the maintenance phase, following the successful completion of the induction phase. When the virtual patients stop their daily application of emollients for the first three maintenance cycles (figure 5c), it results in an immediate worsening of the symptoms (a decrease in B). The calculated optimal strategy suggests application of more potent corticosteroids (by 1.4-fold) than the nominal case for 3 days per week in the subsequent three weeks to maintain the remission. When they also stop applying corticosteroids during these three weeks (figure 5d), a severe worsening of the symptoms (a sharp   Figure 5. Effects of poor adherence to the calculated optimal treatment strategies. Non-adherence to the suggested corticosteroid treatment after 5 days in the induction phase, with continual use of emollients (a) and without use of emollients (b). Non-adherence to the suggested corticosteroid treatment during the maintenance phase, without emollient application for the first three cycles (c), and without both corticosteroid and emollient application for the first three cycles (d) and four cycles (e). Green lines represent the period during the non-adherence, and the dotted lines on the right column demonstrate the optimal strategy calculated for the nominal case where patients adhere to the suggested optimal strategy. decrease in B and a dramatic increase in P and D) is observed. The calculated optimal strategy then suggests application of an increased amount of corticosteroid (1.5-fold more potent for two more days) than the nominal case to maintain the remission in the subsequent weeks. If they miss the treatments for four weeks (figure 5e), the AD flare reoccurs and the optimal strategy suggests the use of a much more potent corticosteroid (by 1.35-fold) than that used during the induction phase, for the duration of 5 days.

Discussion
In this paper, we proposed a computational method to inform the design of patient-specific optimal treatment strategies for moderate to severe AD patients who require constant treatment for stabilization of their AD symptoms. Our proposed framework solves optimal control problems recursively to design treatment schedules for proactive therapy that aims to prevent AD flares and to achieve skin barrier stabilization. The proactive therapy consists of intermittent and scheduled use of low-dose corticosteroids, in addition to a constant application of emollients, once initial induction of remission has been achieved. The objective functions to be minimized correspond to the penalties on the duration and the potency of the treatments applied, as well as the deviation from the target states we specify.
One of the main difficulties in formulating an optimal control problem is to identify appropriate objective functions. Here, we systematically explored different possible target levels for the clinically relevant variables to be controlled (P) and found the most adequate maintenance target level (figure 3), which we used to derive robust treatment strategies even with poor adherence (figure 5), for our nominal patient cohorts and for patient cohorts with severer genetic risk factors (figure 4). These results suggest the importance of choosing appropriate target states to successfully maintain remission, and that our proposed mathematical framework can be used to investigate the effects of poor adherence to the optimal treatment strategies systematically. We could also identify those virtual patient cohorts that would require stronger treatments, with even higher doses or additional pharmacological substances such as antibiotics, phototherapy or systemic immunosuppressant treatment to achieve adequate disease control (figure 4). We evaluated the sensitivity of our results to the changes in the values of k 1 , k 2 , k 3 and k 4 , and the model parameters, and confirmed that the optimal strategies calculated for the nominal case are robust (electronic supplementary material, figures S1 and S2).
Our systematic and computational approach could be effective in informing the design of personalized optimal treatment strategies, and help to solve the issues with the current lack of a clear guidance and consensus for effective treatment strategies, and to minimize the potential side effects of long-term use of corticosteroids. Moderate to severe AD patients require repeated treatment to stabilize their pathological state that would naturally remain unresolved if treatment is not applied. In addition, moderate to severe AD patients usually require a combination of treatments (such as corticosteroids and emollients) to be applied to achieve successful maintenance of remission in two phases (induction and maintenance phases). Accordingly, designing the optimal treatment strategies to stabilize the AD symptoms for these patients may benefit from an advanced optimization technique such as the one proposed in this paper.
In this paper, we computationally obtained the optimal treatment strategy by recursively solving the optimal control problems. The proposed computational framework could be easily extended to the application of model predictive control (MPC), which uses the measured states of the system to predict and optimize the control input that minimizes the objective function over a future time horizon [24]. MPC has been already successfully applied to design treatment profiles for diabetes [25,26], prostate cancer [27] and leukaemia [28,29]. Application of MPC, i.e. inclusion of the receding horizon, will allow us to obtain smoother graded therapy, as it will enable us to find the optimal treatment schedule that does not necessarily achieve the target level within each maintenance cycle but achieves it gradually over a longer period.
Our simulation results in this paper show that poor adherence to the suggested optimal treatment schedule inevitably leads to higher treatment doses in subsequent cycles. This indicates a potential benefit of using our approach under more realistic clinical scenarios to provide a theoretical argument to recommend patients to adhere to the suggested treatment strategies. Our results are also consistent with the current clinical recommendations, for example, the weekend therapy with application of corticosteroids for two consecutive days per week, in addition to the daily application of emollients. Our investigation on the effects of the AD flare that occurs during maintenance therapy demonstrated that resuming continuous use of a much stronger corticosteroid can successfully achieve remission but result in an increase in the total amount of corticosteroids applied. As an important next step, we need to compare the computationally obtained treatment strategies with treatment options that are currently used in the clinical setting.
We also need to develop robust ways to identify the simulation parameters and the model parameters from each patient's clinical data such as initial skin thickness and pattern of eczema, in order to effectively calculate the optimal treatment strategies. If the temporal data of patients become available, the information on the discrepancy between the measured values and the model prediction will be used to identify the model parameters online. The approach we proposed in this paper could be a first step towards designing personalized effective treatment strategies for prevention, and adequate control, of AD symptoms. Exploring the personalized optimal treatment strategies in the clinical setting would be challenging if we do not apply a systematic and computational approach, because of the combinatorial explosion of treatment types, durations, potencies and each patient's information. As the model was developed based on the understanding of the pathological mechanisms, the obtained treatment strategies could be readily interpreted and the framework is applicable to different patient cohorts and different scenarios. For example, it will help to identify a way to reduce the frequency of clinic visits by placing control back in the hands of parents and children, evaluate the effects of reduced visits and whether we can still achieve the optimal treatment strategies by monthly or bimonthly clinical visits.
This paper demonstrated the proof of concept of the computational design of optimal treatment strategies, using a mathematical model that describes the treatment effects in a simple form. We will further investigate the appropriateness of the model description of the treatment effects, using dynamic data of AD patients after application of corticosteroids and emollients.

Material and methods
All simulations were conducted using Matlab v. R2016a (The MathWorks, Inc., Natick, MA, USA). We used ode45 for the numerical integration of the system and events location functionality of Matlab to identify the switching boundaries of the hybrid system. We identified the optimal treatment strategy for the induction phase by applying a DE algorithm with 1000 generations. Starting from 30 randomly chosen initial vectors, (C, T r ), with 0 ≤C ≤C max and 0 ≤ T r ≤ T max r , we found the optimal solution by evolving a population of 30 vectors at each generation, using the mutation strategy DE/rand-to-best/1 [30] with a differential weight of 0.6 and recombination with a crossover rate of 0.5. The same procedures were repeated for each maintenance cycle. The global sensitivity analysis was conducted by simultaneously varying the values of the model parameters or the weights for the objective functions from their nominal values by ±50% for 529 and 400 simulations, respectively.