Front. Physiol. Frontiers in Physiology Front. Physiol. 1664-042X Frontiers Media S.A. 10.3389/fphys.2017.00028 Physiology Original Research A Dynamic Model for Prediction of Psoriasis Management by Blue Light Irradiation Félix Garza Zandra C. 1 * Liebmann Joerg 2 Born Matthias 2 Hilbers Peter A. J. 1 van Riel Natal A. W. 1 1Department of Biomedical Engineering, Eindhoven University of Technology Eindhoven, Netherlands 2Philips Technologie GmbH, Innovative Technologies Aachen, Germany

Edited by: Matteo Barberis, University of Amsterdam, Netherlands

Reviewed by: Marcel Schilling, German Cancer Research Center, Germany; Tianhui Niu, The General Hospital of Air Force, China

*Correspondence: Zandra C. Félix Garza z.c.felix.garza@tue.nl

This article was submitted to Systems Biology, a section of the journal Frontiers in Physiology

 26 01 2017 2017 8 28 27 10 2016 11 01 2017 Copyright © 2017 Félix Garza, Liebmann, Born, Hilbers and van Riel. 2017 Félix Garza, Liebmann, Born, Hilbers and van Riel

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

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.

 inflammatory skin conditions epidermis keratinocytes computational model phototherapy visible light Philips10.13039/100004320 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 light-based 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.

Schematic description of the mechanistic model for blue light treatment of psoriasis. (A) The model considers the 4 sub-layers of the epidermis, the 6 stages of differentiation for keratinocytes across the sub-layers. (B) BLISS accounts for the cellular processes of proliferation (γ), differentiation (k), apoptosis (β), and desquamation (α).

 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 (Psch, d), transit amplifying cells (Pah, d), growth arrested cells (Pgah, d), spinous cells (Psph, d), granular cells (Pgch, d), corneocytes (Pcch, 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 (kh, 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 (kah, d)], or induce two daughter cells where both correspond to the next stage of differentiation [symmetric division (ksh, 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).

dPschdt=[γ1hθBLγ(1-Psch+PscdλPscmax)-k1shθBLk-β1h]Psch+k-1h Ptah dPtahdt=(γ2hθBLγ-k2shθBLk-β2h-k-1h)Ptah-(k1ah+2k1sh) θBLkPsch+k-2hPgah dPgahdt=(k2ah+k2sh)θBLkPtah-(k-2h+k3hθBLk+β3h)Pgah dPsphdt=k3hθBLkPgah-(k4hθBLk+β4h)Psph dPgchdt=k4hθBLkPsph-(k5hθBLk+β5h)Pgch dPcchdt=k5hθBLkPgch-ahPcch dPscddt=[γ1dθBLγ(1-Psch+PscdλPscmax)-k1sdθBLk-β1d]Pscd-KpPscd2Ka2+Pscd2+k-1dPtad dPtaddt=(γ2dθBLγ-k2sdθBLk-β2d-k-1d)Ptad-(k1ad+2k1sd)θBLkPscd+k-2dPgad dPgaddt=(k2ad+k2sd)θBLkPtad-(k-2d+k3dθBLk+β3d)Pgad dPspddt=k3dθBLkPgad-(k4hθBLk+β4d)Pspd dPgcddt=0 dPccddt=k5dθBLkPgcd-adPccd

BLISS model parameterization.

Parameter Description Value Source
Pscmax Growth capacity of stem cells 4.5 × 103 mm−2 Zhang et al., 2015
γ1hom Minimal stem cell self-proliferation rate constant 3.30 × 10−3 d−1 Zhang et al., 2015
k1shom Minimal symmetric healthy stem cell division rate constant 1.64 × 10−3 d−1 Clayton et al., 2007; Zhang et al., 2015
k1ahom 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
k2sh 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
k2ah 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
k3h 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
k4h 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
k5h 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−1h Back conversion rate constant of healthy cells (transit amplifying to stem cells) 1.00 × 10−6 d−1 Zhang et al., 2015
k−2h Back conversion rate constant of healthy cells (growth arrested to transit amplifying cells) 1.00 × 10−6 d−1 Zhang et al., 2015
ω Maximum fold increase of stem cells proliferation rate 100 Heenen et al., 1998; Zhang et al., 2015
n Stem cells proliferation rate regulation by transit amplifying cells 3 Zhang et al., 2015
AIh Epidermal apoptosis index for healthy skin 0.12% Bauer et al., 2001; Hoath and Leahy, 2003; Zhang et al., 2015
β1h Apoptosis rate of healthy epidermal stem cells 1.97 × 10−6 d−1 Calculated as described by Equation 16
β2h Apoptosis rate of healthy transit amplifying cells 2.08 × 10−5 d−1 Calculated as described by Equation 16
β3h Apoptosis rate of healthy growth arrested cells 2.60 × 10−4 d−1 Calculated as described by Equation 16
β4h Apoptosis rate of healthy spinous cells 6.68 × 10−5 d−1 Calculated as described by equation 16
β5h Apoptosis rate of healthy granular cells 1.33 × 10−4 d−1 Calculated as described by Equation 16
αh Healthy corneocytes desquamation rate constant 7.14 × 10−2 d−1 Weinstein et al., 1985; Bauer et al., 2001; Hoath and Leahy, 2003; Zhang et al., 2015
ρsc Fold change of stem cells proliferation in psoriasis 4 Weinstein et al., 1985; Zhang et al., 2015
ρta Fold change of transit amplifying cells proliferation in psoriasis 4 Weinstein et al., 1985; Zhang et al., 2015
ρtr Fold change of psoriatic cells transit rate 5 Weatherhead et al., 2011; Zhang et al., 2015
ρde Fold change of psoriatic corneocytes desquamation 4 Weinstein and Van Scott, 1965; Zhang et al., 2015
λ Fold change of stem cells growth capacity in psoriasis 3.5 Heenen et al., 1987; Simonart et al., 2010; Zhang et al., 2015
Kp Maximum immune response rate 6 Zhang et al., 2015
Ka Half-activation of immune system by psoriatic stem cells density 380 Zhang et al., 2015
γ1d Diseased stem cell self-proliferation rate constant 1.16 × 10−2 d−1 Calculated as the product of γ1hom and ρsc
k1sd Symmetric diseased stem cell division rate constant 5.70 × 10−3 d−1 Calculated as the product of k1shom and ρsc
k1ad Asymmetric diseased stem cell division rate constant 4.59 × 10−2 d−1 Calculated as the product of k1shom and ρsc
γ2d Diseased transit amplifying cells self-proliferation rate constant 4.90 × 10−2 Calculated as the product of γ2h and ρta
k2sd Diseased transit amplifying cells symmetric division rate constant 6.06 × 10−2 d−1 Calculated as the product of k2sh and ρta
k2ad Diseased transit amplifying cells asymmetric division rate constant 4.83 × 10−1 d−1 Calculated as the product of k2ah and ρta
k3d Diseased growth arrested to spinous cells differentiation rate constant 9.72 × 10−1 d−1 Calculated as the product of k3h and ρtr
k4d Diseased spinous to granular cells differentiation rate constant 2.50 × 10−1 d−1 Calculated as the product of k4h and ρtr
k5d Diseased granular cells to corneocytes differentiation rate constant 3.89 × 10−1 d−1 Calculated as the product of k5h and ρtr
AId Epidermal apoptosis index for psoriatic skin 0.035% Bauer et al., 2001; Hoath and Leahy, 2003; Zhang et al., 2015
β1d Apoptosis rate of diseased epidermal stem cells 2.01 × 10−6 d−1 Calculated as described by Equation 17
β2d Apoptosis rate of diseased transit amplifying cells 2.12 × 10−5 d−1 Calculated as described by Equation 17
β3d Apoptosis rate of diseased growth arrested cells 3.40 × 10−4 d−1 Calculated as described by Equation 17
β4d Apoptosis rate of diseased spinous cells 8.76 × 10−5 d−1 Calculated as described by Equation 17
β5d Apoptosis rate of diseased granular cells 1.36 × 10−4 d−1 Calculated as described by Equation 17
αd Diseased corneocytes desquamation rate constant 2.50 × 10−1 d−1 Calculated as the product of αh and ρde
aγ Blue light coefficient for proliferation factor 1 Estimated from Liebmann et al., 2010
bγ Blue light coefficient for proliferation factor −3.40 × 10−3 Estimated from Liebmann et al., 2010
ak Blue light coefficient for differentiation factor 2.46 Estimated from Liebmann et al., 2010
bk Blue light coefficient for differentiation factor 1.94 × 10−2 Estimated from Liebmann et al., 2010
ck Blue light coefficient for differentiation factor 3.46 Estimated from Liebmann et al., 2010
θBL β500−750 Blue light factor increasing the apoptosis rate at fluences higher than 500 Jcm−2 and lower than 750 Jcm−2 3.9 × 10−2 Awakowicz et al., 2009
θBL β>750 Blue light factor increasing the apoptosis rate at fluences higher than 750 Jcm−2 5 × 10−2 Awakowicz et al., 2009
ξabs Energy absorbance of the epidermis for a low perfused Caucasian skin 57.9% Calculated from optical model

Where

γ1hγ ̀1,hom=k1ahk1a,hom=k1shk1s,hom=ω1+(ω-1)(Ptah + PtadPta,hom)n

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 hyper-proliferative 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 Pscmax available in the epidermis, which limits the numbers of healthy and diseased stem cells Psch,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 Kp and the half activation of the immune system Ka due to psoriatic stem cells. The immune response is ample when the cell density of diseased cells exceeds the threshold defined by Ka.

In Equation (13), γ1,hom, k1a,hom, k1s,hom are the homeostatic rate constants for each division process, and Pta,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, k2sd, and k2ad) are ρta fold changes higher than the rates of healthy transit amplifying cells. Similarly the diseased cells in the non-proliferative compartment, i.e., Pgad and Pspd, 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).

θBLγ = aγe(bγF)                    with aγ>0 and bγ<0

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 proliferation and increase of differentiation rates were initially estimated by fitting in vitro data from Liebmann et al. (2010).

θBLk = ck-ake(bkF)                      with ak,bk and ck >0

Blue light phototherapy is divided into sessions with a certain fluence F, defined by the product of the irradiation time and the average intensity Iav 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 (k1sh, k2sh, k3h, k4h, k5h) and the probability of a cell undergoing apoptosis (apoptotic index AIh). 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.

βh= AIhkj1-AIh+θBLβ               j∈{1sh,2sh,3h,4h,5h}

The diseased keratinocytes have an apoptosis rate defined by the apoptotic index of psoriatic cells AId, the differentiation rate of each population (k1sd, k2sd, k3d, k4d, k5d) and qBLβ. (Equation 17).

βd= AIdkj1-AId+θBLβ               j∈{1sd,2sd,3d,4d,5d}

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 five-layered 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 mm2, which is 2 times higher than that of healthy skin (~100,000 cells per mm2; Hoath and Leahy, 2003). BLISS is available at GitHub1 and the BioModels Database2 (Chelliah et al., 2015) with identifier MODEL1701090001.

Blue light irradiation parameters.

Parameter Value
Fluence 90 Jcm−2
Irradiation mode Continuous or pulsed
Irradiation time 1800 s (30 min)
Duration of treatment 84 days

Initial cell density distribution (cells per mm2).

Keratinocytes population Healthy Psoriatic
Corneocytes 185 77,633
Granular cells 119 0
Spinous cells 238 79,788
Growth arrested cells 61 20,536
Transit amplifying cells 77 32,098
Stem cells 362 6459
Total 1042 + 216514 = 217,556

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 Ip, which is equal to the average power density Iav. In PW, the epidermis is irradiated in short pulses at an Ip higher than Iav. 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 Ip 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).

LPSIi=LPSI0-LPSI0(PToti - PTot0PTot0)             where i≥1

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 yref(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 ysim and reference output yref (Equation 19).

v(i) = ∑k = 1N(yref(k)−ysim(k,θ(i))​)2

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 Sa(θ) and Su(θ) (Equation 20).

D=max|Sa(θ)−Sq(θ)|

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.

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

Possible cases describing the blue light effects on the cellular processes of keratinocytes.

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 (k1a and k2a) and symmetric (k1s and k2s) division rates of proliferative cells are increased. θBLγ = 1 and θBLk > 1 for proliferative cells transit amplifying cells.
3 The self-proliferation rates (γ1 and γ2) are decreased and the asymmetric and symmetric division rates (k1a, k2a, k1s, and k2s) are increased. θBLγ < 1 and θBLk > 1 for proliferative cells
4 Asymmetric (k1a and k2a) division rates are increased. θBLγ = 1 and θBLk > 1 (for the asymmetric division rates of proliferative cells)
5 Symmetric (k1s and k2s) division rates are increased. θBLγ = 1 and θBLk > 1 (for the symmetric division rates of proliferative cells)
6 The self-proliferation rates (γ1 and γ2) are decreased and the asymmetric division rates (k1a and k2a) 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 (k1s and k2s) are increased. θBLγ < 1 and θBLk > 1 (for the symmetric division rates of proliferative cells)
8 Differentiation rates (k3, k4, and k5) of non-proliferative cells are increased. θBLγ = 1 and θBLk > 1 (for division rates of non-proliferative cells) Non-proliferative cells, i.e. growth arrested cells, spinous cells, and granular cells.
9 Differentiation rates (k1a, k1s, k2a, k2s, k3, k4, and k5) of all cells 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 (k3, k4, and k5) of non-proliferative cells are increased. θBLγ < 1 and θBLk > 1 (for division rates of non-proliferative cells) cells, i.e. growth arrested cells, spinous cells, and granular cells.
11 The asymmetric division and differentiation rates (k1a, k2a, k3, k4, and k5) are increased. θBLγ = 1 and θBLk > 1 (for the asymmetric division rates of proliferative cells and the division rates of non-proliferative cells)
12 The symmetric division and differentiation rates (k1s, k2s, k3, k4, and k5) are increased. θ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 (k1a, k1s, k2a, k2s, k3, k4, and k5) are increased. θBLγ > 1 and θBLk > 1 (for division rates of all cells)

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

Representing in silico the blue light effect by decreasing the proliferation rate of keratinocytes yields the behavior observed experimentally. The panels show the changes on the proliferation and differentiation of healthy and diseased keratinocytes due to different approaches of blue light in silico representation. Thirteen different cases are depicted in this figure. A detailed description of each case can be found in Table 4. These cases are grouped according to the parameters affected by blue light. Panels (A–G) show cases where only the proliferation and division parameters from the proliferative compartment are modified during blue light irradiation. Panel (H) presents the case where only the differentiation parameters from the non-proliferative compartment are affected. Panels (I–M) described the cases where blue light irradiation impacts both compartments. The upper plots of each panel show the impact of the blue light factor on the proliferation and differentiation rates of all keratinocyte populations as a consequence of the parameters affected by blue light. The lower panels display the impact of each potential representation on the proliferative (green), differentiation (blue), and immune system-cytotoxic (red) capacities.

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

BLISS accurately describes in vitro and in vivo data. (A) The relative keratinocyte density is computed (solid black line) for irradiation schemes with fluences between 0 and 100 Jcm−2 and compared to the experimental data (Liebmann et al., 2010) of keratinocytes irradiated with three fluences of continuous wave blue light (gray dot). (B) The changes in the severity of a psoriasis lesion during (blue bar) and after treatment (black bar) are predicted (solid lines) for a pulsed wave treatment at high (red) and low (blue) intensities. These changes in the severity of the disease are shown in terms of the local psoriasis severity index (LPSI). The in silico results are compared to the clinical data (dots) of Pfaff et al. (2015).

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.

Multiple parameter sensitivity shows that variations in the proliferation parameters have a high impact on the final keratinocyte cell densities and the local severity of psoriasis at the end of treatment. (A) The sensitivity of healthy and diseased keratinocytes (x axis) to variations in the 58 model parameters (y axis) is depicted in the heat map. Only those with a Kolmogorov–Smirnov statistic value higher than 0.3 are shown in the y axis. (B) The sensitivity of the treatment outcome to changes in the same model parameters of Panel (A) are shown in the bar plot in terms of the local psoriasis severity index. Only those parameters with a Kolmogorov-Smirnov statistic value higher than 0.3 are shown in the x axis.

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 mm2 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 diseased keratinocytes is depicted in 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 induces a temporary reduction of all keratinocyte populations during the length of the treatment. (A) The time evolution of the healthy keratinocyte cell densities (dashed lines) comprised in a psoriatic epidermis are described during (blue bar) and after (black bar) blue light-based therapy. (B) Similarly, the time evolution of the diseased keratinocytes (continuous lines) is shown during and after treatment. Note the difference in scales of the y-axes between panels (A) and (B). In this simulation, the treatment periods divided into two sections. The first one consists of 28 days, with irradiation sessions occurring every day. The second section comprises 56 days, with treatment sessions happening three times per week. On each session a fluence of 90 Jcm−2 is applied for 30 min.

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

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

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

High fluences of blue light yield higher management of psoriasis than low fluences. (A) The LPSI over time is computed for a 12 weeks continuous wave treatment with an initial LPSI of 5.17 (Pfaff et al., 2015) using 11 fluences. (B) The behavior of the blue light factor θBLγ and the proliferation rate are visualized as functions of the fluence. (C) The bi-stability of the model is evaluated in the logarithmic scale for a wide range of fluences in a phase diagram of healthy and diseased keratinocytes. Only fluences above 500 Jcm−2 yield a shift the dominant keratinocyte population and the phenotype of the skin. (D) Fluences below this level remain in the disease state.

 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, 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. Self-proliferation 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 self-proliferation, 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.

 Author contributions

Wrote the paper: ZCFG. Conceived and designed the computational framework: ZCFG, JL, MB, PAJH, NAWvR. Developed the model: ZCFG. Performed the computational analysis: ZCFG. Revised the paper: JL, MB, PAJH, NAWvR. Supervised the study: JL, MB, PAJH, NAWvR.

 Funding

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

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

We are grateful to Elvira Paulussen for her contribution to the inclusion of optical parameters in the model. The model described here was deposited in the online repository BioModels Database (Chelliah et al., 2015) with identifier MODEL1701090001.

 References Aasen T. Belmonte J. C. I. (2010). Isolation and cultivation of human keratinocytes from skin or plucked hair for the generation of induced pluripotent stem cells. Nat. Protoc. 5, 371–382. 10.1038/nprot.2009.24120134422 Awakowicz P. Bibinov N. Born M. Busse B. Gesche R. Helmke A. . (2009). Biological stimulation of the human skin applying healthpromoting light and plasma sources. Contrib. Plasma Phys. 49, 641–647. 10.1002/ctpp.200910068 Bauer J. Bahmer F. A. Wörl J. Neuhuber W. Schuler G. Fartasch M. (2001). A strikingly constant ratio exists between Langerhans cells and other epidermal cells in human skin. A stereologic study using the optical disector method and the confocal laser scanning microscope. J. Invest. Dermatol. 116, 313–318. 10.1046/j.1523-1747.2001.01247.x11180009 Benedix F. Berneburg M. Röcken M. (2005). Phototherapy with narrowband vs broadband UVB. Acta Derm. Venereol. 85, 98–108. 10.1080/0001555051002557915823900 Bouly J.-P. Schleicher E. Dionisio-Sese M. Vandenbussche F. Van Der Straeten D. Bakrim N. . (2007). Cryptochrome blue light photoreceptors are activated through interconversion of flavin redox states. J. Biol. Chem. 282, 9383–9391. 10.1074/jbc.M60984220017237227 Chelliah V. Juty N. Ajmera I. Ali R. Dumousseau M. Glont M. . (2015). BioModels: ten-year anniversary. Nucleic Acids Res. 43, D542–D548. 10.1093/nar/gku118125414348 Clayton E. Doupé D. P. Klein A. M. Winton D. J. Simons B. D. Jones P. H. (2007). A single type of progenitor cell maintains normal epidermis. Nature 446, 185–189. 10.1038/nature0557417330052 Dai T. Gupta A. Murray C. K. Vrahas M. S. Tegos G. P. Hamblin M. R. (2012). Blue light for infectious diseases: propionibacterium acnes, Helicobacter pylori, and beyond? Drug Resist. Updat. 15, 223–236. 10.1016/j.drup.2012.07.00122846406 Eichler M. Lavi R. Shainberg A. Lubart R. (2005). Flavins are source of visible-light-induced free radical formation in cells. Lasers Surg. Med. 37, 314–319. 10.1002/lsm.2023916196041 Fischer M. R. Abel M. Lopez Kostka S. Rudolph B. Becker D. von Stebut E. (2013). Blue light irradiation suppresses dendritic cells activation in vitro. Exp. Dermatol. 22, 558–560. 10.1111/exd.1219323879817 Gandolfi A. Iannelli M. Marinoschi G. (2011). An age-structured model of epidermis growth. J. Math. Biol. 62, 111–141. 10.1007/s00285-010-0330-320177685 Grabe N. Neuber K. (2007). Simulating psoriasis by altering transit amplifying cells. Bioinformatics 23, 1309–1312. 10.1093/bioinformatics/btm04217308343 Hashmi J. T. Huang Y.-Y. Sharma S. K. Kurup D. B. De Taboada L. Carroll J. D. . (2010). Effect of pulsing in low-level light therapy. Lasers Surg. Med. 42, 450–466. 10.1002/lsm.2095020662021 Heenen M. Galand P. de Maertelaer V. Heenen P. H. (1987). Psoriasis: hyperproliferation cannot induce characteristic epidermal morphology. Cell Tissue Kinet. 20, 561–570. 10.1111/j.1365-2184.1987.tb01365.x3502926 Heenen M. Thiriar S. Noël J. C. Galand P. (1998). Ki-67 immunostaining of normal human epidermis: comparison with 3H-thymidine labelling and PCNA immunostaining. Dermatol. Basel Switz. 197, 123–126. 10.1159/0000179829732159 Hoath S. B. Leahy D. G. (2003). The organization of human epidermis: functional epidermal units and phi proportionality. J. Invest. Dermatol. 121, 1440–1446. 10.1046/j.1523-1747.2003.12606.x14675195 Kessel K. A. Habermehl D. Jäger A. Floca R. O. Zhang L. Bendl R. . (2013). Development and validation of automatic tools for interactive recurrence analysis in radiation therapy: optimization of treatment algorithms for locally advanced pancreatic cancer. Radiat. Oncol. Lond. Engl. 8:138. 10.1186/1748-717X-8-13824499557 Klein A. M. Nikolaidou-Neokosmidou V. Doupé D. P. Jones P. H. Simons B. D. (2011). Patterning as a signature of human epidermal stem cell regulation. J. R. Soc. Interface 8, 1815–1824. 10.1098/rsif.2011.024021632613 Kleinpenning M. M. Otero M. E. van Erp P. E. J. Gerritsen M. J. P. van de Kerkhof P. C. M. (2012). Efficacy of blue light vs. red light in the treatment of psoriasis: a double-blind, randomized comparative study: blue light versus red light in psoriasis. J. Eur. Acad. Dermatol. Venereol. 26, 219–225. 10.1111/j.1468-3083.2011.04039.x21435024 Kleinpenning M. M. Smits T. Frunt M. H. van Erp P. E. van de Kerkhof P. Gerritsen R. M. (2010). Clinical and histological effects of blue light on normal skin. Photodermatol. Photoimmunol. Photomed. 26, 16–21. 10.1111/j.1600-0781.2009.00474.x20070834 Lapchak P. A. Salgado K. F. Chao C. H. Zivin J. A. (2007). Transcranial near-infrared light therapy improves motor function following embolic strokes in rabbits: an extended therapeutic window study using continuous and pulse frequency delivery modes. Neuroscience 148, 907–914. 10.1016/j.neuroscience.2007.07.00217693028 Lee C.-H. Wu S.-B. Hong C.-H. Yu H.-S. Wei Y.-H. (2013). Molecular mechanisms of UV-induced apoptosis and its effects on skin residential cells: the implication in UV-based phototherapy. Int. J. Mol. Sci. 14, 6414–6435. 10.3390/ijms1403641423519108 Liebmann J. Born M. Kolb-Bachofen V. (2010). Blue-light irradiation regulates proliferation and differentiation in human skin cells. J. Invest. Dermatol. 130, 259–269. 10.1038/jid.2009.19419675580 Maari C. Viau G. Bissonnette R. (2003). Repeated exposure to blue light does not improve psoriasis. J. Am. Acad. Dermatol. 49, 55–58. 10.1067/mjd.2003.44512833008 Marino S. Hogue I. B. Ray C. J. Kirschner D. E. (2008). A methodology for performing global uncertainty and sensitivity analysis in systems biology. J. Theor. Biol. 254, 178–196. 10.1016/j.jtbi.2008.04.01118572196 Markovitsi D. (2016). UV-induced DNA damage: the role of electronic excited states. Photochem. Photobiol. 92, 45–51. 10.1111/php.1253326436855 Monfrecola G. Lembo S. Cantelli M. Ciaglia E. Scarpato L. Fabbrocini G. . (2014). The effect of visible blue light on the differentiation of dendritic cells in vitro. Biochimie 101, 252–255. 10.1016/j.biochi.2014.02.00124530862 Ng C. M. Joshi A. Dedrick R. L. Garovoy M. R. Bauer R. J. (2005). Pharmacokinetic-pharmacodynamic-efficacy analysis of efalizumab in patients with moderate to severe psoriasis. Pharm. Res. 22, 1088–1100. 10.1007/s11095-005-5642-416028009 Oh P.-S. Hwang H. Jeong H.-S. Kwon J. Kim H.-S. Kim M. . (2016). Blue light emitting diode induces apoptosis in lymphoid cells by stimulating autophagy. Int. J. Biochem. Cell Biol. 70, 13–22. 10.1016/j.biocel.2015.11.00426555754 Opländer C. Deck A. Volkmar C. M. Kirsch M. Liebmann J. Born M. . (2013). Mechanism and biological relevance of blue-light (420–453nm)-induced nonenzymatic nitric oxide generation from photolabile nitric oxide derivates in human skin in vitro and in vivo. Free Radic. Biol. Med. 65, 1363–1377. 10.1016/j.freeradbiomed.2013.09.022 Parisi R. Symmons D. P. Griffiths C. E. Ashcroft D. M. Identification Management of Psoriasis Associated ComorbidiTy (IMPACT) project team (2013). Global epidemiology of psoriasis: a systematic review of incidence and prevalence. J. Invest. Dermatol. 133, 377–385. 10.1038/jid.2012.33923014338 Pathak M. A. (1991). Ultraviolet radiation and the development of non-melanoma and melanoma skin cancer: clinical and experimental evidence. Skin Pharmacol. 4(Suppl. 1), 85–94. 1764252 Perera G. K. Di Meglio P. Nestle F. O. (2012). Psoriasis. Annu. Rev. Pathol. 7, 385–422. 10.1146/annurev-pathol-011811-13244822054142 Pfaff S. Liebmann J. Born M. Merk H. F. von Felbert V. (2015). Prospective randomized long-term study on the efficacy and safety of uv-free blue light for treating mild Psoriasis Vulgaris. Dermatology 231, 24–34. 10.1159/00043049526044167 Sadeghian K. Bocola M. Schütz M. (2008). A conclusive mechanism of the photoinduced reaction cascade in blue light using flavin photoreceptors. J. Am. Chem. Soc. 130, 12501–12513. 10.1021/ja803726a18722438 Savill N. J. (2003). Mathematical models of hierarchically structured cell populations under equilibrium with application to the epidermis. Cell Prolif. 36, 1–26. 10.1046/j.1365-2184.2003.00257.x12558657 Schneider L. A. Hinrichs R. Scharffetter-Kochanek K. (2008). Phototherapy and photochemotherapy. Clin. Dermatol. 26, 464–476. 10.1016/j.clindermatol.2007.11.00418755365 Shnitkind E. Yaping E. Geen S. Shalita A. R. Lee W.-L. (2006). Anti-inflammatory properties of narrow-band blue light. J. Drugs Dermatol. 5, 605–610. 16865864 Simonart T. Heenen M. Lejeune O. (2010). Epidermal kinetic alterations required to generate the psoriatic phenotype: a reappraisal: kinetic alterations in psoriasis. Cell Prolif. 43, 321–325. 10.1111/j.1365-2184.2010.00672.x20546247 Sparsa A. Faucher K. Sol V. Durox H. Boulinguez S. Doffoel-Hantz V. . (2010). Blue light is phototoxic for B16F10 murine melanoma and bovine endothelial cell lines by direct cytocidal effect. Anticancer Res. 30, 143–147. 20150629 Taoufik K. Mavrogonatou E. Eliades T. Papagiannoulis L. Eliades G. Kletsas D. (2008). Effect of blue light on the proliferation of human gingival fibroblasts. Dent. Mater. 24, 895–900. 10.1016/j.dental.2007.10.00618164382 van Gemert M. J. Jacques S. L. Sterenborg H. J. Star W. M. (1989). Skin optics. IEEE Trans. Biomed. Eng. 36, 1146–1154. 10.1109/10.421082606488 Wataha J. C. Lewis J. B. Lockwood P. E. Noda M. Messer R. L. Hsu S. (2008). Response of THP-1 monocytes to blue light from dental curing lights. J. Oral Rehabil. 35, 105–110. 10.1111/j.1365-2842.2007.01806.x18197843 Watt F. M. (2014). Mammalian skin cell biology: at the interface between laboratory and clinic. Science 346, 937–940. 10.1126/science.125373425414300 Weatherhead S. C. Farr P. M. Jamieson D. Hallinan J. S. Lloyd J. J. Wipat A. . (2011). Keratinocyte apoptosis in epidermal remodeling and clearance of psoriasis induced by UV radiation. J. Invest. Dermatol. 131, 1916–1926. 10.1038/jid.2011.13421614017 Weatherhead S. C. Farr P. M. Reynolds N. J. (2013). Spectral effects of UV on psoriasis. Photochem. Photobiol. Sci. 12, 47–53. 10.1039/C2PP25116G23023652 Weinstabl A. Hoff-Lesch S. Merk H. F. von Felbert V. (2011). Prospective randomized study on the efficacy of blue light in the treatment of psoriasis vulgaris. Dermatology 223, 251–259. 10.1159/00033336422105015 Weinstein G. D. McCullough J. L. Ross P. A. (1985). Cell kinetic basis for pathophysiology of psoriasis. J. Invest. Dermatol. 85, 579–583. 10.1111/1523-1747.ep122835944067329 Weinstein G. D. Van Scott E. J. (1965). Autoradiographic analysis of turnover times of normal and psoriatic epidermis. J. Invest. Dermatol. 45, 257–262. 10.1038/jid.1965.1265837896 Yoshida A. Yoshino F. Makita T. Maehata Y. Higashi K. Miyamoto C. . (2013). Reactive oxygen species production in mitochondria of human gingival fibroblast induced by blue light irradiation. J. Photochem. Photobiol. B Biol. 129, 1–5. 10.1016/j.jphotobiol.2013.09.00324141287 Zhang H. Hou W. Henrot L. Schnebert S. Dumas M. Heusèle C. . (2015). Modelling epidermis homoeostasis and psoriasis pathogenesis. J. R. Soc. Interface 12. 10.1098/rsif.2014.107125566881 Zhou H. Hu C. Zhu Y. Lu M. Liao S. Yeilding N. . (2010). Population-based exposure-efficacy modeling of ustekinumab in patients with moderate to severe plaque psoriasis. J. Clin. Pharmacol. 50, 257–267. 10.1177/009127000934369519934030 Zonios G. Bykowski J. Kollias N. (2001). Skin melanin, hemoglobin, and light scattering properties can be quantitatively assessed in vivo using diffuse reflectance spectroscopy. J. Invest. Dermatol. 117, 1452–1457. 10.1046/j.0022-202x.2001.01577.x11886508 Abbreviations BL

Blue Light

 Pv

Psoriasis vulgaris

 F

Fluence

 UV

Ultraviolet light

 BLISS

Computational model for Blue Light Irradiation of Psoriatic Skin

 Psc

Stem Cells

 Pta

Transit Amplifying cells

 Pga

Growth Arrested cells

 Psp

Spinous cells

 Pgc

Granular Cells

 Pcc

Corneocytes

 Iav

Average power density

 CW

Continuous Wave mode

 PW

Pulsed Wave mode

 Ip

Peak power density

 PASI

Psoriasis Area and Severity Index

 LPSI

Local Psoriasis Severity Index

 MPSA

Multiple parameter sensitivity analysis

 LPSA

Local sensitivity analysis

 LI

Low Intensity

 HI

High Intensity.

1https://github.com/ZFelixGarza/BLISS

2http://www.ebi.ac.uk/biomodels-main/