A Dynamic Model for Prediction of Psoriasis Management by Blue Light Irradiation

Clinical investigations prove that blue light irradiation reduces the severity of psoriasis vulgaris. Nevertheless, the mechanisms involved in the management of this condition remain poorly defined. Despite the encouraging results of the clinical studies, no clear guidelines are specified in the literature for the irradiation scheme regime of blue light-based therapy for psoriasis. We investigated the underlying mechanism of blue light irradiation of psoriatic skin, and tested the hypothesis that regulation of proliferation is a key process. We implemented a mechanistic model of cellular epidermal dynamics to analyze whether a temporary decrease of keratinocytes hyper-proliferation can explain the outcome of phototherapy with blue light. Our results suggest that the main effect of blue light on keratinocytes impacts the proliferative cells. They show that the decrease in the keratinocytes proliferative capacity is sufficient to induce a transient decrease in the severity of psoriasis. To study the impact of the therapeutic regime on the efficacy of psoriasis treatment, we performed simulations for different combinations of the treatment parameters, i.e., length of treatment, fluence (also referred to as dose), and intensity. These simulations indicate that high efficacy is achieved by regimes with long duration and high fluence levels, regardless of the chosen intensity. Our modeling approach constitutes a framework for testing diverse hypotheses on the underlying mechanism of blue light-based phototherapy, and for designing effective strategies for the treatment of psoriasis.


INTRODUCTION
Blue light (BL) at 453 nm is non-toxic (Awakowicz et al., 2009;Kleinpenning et al., 2010) and decreases the symptoms of psoriasis vulgaris (Pv) (Weinstabl et al., 2011;Pfaff et al., 2015), a chronic, inflammatory skin condition that affects 2-3% of the world's population (Parisi et al., 2013). Psoriasis vulgaris is characterized by hyper proliferation and lowered differentiation of skin keratinocytes (Weinstein et al., 1985), evident in lesional areas of thick skin (Perera et al., 2012). These areas also exhibit sustained inflammation caused by immune cells such as T cells and dendritic cells (Perera et al., 2012). Blue light minimizes the proliferation of keratinocytes and induces their differentiation in a wavelength and fluence dependent manner (Liebmann et al., 2010). Additionally, blue light irradiation suppresses dendritic cell activation (Fischer et al., 2013). These effects of blue light on keratinocytes and immune cells may explain the reduced inflammation and diminished epidermal thickness of lesional psoriatic skin after the treatment. Nevertheless, the underlying mechanism of this therapeutic approach is not fully understood. It is not clear how the cellular processes of proliferation and differentiation are modified in the cells after irradiation with blue light. Further, it is uncertain whether blue light affects only proliferative cells or both proliferative and non-proliferative cells. It is hypothesized that BL improves psoriatic skin by decreasing the proliferative capacity of keratinocytes.
Large variations are observed in the results from the available clinical investigations (Maari et al., 2003;Weinstabl et al., 2011;Kleinpenning et al., 2012;Pfaff et al., 2015). These discrepancies in the reported outcome may be due to differences in the main treatment parameters, i.e. length of treatment (days), fluence F (Jcm −2 ) also referred to as dose, and power density (mWcm −2 ) also denoted as intensity. Additionally, a BL treatment session may occur one or more times per week for a certain number of weeks (Maari et al., 2003;Weinstabl et al., 2011). Despite all the potential combinations of treatment parameters and their impact on the effectiveness of the treatment, the effect of variations in the irradiation parameters of BL on the dynamics of keratinocytes is unknown. Additionally, no clear guidelines have been defined in the literature for BL treatment of psoriasis. The current treatments for psoriasis include phototherapy with ultraviolet (UV) light (Schneider et al., 2008). UV phototherapy is effective but only indicated for severely affected patients due to the risk of skin cancer (Pathak, 1991) caused by DNA damage. In contrast with UV, blue light is not toxic to skin cells (Awakowicz et al., 2009). The difference in the therapeutic effects of both spectral ranges is the result of the interactions between these spectral ranges and their molecular photoacceptors. UV is absorbed by DNA (Markovitsi, 2016) and yields DNA damage (Lee et al., 2013), cell death, and remission of psoriasis. Blue light is accepted by flavins (Eichler et al., 2005) and porphyrins (Dai et al., 2012), resulting in a decreased proliferative capacity and management of psoriasis. The non-toxicity and beneficial effects of blue light make it an attractive alternative to UV phototherapy. Emphasizing the need for establishing clear treatment recommendations that lead to an effective therapeutic regime using blue light.
Considering that the details of the mechanism of blue lightbased therapy for psoriasis are still not fully elucidated, it seems suitable to use computational methods for their analysis. Computational models have previously been used to predict cellular behaviors (Savill, 2003;Gandolfi et al., 2011) and the effects of ultraviolet irradiation on the skin (Weatherhead et al., 2011;Zhang et al., 2015). This approach may contribute to the investigation of processes occurring in the epidermis after blue light irradiation, and extend our knowledge on the principles of regulatory effects induced by this spectral range on the keratinocytes dynamics.
Here, we study the blue light treatment of psoriasis using an in silico approach. We first used the model to explore whether a temporary decrease of keratinocytes hyper-proliferation can explain the outcome of phototherapy with blue light. The model accurately described the response to blue light therapy. The simulations suggested that the observed decrease in the keratinocytes proliferation rate is sufficient to reduce the epidermal thickness and severity of psoriasis. However, it was not sufficient to allow psoriatic epidermis to completely remodel back to a healthy phenotype, regardless the treatment scheme. Then, we analyzed the effect of length of treatment, fluence, and intensity on the management of this inflammatory skin condition. The model predicted that high efficacy is achieved by treatment schemes with long duration and high fluence levels.

MATERIALS AND METHODS
To describe the management of psoriasis by blue light, we implemented a computational model for BL irradiation of psoriatic skin (BLISS; Figure 1). This model is defined by a set of 12 ordinary differential equations (ODEs) describing the time evolution of keratinocytes as they move vertically through the layers of the epidermis while blue light is shined upon them. BLISS was developed based on the phenomenological observations of decreased proliferation and increased differentiation of keratinocytes due to blue light, particularly at a wavelength of 453 nm. In this section we describe the general structure of the model, its implementation, and analysis.

Computational Model
In contrast with UV, BL does not lead to cell death below fluences of 500 Jcm −2 (Awakowicz et al., 2009). Instead, it affects the proliferation and differentiation of the keratinocytes, which are key processes in the model we propose. The general structure and assumptions considered in the model regarding the kinetics of keratinocytes in psoriasis are based on the work of Zhang et al. (2015) for UV phototherapy of psoriasis. However, the main difference with their model is the implementation of the underlying mechanism by which the kinetics of the keratinocytes are affected upon BL irradiation.
In the model (Figure 1), we represent psoriasis as a bi-stable system and assume that keratinocytes show either a healthy or diseased phenotype. We consider that both populations coexist and interact within the epidermis. This assumption allows the description of the effects induced by blue light therapy. The set of ODEs describing the interactions between healthy and diseased keratinocytes under the influence of blue light are presented by Equations (1-12). This set of equations describes the kinetics of keratinocytes at six stages of differentiation in both their healthy (Equations 1-6) and diseased (Equations 7-12) state, i.e., stem cells (P sch,d ), transit amplifying cells (P ah,d ), growth arrested cells (P gah,d ), spinous cells (P sph,d ), granular cells (P gch,d ), corneocytes (P cch,d ). Stem cells and transit amplifying cells form the proliferative compartment of the epidermis, and corneocytes are the end point of the differentiation process. The ODEs account for the cellular processes of proliferation (γ h,d ), differentiation (k h,d ), apoptosis (β h,d ), and desquamation (α h,d ). The cells in the proliferative compartment may divide in one of three modes, i.e., produce to daughter cells equal to the progenitor [self-proliferation (γ h,d )], generate two daughter cells where one corresponds to the next stage of differentiation [asymmetric division (k ah,d )], or induce two daughter cells where both correspond to the next stage of differentiation [symmetric division (k sh,d )]. The description of each parameter considered in these equations is presented in Table 1. These parameters were derived from one of three sources, i.e., the literature, calculated from other model parameters as specified in Table 1, or estimated by fitting the model to experimental data of Liebmann et al. (2010).  Clayton et al., 2007;Zhang et al., 2015 k 1ahom Minimal asymmetric healthy stem cell division rate constant 1.31 × 10 −2 d −1 Clayton et al., 2007;Zhang et al., 2015 γ 2h Healthy transit amplifying cells self-proliferation rate constant 1.40 × 10 −2 d −1 Bauer et al., 2001;Hoath and Leahy, 2003;Zhang et al., 2015 k 2sh Healthy transit amplifying cells symmetric division rate constant 1.73 × 10 −2 d −1 Bauer et al., 2001;Hoath and Leahy, 2003;Zhang et al., 2015 k 2ah Healthy transit amplifying cells asymmetric division rate constant 1.38 × 10 −1 d −1 Bauer et al., 2001;Hoath and Leahy, 2003;Zhang et al., 2015 k 3h Healthy growth arrested to spinous cells differentiation rate constant 2.16 × 10 −1 d −1 Bauer et al., 2001;Hoath and Leahy, 2003;Zhang et al., 2015 k 4h Healthy spinous to granular cells differentiation rate constant 5.56 × 10 −2 d −1 Bauer et al., 2001;Hoath and Leahy, 2003;Zhang et al., 2015 k 5h Healthy granular cells to corneocytes differentiation rate constant 1.11 × 10 −1 d −1 Bauer et al., 2001;Hoath and Leahy, 2003;Zhang et al., 2015 k In psoriasis, the increased number of keratinocytes is the result of a hyper-proliferative population of basal cells and a sustained activation of the immune system. In the model, the hyperproliferative population is represented by the diseased stem cells and transit amplifying cells proliferating at a rate γ d , which is faster than the proliferative rate of their healthy counterparts. It is assumed that there are a maximum number of stem cells P max sc available in the epidermis, which limits the numbers of healthy and diseased stem cells P sch,d. The diseased stem cells have a larger growth capacity defined by fold increase λ. The sustained activation of the immune system is considered through the removal of diseased stem cells in Equation (7). The sustained immune system response is regulated by the density of diseased stem cells and defined by the maximum killing rate K p and the half activation of the immune system K a due to psoriatic stem cells. The immune response is ample when the cell density of diseased cells exceeds the threshold defined by K a .
In Equation (13), γ 1,hom, k 1a,hom, k 1s,hom are the homeostatic rate constants for each division process, and P ta,hom is the homeostatic density of transit amplifying cells (Zhang et al., 2015). The density of the total transit amplifying cells population is the sum of the healthy and diseased groups of transit amplifying cells. The maximum increase in the growth fraction and decrease in the cell cycle time of fast proliferating stem cells is indicated by ω. n represents the steepness of the stem cells proliferation rate regulated by the transit amplifying cells population. When the number of healthy transit amplifying cells is low, the proliferation and division rates of the healthy stem cells increase. Alternatively, the healthy stem cell parameter rates are at their minimum values when the population of healthy transit amplifying cells equals its homeostatic population. It is assumed that psoriatic stem cells are not regulated by the population of transit amplifying cells. Thus, they proliferate with rates ρ sc fold higher than the healthy stem cells rate constants in homeostasis. The proliferation and differentiation rates of diseased transit amplifying cells (γ 2d , k 2sd , and k 2ad ) are ρ ta fold changes higher than the rates of healthy transit amplifying cells. Similarly the diseased cells in the non-proliferative compartment, i.e., P gad and P spd , differentiate ρ tr times faster than their healthy counterparts. The desquamation (α d ) of psoriatic corneocytes is affected by ρ de . Note that Equation (11) equals to zero, given that granular cells are lost in psoriasis due to abnormal differentiation.
Irradiation with blue light exponentially decreases the proliferation rate of all healthy and diseased proliferative keratinocytes by a fluence dependent factor θ BLγ with a value between 0 and 1 (Equation 14).
with a γ > 0 and b γ < 0 (14) Based on the keratinocytes behavior experimentally observed (Liebmann et al., 2010), blue light might also directly affect the differentiation rate of the keratinocytes in an exponential manner by a fluence dependent factor θ BLk with a value between 0 and 1 ( Equation 15). Parameters related to the decrease of the Frontiers in Physiology | www.frontiersin.org proliferation and increase of differentiation rates were initially estimated by fitting in vitro data from Liebmann et al. (2010).
with a k , b k and c k > 0 (15) Blue light phototherapy is divided into sessions with a certain fluence F, defined by the product of the irradiation time and the average intensity I av at which light is shined on the skin. The treatment sessions can occur one time per day in a daily or spread out basis. In the model we consider that each day of treatment affects the proliferation of the keratinocytes, based on F >0. In days where no treatment is provided F = 0 and θ BLγ = θ BLk = 1. The apoptosis rates of the healthy keratinocytes β h (Equation 16) are determined by the differentiation rate of each population (k 1sh , k 2sh , k 3h , k 4h , k 5h ) and the probability of a cell undergoing apoptosis (apoptotic index AI h ). The apoptotic index is a common measurement for quantifying the extent of apoptosis for a given cell population. It is defined as the ratio of apoptosis rate to the total out flux. At fluences higher than 500 Jcm −2 BL increases the apoptosis of keratinocytes by factor θ BLβ . Thus, accounting for the cytotoxicity observed only at high fluences of blue light (Awakowicz et al., 2009). Below 500 Jcm −2 , θ BLβ has a value of 0, for F between 500 and 750 Jcm −2 its value is 0.039, and 0.05 above 750 Jcm −2 .
The diseased keratinocytes have an apoptosis rate defined by the apoptotic index of psoriatic cells AI d , the differentiation rate of each population (k 1sd , k 2sd , k 3d , k 4d , k 5d ) and q BLβ . (Equation 17).
The skin has specific wavelength dependent optical properties that must be considered in the model, i.e., refractive index, absorption coefficient, scattering coefficient, and anisotropy factor (van Gemert et al., 1989). These properties vary per skin layer and skin type (Zonios et al., 2001). BLISS accounts for a low perfused type 1-3 skin with an epidermal energy absorbance level of 57.9%. This value was computed using an unpublished fivelayered optical model developed in LightTools (Synopsis). The energy absorbance value defines the amount of BL energy that is absorbed by the epidermis. Skin types 4-6 have a higher energy absorbance, however are not considered in the model.

Model Implementation and Analysis
The model was implemented in Matlab (The Mathworks Inc.).
The ODE system was solved with ode-solver ode15s. The inputs of the model were the BL irradiation parameters ( Table 2) and the initial cell density of the 12 keratinocyte populations (Table 3) representing an psoriatic epidermis (Simonart et al., 2010;Zhang et al., 2015). The total keratinocytes cell density in the simulated epidermis at t = 0 was of 217,556 cells per mm 2 , which is 2 times higher than that of healthy skin (∼100,000 cells per mm 2 ; Hoath and Leahy, 2003). BLISS is available at GitHub 1 and 1 https://github.com/ZFelixGarza/BLISS  the BioModels Database 2 (Chelliah et al., 2015) with identifier MODEL1701090001.
A full treatment of blue light was simulated by modifying the proliferation, differentiation, and apoptotic rates on the day of a BL session during the whole treatment. For example, if the treatment sessions occurred three times per week for 8 weeks, the rates were only varying on those specific days during the 8 weeks of treatment. During a treatment session, light can be shined in a continuous (CW) or pulsed (PW) mode. In CW the fluence of blue light, defined as power density I multiplied by time t, is provided at a constant peak power density I p , which is equal to the average power density I av . In PW, the epidermis is irradiated in short pulses at an I p higher than I av . The model allows for the simulation of both modes based on a given duty cycle, fluence, and treatment time. These inputs are used to compute the I p and duration of each peak and define the instants in which the proliferation rates are decreased during the treatment time.
BLISS calculates the cell densities over time, during and after blue light phototherapy for each healthy and diseased keratinocyte population. The model's results can be directly compared to in vitro data. However, in psoriasis most clinical data is given in terms of the psoriasis area and severity index (PASI), or its local form, i.e., the local psoriasis severity index (LPSI; Maari et al., 2003;Weinstabl et al., 2011;Pfaff et al., 2015). PASI is a quantitative rating score for measuring the severity of a psoriatic lesion based on the degree of erythema, scaling and induration per anatomic area, i.e., head, trunk, upper, and lower limbs. Each of the three characteristics is graded on a scale from 0 to 4, giving a maximum score of 12 per area based on the percentage of coverage. The maximum PASI is 72, mild psoriasis which is the target of blue light treatment corresponds to a PASI of 0-10. Induration is related to the thickness of the epidermis, and scaliness is related to the status of the Stratum Corneum.
Both features are comprised within the model. Thus, assuming an initial LPSI value and a positive correlation between thickness and scaliness with the cell density of keratinocytes, the LPSI is derived from the relative change in the total cell density of keratinocytes as indicated by Equation (18). A similar approach has been previously described in the literature (Ng et al., 2005;Zhou et al., 2010).
Data from the clinical investigation of Pfaff et al. (2015) was used to verify the adequate description of BL treatment of psoriasis by the model.

Multiple Parameter Sensitivity Analysis
A multiple parameter sensitivity analysis (MPSA) was performed to understand how the system may be affected by uncertainty in the model's parameters. An initial simulation for blue light treatment was executed with the nominal model parameters θ ref presented in Table 1. From this first simulation the reference output response y ref (k) was computed, where k represents the cell density of each keratinocyte cell densities derived by the model. 500,000 random parameter value sets were generated with the lhsdesign MATLAB function. The parameter sets θ (i) were used as input to the model and simulations were implemented. The cost function v(i) for the simulation results obtained with each parameter set was calculated as the sum of squared errors between the perturbed y sim and reference output y ref (Equation  19).
The parameter sets were then classified as acceptable or unacceptable and the cumulative frequency functions were computed for the acceptable and unacceptable values of each parameter. Then, the Kolmogorov-Smirnov statistic was calculated as the separation between S a (θ ) and S u (θ ) (Equation 20).
Finally, this value was used to derive the sensitivity of each parameter on the cell density of all keratinocyte populations at each of the six differentiation states and the local psoriasis severity index. Additionally, a local sensitivity analysis (LPSA; Marino et al., 2008) was performed for the parameters identified by the MPSA for further understanding of how small perturbations in these parameters reflect on the output of the model.

The Regulation of the Keratinocytes Proliferation Is Key in the Mechanism of Blue Light Irradiation
Given the phenomenological observations reported in the literature and the flexibility offered by the model to describe them, we first analyzed all the possible in silico representations of the changes induced on the proliferation and differentiation of the keratinocytes. Table 4 presents a summary of all potential cases. The general set of ODEs used in this analysis is presented by (Equations 1-17). Note that for cases where only θ BLγ is considered, θ BLk is equal to 1. Similarly when only θ BLk is considered, θ BLγ equals 1. For cases where only the stem cells and transit amplifying cells are affected, θ BLk for the differentiation rates of growth arrested, spinous, and granular cells equals 1. For each case, the time evolution of proliferation and differentiation parameters, the net proliferative and differentiation capacity of the keratinocytes, and the regulatory capacity of the immune system is computed (Figure 2). From this analysis, it is clear that shifting only parameters related to the symmetric and asymmetric division of proliferative keratinocytes does not yield a decrease on the proliferative compartment (Figures 2B,D,E). Further, changing the division parameters causes a marked increase in the differentiation rates (Figures 2B-D,F, H-M), which is reflected in an increased cell density of the transit amplifying cells in the proliferative compartment and all cells in the non-proliferative compartment. Modifying the selfproliferation rate of the proliferative keratinocytes leads to a decrease on the proliferation and increase on the differentiation of the cells in the proliferative compartment, without unrealistic surge of the keratinocytes differentiation capacity (Figure 2A). Thus, we conclude that the effect of blue light is best described in silico when θ BLγ is lower than 1 and θ BLk is equal to 1. The simulations and results presented in the following sections of this work are based on this conclusion.

A Transient Decrease in the Proliferative Capacity of Keratinocytes Lowers the Severity of Psoriatic Lesions
The results from Section The Regulation of the Keratinocytes Proliferation is Key in the Mechanism of Blue Light Irradiation suggested that the experimental observations for keratinocytes (Liebmann et al., 2010) are best described by the model when the self-proliferation rate of all stem cells and transit amplifying cells decreases during blue light irradiation. Simulations performed for the same conditions of that in vitro study show a similar reduction in the relative cell number achieved after BL irradiation at fluence of 33, 66, and 100 Jcm−2 ( Figure 3A). To verify that the model was also able to describe the reduced severity achieved by blue light treatment of psoriatic skin, simulations ( Figure 3B) were performed and compared to data from the clinical investigation of Pfaff et al. (2015). Their clinical study measured the LPSI of psoriatic patients after irradiating the skin with 90 Jcm−2 of pulsed blue light for 12 weeks at either a low (LI), 100 mWcm−2, or high (HI), 200 mWcm−2, peak intensities. During the first 4 weeks the patients were treated every day, the next 8 weeks the therapy occurred three times per week. In the simulations, the initial LPSI values are equal to those of the clinical investigation (Pfaff et al., 2015; i.e., 5.52 for LI and 5.17 for HI). The model achieved an accurate representation of the applied treatments at low and high peak

Case
Rate modified Value of blue light factors θ BLγ and θ BLk Population affected 1 Self-proliferation rates γ 1 and γ 2 are decreased. θ BLγ < 1 and θ BLk = 1 for proliferative cells Proliferative cells, i.e. stem cells and 2 Asymmetric (k 1a and k 2a ) and symmetric (k 1s and k 2s ) division rates of proliferative cells are increased.
θ BLγ < 1 and θ BLk > 1 (for the asymmetric division rates of proliferative cells) 7 The self-proliferation rates (γ 1 and γ 2 ) are decreased and the symmetric division rates (k 1s and k 2s ) are increased.
θ BLγ = 1 and θ BLk > 1 (for division rates of all cells) Proliferative cells, i.e. stem cells and transit amplifying cells. Non-proliferative 10 The self-proliferation rates (γ 1 and γ 2 ) are decreased and the differentiation rates (k 3 , k 4 , and k 5 ) of non-proliferative cells are increased.
θ BLγ = 1 and θ BLk > 1 (for the asymmetric division rates of proliferative cells and the division rates of non-proliferative cells)
θ BLγ = 1 and θ BLk > 1 (for the symmetric division rates of proliferative cells and the division rates of non-proliferative cells) 13 The self-proliferation rates (γ 1 and γ 2 ) are decreased and the division and differentiation rates (k 1a , k 1s, k 2a , k 2s , k 3 , k 4 , and k 5 ) are increased.
θ BLγ > 1 and θ BLk > 1 (for division rates of all cells) intensities, particularly during the first 4 weeks. In the next 8 weeks, the in silico results are less pronounced compared to the clinical study. The final values predicted by the model are within the error margins of the experimental results. Based on the results of the simulations performed for the experimental and clinical studies, it is evident that the transient decrease in proliferation has a pivotal role in the underlying mechanism of blue light treatment.
To assess the impact of changes in the proliferation rate and other model parameters on the outcome of the treatment, a multiple parameter sensitivity analysis was implemented. Figure 4 presents the results of this analysis for all model parameters in relation to the keratinocyte cell densities ( Figure 4A) and the LPSI (Figure 4B). According to Figure 4A, the model output is mainly affected by 12 parameters related to the proliferation of stem cells and transit amplifying cells. However, not all keratinocyte populations are evenly altered by variations in these parameters, some affect only the healthy populations while others have an impact on the diseased populations. The cell density of psoriatic keratinocytes is mainly affected by the proliferation rate of healthy and diseased transit amplifying cells. Conversely, the proliferation rate of psoriatic stem cells only impacts the healthy populations and the psoriatic stem cells. The LPSA performed for these 12 sensitive parameters indicated that diseased transit amplifying cells are most sensitive to small perturbations on the proliferation rate of healthy transit amplifying cells. This parameter as well as the proliferation rate of diseased stem cells and transit amplifying cells are directly related to the fluence, intensity and irradiation time used by the implemented treatment. From the 12 parameters identified in the MPSA for keratinocytes cell density, only seven have a strong impact on the local severity of psoriasis at the end of treatment (Figure 4B). These seven parameters correspond to those affecting the cell densities of diseased keratinocyte populations. The impact of these parameters on the LPSI is consistent with the effect on the cell density of diseased keratinocytes, but lower in comparison to the latter.
Based on the sensitivity analysis, it was confirmed that the decreased proliferation capacity induced by blue light irradiation directly impacts all keratinocytes populations. Despite the insightful information this analysis provided, it was yet unclear how each keratinocyte population varied before, during, and after shining blue light on the skin. Thus, a simulation ( Figure 5) was executed for the blue light treatment of a lesional skin area of 1 mm 2 for 84 days, using a fluence of 90 Jcm −2 . The simulated scheme was divided into two sections, the first one with everyday treatment, and the second one with three times per week irradiations. The time evolution of both healthy and   Figures 5A,B respectively. The model predicted that the cell density of all keratinocyte populations decreases due to the repetitive irradiation with blue light. However, the cell density of diseased keratinocytes remained considerably high after the treatment compared to the cell density of healthy keratinocytes. Thus, no shift to a healthy phenotype was achieved by this therapeutic approach. The model suggested that a lesional state may prevail after the treatment, but some improvements might be observed within the period of phototherapy. Additionally, it predicted that a less pronounced decrease in the cell density of all keratinocytes is achieved when the treatment occurs every other day compared to daily treatment. This observation suggests that the longer the period of daily treatment the lower the final LPSI value.

Blue Light Treatment Schemes with Long Duration and High Fluence Yield High Efficacy in the Management of Psoriasis
Hitherto, the model predicted that regulating the hyperproliferative populations of keratinocytes induces a transient management of psoriasis. However, it also showed that the treatment scheme used in the simulations has a direct impact in the predicted outcome. This observation agrees with the large variations perceived in the results from the available clinical investigations (Maari et al., 2003;Weinstabl et al., 2011;Kleinpenning et al., 2012;Pfaff et al., 2015). Therefore, we studied the impact of each treatment parameter on the efficacy of the therapeutic approach. Figure 5B shows a prominent decrease in the cell densities of diseased keratinocytes when the blue light is constantly shinned upon the skin compared to when the irradiation is less frequent. To further analyze this observation, simulations were executed for therapy schemes with a total duration of 4-28 weeks, using either a daily or every other day treatment (Figure 6A). The results showed that the longer the period of daily treatment the lower the final LPSI value. Nevertheless, the relative change in the LPSI value decreased for treatment periods longer than 20 weeks. Further, BLISS predicted that for a blue light phototherapy scheme where both daily and every other day treatment sessions are used, higher efficacy is achieved with an increasing number of daily treatment sessions.
Given that repeated treatment sessions seemed to yield higher efficacy, we then tested whether the same trend applied to the irradiation mode of treatment scheme. As previously mentioned a given fluence of blue light can be shinned on the skin in either a continuous or pulsed mode. Both cases comprise the same average intensity, however the peak intensity differs. Pfaff et al. studied in vivo the impact of low and high peak intensities in the pulsed irradiation of psoriatic skin (Pfaff et al., 2015). They observed minor differences between low and high peak intensity groups; however this trend was not seen in silico (Figure 3B solid lines vs. dots). In the model, both low and high intensity pulsed treatment led to the same trend. These results led to questioning whether the model could predict a different outcome for treatment strategies with pulsed and continuous irradiation. Thus, simulations were implemented for low and high power densities in pulsed and continuous modes (Figures 6B-D), considering a 12 weeks treatment and 4 weeks follow-up with an average fluence of 90 Jcm −2 and an average intensity of 50 mWcm −2 for both irradiation modes. Figure 6B compares the efficacy of pulsed and continuous mode in terms of the local psoriasis severity index. In the simulations for the pulsed mode, low (100 mWcm −2 ) and high (200mWcm −2 ) peak intensities were used. No differences were observed between both modes, or between low and high peak power densities. Hence, the model suggested that duration rather than intensity determines the efficacy of the treatment. This prediction was further tested by simulating the treatment of a psoriatic epidermis using power  densities of 100, 200, and 600 mWcm −2 in the continuous mode (Figures 6C,D).
In addition to time and intensity, fluence is another important parameter of the treatment scheme. Simulations performed for fluences between 0 and 750 Jcm −2 , considering a treatment of 12 weeks ( Figure 7A) predicted that fluence levels as low as 45 Jcm −2 lead to a decrease in the LPSI value. Moreover, at fluence levels between 200 and 500 Jcm −2 a saturation point was observed. At fluences equal or higher than 500 Jcm −2 the decrease in the LPSI value was considerably larger due to the cytotoxicity induced by BL at these levels. According to these results, low fluences would lead to a minor decrease in the LPSI compared to higher fluences. This observation is better understood when analyzing the relation between the fluence, the blue light factor θ BLγ , and the consequent proliferation rate derived by the model (Figure 7B). Here, we show only the analysis for the proliferation rate of diseased transit amplifying cells, but the trends are consistent for the other proliferative populations. A clear exponential decrease is seen on both the blue light factor and the proliferation rate with increasing fluences. However, the impact of fluence on the proliferation rate is steeper and faster compared to the effect on the blue light factor. Further, for fluences below 500 Jcm −2 both the blue light factor and the proliferation rate are reduce in an almost linear manner. But it is less abrupt at fluences above 200 Jcm −2 . The fluence effect translates into a linear reduction of the proliferation in relation to the blue light factor, which yields a lower cell density and the consequent improvement on the psoriatic skin. The model was designed to allow bi-stability at either a diseased or healthy phenotype. The healthy state can only be reached when the total cell density at the end of treatment is dominated by healthy keratinocyte populations. To achieve the switch in the ruling type of population, the number of diseased cells must be considerably lowered during the treatment. Our results showed that this shift in phenotypes could not be acquired with blue light fluences where no cytotoxicity was induced (Figures 7C,D). Figure 7C presents the relation between healthy and diseased stem cells with increasing fluence levels. From this figure is clear that only fluences above 500 Jcm-2 achieve a majority of healthy stem cells. Fluences below that level reduce the cell densities but remain within the diseased state due to a high predominance of diseased keratinocytes ( Figure 7D).

DISCUSSION
Blue light at 420 and 453 nm is highly effective in the management of psoriasis (Weinstabl et al., 2011;Pfaff et al., 2015) without any adverse effects (Awakowicz et al., 2009;Kleinpenning et al., 2010). We have explored the underlying mechanism of this therapeutic method and the impact of the treatment scheme on the efficacy of the treatment using an in silico approach. BLISS suggests that the decrease in the keratinocytes proliferative capacity is sufficient to induce a transient decrease in the severity of psoriasis. Further, FIGURE 6 | The length of treatment is a key factor to achieve high treatment efficacy. (A) The LPSI at the end of treatment is predicted for treatment protocols with a total duration of 4-28 weeks. The treatment sessions in the simulated protocols occur either on a daily (blue dots) or every other day (green dots) basis. (B) A 12 weeks treatment with fluence of 90 Jcm −2 is simulated for continuous (black) and pulsed irradiation with low (blue) and high (red) peak intensities. (C) The variation in the LPSI value is predicted for irradiation schemes with low and high intensities of continuous wave blue light. (D) The evolution of the total keratinocyte cell densities during (blue bar) and after (black bar) treatment are derived using the same low and high intensities of continuous wave blue light presented in panel (C).
it predicts that treatment schemes with long duration and high fluence yield the most optimal outcome. The model constitutes a fast and flexible tool that considers not only the properties of the irradiated epidermis but also the interactions between the structural cells within it. Our in silico method provides a framework for further insight on the underlying mechanism of this blue light-based phototherapy and the optimal treatment of psoriasis. To the extent of our knowledge, this is the first model that studies the changes induced by blue light irradiation on the cell dynamics of keratinocytes. The general concepts of epidermal kinetics used here have been applied to UV phototherapy (Grabe and Neuber, 2007;Weatherhead et al., 2011;Zhang et al., 2015) before. Nevertheless, the differences on the underlying mechanism make the previous models and ours two utterly different applications.
On this work, we first addressed how to best represent these phenomenological observations in silico while accounting for the characteristic phenotype of psoriasis (Perera et al., 2012). From experimental studies, it is known that blue light decreases the proliferative capacity of keratinocytes and other cell types (Taoufik et al., 2008;Wataha et al., 2008), and that it reduces their proliferation by inducing their differentiation (Liebmann et al., 2010). However, it is not clear how exactly these changes are exerted on the cellular processes of keratinocytes. BLISS accounts for the three main theories of cell division (Clayton et al., 2007;Klein et al., 2011;Watt, 2014), i.e., self-proliferation, asymmetric, and symmetric division. Selfproliferation refers to the ability of a cell dividing into two daughter cells of the same kind. Asymmetric division considers that one proliferative cell divides intro one daughter cell of the same kind and one of the next differentiation stage. Finally, symmetric division comprises the division of a proliferative cell into two cells of the next differentiation stage. In the model the proliferation and differentiation of keratinocytes in the proliferative compartment is represented by their selfproliferation, symmetric, and asymmetric division parameters. The keen reader could think that in order to best describe the effect of blue light on proliferation and differentiation, increasing the rates of asymmetric and symmetric division should be enough to achieve similar observations as those reported in vitro. Alternatively, the direct increase on the division rates of proliferative cells should be coupled to the decrease of the self-proliferation rate of all keratinocytes in the proliferative compartment. Thus, capturing the increase in differentiation and the decrease in proliferation of the irradiated keratinocytes. Our analysis (Figure 2) showed that modifying the division rates alone does not appropriately reflect on the proliferative capacity of the cells. Additionally, it translates into high peaks in the cell density of transit amplifying and non-proliferative cells, which are biologically unlikely. However, decreasing only the proliferation rate does yield the effect on the differentiation and proliferative capacity of keratinocytes observed experimentally. The outcome of this approach is a transient decline in the proliferative rates of the cells and consequent higher differentiation, translated into a temporary drop local severity of the disease. BLISS comprises a fine tuning between the proliferation and division rates of all keratinocytes, particularly of healthy stem cells. It considers a strong connection between the proliferation and division rates of healthy stem cells and their regulation by the total population of healthy and diseased transit amplifying cells. This fact might explain the outcome of our analysis. From the sensitivity analysis (Figure 4), it was clear that proliferation-related parameters had the highest influence on the final cell densities and the improvement of psoriasis due to blue light treatment. The most sensitive parameters were directly related to the treatment scheme. Their impact was reflected on the cell density of diseased keratinocyte populations and local psoriasis severity index. These observations may explain the large variations observed among clinical investigations of BL treatment (Maari et al., 2003;Weinstabl et al., 2011;Pfaff et al., 2015).
Previous models of phototherapy, particularly UV-based approaches represented the therapeutic effect by inducing strong apoptosis on the keratinocytes (Weatherhead et al., 2013;Zhang et al., 2015). This temporal increase in the apoptotic rate leads to a remission phase, characterized by a lower total cell density after the treatment compared to the initial one. In our model, apoptosis is only affected at fluences above 500 Jcm −2 , instead the proliferation and differentiation capacity of keratinocytes are altered during the treatment sessions. The key difference in the approaches of representing UV and BL phototherapy in silico is defined by the difference in the underlying mechanisms of both therapeutic approaches. Despite both methods affecting the cells behavior by absorption of light at a certain wavelength, the photo-acceptors and consequent events differ. UV is mainly absorbed by the cell's DNA (Benedix et al., 2005), consequently inducing direct DNA damage, apoptosis (Weatherhead et al., 2011), and a decrease in the total cell density of the lesional skin. Blue light irradiation does not induce DNA damage (Awakowicz et al., 2009), it is absorbed by other photo-acceptors in the cell, e.g., flavins (Sadeghian et al., 2008) and cryptochromes (Bouly et al., 2007). The consequent effect is the release of nitric oxide (Opländer et al., 2013) and reactive oxygen species (Yoshida et al., 2013), which alter the cell's proliferation and differentiation rates.
Clinically, this transient effect on the keratinocytes proliferation and differentiation is observed as the improvement of a psoriatic lesion (Weinstabl et al., 2011;Kleinpenning et al., 2012). The model is bi-stable, thus it can reproduce steady states for both healthy and diseased skin phenotypes. In our simulations, we account for mild psoriasis through a specific set of initial cell densities and show that although the plaque is controlled during the treatment the underlying state is still diseased. After the treatment, the model eventually goes back to the initial diseased state. In the UV phototherapy model of Zhang et al. (2015), they demonstrate that UV treatment with consecutive episodes can theoretically control the disease and achieve a remission of the psoriatic phenotype. In our model, we do not observed this remission phase for mild psoriasis considered with fluences below 500 Jcm −2 .The lack of remission phase in our simulations for blue light treatment in contrast with UV phototherapy may be due to the underlying mechanism. Increasing the apoptosis of keratinocytes yields a lower cell density regardless of their proliferation and differentiation rate. This effect results in a shift of the stem cell density from the disease state to the healthy state ( Figure 7C). Altering the proliferation and differentiation rates of the proliferative keratinocytes due to blue light irradiation only impacts their apoptosis rates at high fluence levels. Thus, the effect of blue light is mild compared to that of UV (Weatherhead et al., 2013). Nevertheless, it is possible that other factors are involved in the underlying mechanism of blue light management of psoriasis. The inclusion of these factors in the model could potentially yield to the shift in phenotypes and full remission of psoriatic skin. In our method, the immune system regulates the dynamics of the keratinocytes through the removal of disease stem cells at a rate defined by the maximum immune killing rate and the half activation of the immune system (Equation 7). Blue light also induces apoptosis of T-cells (Liebmann et al., 2010;Oh et al., 2016), alter the immune response (Fischer et al., 2013) and decrease inflammation (Shnitkind et al., 2006). However, the impact of BL on the immune system is not considered in the current model. The inclusion of blue light related effects on the immune system might result in lower keratinocyte cell densities and the consequent shift in the phenotype. The available clinical investigations have not assessed yet the impact of blue light on the immune system in the context of psoriasis. Data from skin biopsies, experimental, and clinical studies on blue light irradiation is needed to determine which aspects could be involved and how they should be included in the model.
BLISS explicitly includes relevant treatment parameters, i.e., fluence, time, intensity, and mode of irradiation. This feature of the model allows us to investigate the influence of these parameters on the management of psoriasis. Recent experimental studies have explored the effect of various wavelengths on the blue range of the electromagnetic spectrum (Liebmann et al., 2010;Opländer et al., 2013), and fluence on various cell types (Awakowicz et al., 2009;Liebmann et al., 2010;Sparsa et al., 2010;Monfrecola et al., 2014). Clinically, the limited number of available studies focuses on the assessment of treatment efficacy with a specific combination of parameters. Only two studies implements a wavelength comparison (Weinstabl et al., 2011;Kleinpenning et al., 2012), and another evaluates high and low pulsed intensities (Pfaff et al., 2015). Nevertheless, no global guidelines for treatment have been defined. Furthermore, the impact on the kinetics of the main cells involved in psoriasis and their role on the management of the disease are yet to be resolved. One of the objectives of this study was to study the impact of treatment parameters on the blue light mediated regulation of keratinocytes in psoriasis and the consequent control of the lesion. Our results suggested that therapy length and the applied fluence rather than intensity of irradiation may potentially determine the efficacy of the treatment. Simulations for various treatment lengths ( Figure 6A) showed that algorithms with a long treatment period and consecutive episodes may be the most effective for psoriasis. This observation could partially explain the higher improvement of Pv symptoms achieved by Pfaff et al. (2015) (12 weeks) in contrast with prior investigations (4 weeks; Maari et al., 2003;Weinstabl et al., 2011). Another reason could be the different fluences used in these studies [10 Jcm −2 (Maari et al., 2003) and 90 Jcm −2 (Weinstabl et al., 2011;Pfaff et al., 2015)], given the fluence dependency associated with blue light irradiation (Liebmann et al., 2010;Dai et al., 2012). Simulations implemented for a wide range of fluences (Figure 7) show a negative correlation between fluence and final LPSI value. However, no threshold value has yet been defined. To define this value additional in vitro and in vivo data are needed. Research on light therapy for different applications has indicated that pulsed light is more effective than continuous light, while others report no effect, or a worsening effect of pulsed irradiation compared to no treatment (Lapchak et al., 2007;Hashmi et al., 2010). In Psoriasis, only one clinical investigation has used a pulsed irradiation mode (Pfaff et al., 2015), while the others employed the continuous mode (Maari et al., 2003;Weinstabl et al., 2011;Kleinpenning et al., 2012). Simulations performed with the model (Figure 6B) indicated that pulsed light is as effective as continuous light in the treatment of psoriasis.
The current model has some addressable limitations: (i) it does not explicitly comprise the photo-activated processes leading to the cellular regulation and management of psoriasis.
(ii) We neglected the additional cytotoxicity induced on the keratinocytes at wavelengths shorter than 453 nm (Liebmann et al., 2010). The model was built assuming irradiation at 453 nm, which is the most common wavelength used for psoriasis (Weinstabl et al., 2011;Pfaff et al., 2015). Despite the minimum differences in treatment efficacy reported in the literature (Kleinpenning et al., 2010;Weinstabl et al., 2011), inclusion of this effect could be worth a thorough analysis. However, additional data is needed on the optical parameters and the skin before this feature can be included in the model.
One definite experimental test of our model should aim at simultaneously tracing psoriasis and healthy proliferative and non-proliferative keratinocytes to verify whether blue light equally affects healthy and diseased proliferative cells. One potential approach to test this hypothesis could be isolating keratinocytes from psoriatic skin (Aasen and Belmonte, 2010) after shinning blue light on it. In this work, we provide general insights and recommendations based on the model results. In the future, more precise recommendations could be derived from the model by implementing an optimization analysis for the treatment (Kessel et al., 2013). Our results showed that length of treatment is vital in decreasing the severity of a psoriasis lesion, however there is little information reported in the literature on the direct effect of irradiation length and cumulative exposure of cells to blue light. Future studies should explore this theory experimentally to define the boundaries of the cellular response and the most optimal outcome. Also, they should assess the long term management of psoriasis and the recurrence of the lesions.

CONCLUSIONS
Overall, this study demonstrates that regulation of the proliferative capacity of keratinocytes is a crucial process in the blue light induced management of psoriasis. Considering the uprising interest in blue light as treatment for inflammatory skin conditions, the availability of irradiation guidelines becomes of great relevance. Our model constitutes a flexible tool for studying the underlying mechanism, formulating valuable recommendations, and designing effective therapeutic regimes with BL in a constraint-free environment.

FUNDING
This project was funded under the public-private partnership between Philips and the Eindhoven University of Technology.