Posted on in the International Journal of Philippine Science and Technology 2015 Volume 8 Issue Number 1 Page 52—57
*Corresponding Author
Email Address: pcnaval@dcs.upd.edu.ph
Submitted: August 24, 2014
Revised: February 6, 2015
Accepted: March 1, 2015
Published: May 30, 2015
Abstract—Ultraviolet radiation (UVR) is proven to be an effective treatment for psoriasis yet it needs to be carefully controlled since it is a widely known carcinogen. Due to the time-consuming, expensive and hazardous nature of clinical tests on actual patients, the range of possible values for the optimal regimen set has not yet been thoroughly explored. An agent-based model was constructed in order to study the complex role of UVR in psoriasis phototherapy, including possible contribution to the formation of skin cancers. The treatment progression of the simulated psoriasis undergoing UVR phototherapy was monitored weekly and compared with the corresponding clinical data. For model calibration, a single-objective optimization via genetic algorithm was employed that aimed to minimize the discrepancy between the model output and the expected clinical result. After fitting and validation, the model was then subjected to multi-objective evolutionary optimization using Non-dominated Sorting Genetic Algorithm-II (NSGA-II) in order to suggest sets of optimal UVR phototherapy regimen for psoriasis taking into account the safety, clearance time as well as aggressiveness of the therapy. Results show a good model fit against clinical data and gave some plausible sets of dosimetry for UVR phototherapy that maximize its therapeutic efficacy while minimizing the associated skin cancer risk.
Keywords—psoriasis, phototherapy optimization, skin model, epidermal kinetics, uvb exposure
Psoriasis is a chronic inflammatory skin disorder generally considered an immune-mediated disease having a prevalence rate of around 2-3% worldwide (Roy et al 2010). Clinically, it manifests itself as scaling, redness, and hardening of the skin (Schon and Boehncke 2005). Although its pathogenesis is not yet fully understood, cellular features such as hyperproliferation of keratinocytes, infiltration of immune cells and overproduction of pro-inflammatory cytokines are believed to be the crucial factors for its clinical manifestation (Gudjonsson et al 2004).
One of the most common treatments for psoriasis is Ultraviolet Radiation (UVR) phototherapy. In this therapy, the psoriatic skin is exposed to UVR whose immunosuppressive effects include the induction of apoptosis on infiltrating immune cells and the production of anti-inflammatory cytokines (Ozawa et al 1999; Weichenthal and Schwarz 2005). Although UVR is proven effective, its use requires delicate balancing since it is a widely known carcinogen (Soehnge et al 1997). UVR induces DNA damage not only to the infiltrating immune cells but also to the normal skin cells. Hence, although such DNA damage causes the elimination of the immune cells via apoptosis, it may also result in epidermal cell mutations which may lead to skin cancers as reported from studies involving ultraviolet-A wavelength phototherapy (Nijsten and Stern 2003).
For this study, the scope is on the use of the more popular ultraviolet-B (UVB) wavelength. Although no evidence exists regarding UVB phototherapy-induced skin cancers (Lee et al 2005), UVB radiation is also known to be highly mutagenic (Kappes et al 2006). To the best of our knowledge, attempts to find the optimal dosimetry for psoriasis UVB phototherapy either by clinical experiment (Asawanonda et al 2000; Boztepe et al 2006; Hallaji et al 2010) or by modeling (Diffey 2004) have not yet obtained a conclusive answer on what set of regimen values maximize the therapeutic effect of UVR and minimize its side effects. In this study, agent-based modeling was employed by representing the interaction of epidermal and immune cells during psoriasis as well as their responses under UVR exposure especially the possible induction of precancerous keratinocyte development. Using such in silico representation, multi-objective optimization of phototherapy regimen was done where sets of values for initial UVR dose, percentage increment, frequency of exposure along with the dose of adjunct topical glucocorticoids (GC) were given.
Extensive literature search and consultations with an expert dermatologist were done in order to identify the crucial factors for psoriasis pathogenesis. A representative 0.35 mm2 cross sectional area of the whole depth of epidermis was modeled. The main cell types represented include keratinocytes, melanocytes, Langerhans cells, cytotoxic T cells (CTL or Tc cells), T helper 1 cells, macrophages, and neutrophils. Intercellular communication was simulated using three types of signals encompassing a broad array of cytokines and chemokines similar to the classification used in Basic Immune Simulator (BIS) (Folcik et al 2007). These include the pro-inflammatory cytokines produced by T cells (CK1 or cytokine-1), and the pro- (PK or parenchymal-kine) and anti-inflammatory cytokines (CK2 or cytokine-2) produced by macrophages and keratinocytes. Similarly, the reception of signals outside the epidermis was abstracted via the so-called portal agents.
Clinicians use the Scaling-Erythema-Induration (SEI) scoring for disease severity measurements which ranks the three psoriatic skin features from 0 (none) to 4 (worst) and then adding the ranks to attain a value from 0-12. Since the overall SEI score is assumed to correlate linearly with the number of infiltrating immune cells (Ozawa et al 1999; Gudjonsson et al 2004), the count of such cells was used as an indication of psoriasis severity in the model.
For UVR measurements, the unit used by clinicians is usually in terms of the so-called Minimum Erythema Dose (MED) in which 1 unit in MED is the minimum amount of UVR needed to induce redness to the skin. In the model, 1 MED UVR exposure inflicts a particular amount of DNA damage to immune cells and keratinocytes giving them a certain level of probability to release cytokines, undergo apoptosis or become precancerous.
An agent-based model written in Java was constructed with the MASON multiagent simulation framework (Luke et al 2004). The framework already provides the base classes needed to create a 3D system. For the stationary keratinocytes and melanocytes composing more than 90% of the cells, the ObjectGrid3D class of the framework was used. For the mobile immune cells, the SparseGrid3D was utilized. The Diffuser implementation included in the framework was used to model the diffusion of signals which were represented using DoubleGrid3D class. One time step is equivalent to 1 hour.
Figure 1. Overview of agent interactions.
Rules governing the behavior of the individual agents can be either qualitative or quantitative. The qualitative rules define the interactions of the different agents shown in Fig. 1. Quantitative rules are usually computed intracellularly as defined in Equations 1 to 8. Table 1 summarizes both the qualitative and quantitative rules used.
The model equations that define the agents’ behavior are the following:
(1)
where U is the value of photoadapted quantity of UVR as a function of the raw UVR dose U_{raw} and the number of exposures X with constants a and b (see Table 2). The equation reliably describes the findings from (Palmer et al 2006; Gange et al 1985; de Winter et al 2001).
(2)
where the value of base DNA damage D_{base} is a function of effective UVR dose U and the DNA damage D_{med} which is the theoretical damage that could be attained at a UVR dose of 1 MED. The equation is based on the dose-response study from (Young et al 1998; de Leeuw et al 1995).
(3)
where the newly inflicted DNA damage D_{new} is a function of the base DNA damage D_{base} of a cell, the skin layer L where it resides and the skin type S of the person. The trends from (Kobayashi et al 1998; Tadokoro et al 2003; Anderson and Parrish 1981) are reliably simplified by the equation.
(4)
where the DNA damage D_{t} at time t after irradiation depends on the sum of the unrepaired DNA damage D_{unrep} and the newly inflicted DNA damage D_{new} with h as the DNA damage half-life constant. This decay equation along with the half-life constant was approximated from (Errico et al 2003; Young et al 1996).
(5)
where the probability P_{t} for a cell to become precancerous depends on the amount of DNA damage D_{t} at time t after irradiation multiplied by a constant c subject to adjustments. This equation was assumed from (Jonason et al 1996).
(6)
where A_{max} is the maximum probability for apoptosis, D_{t} is the DNA damage, H_{a} is the hill factor, M_{a} is the DNA damage that attains half of A_{max}. The trends from (Bertrand-Vallery et al 2010; Teunissen et al 1993) were fitted into a sigmoid equation known as the E_{max} equation (MacDougall 2006). The equation has few parameters which can be easily interpreted clinically. Also, the values of these parameters can be reliably estimated from literature findings or can be parameterized to fit a particular observation. Such equation is also used in Items 7 and 8.
(7)
where E_{max} is the maximum probability to induce the two GC effects mentioned, G is the GC dose, H_{e} is the hill factor, M_{e} is the GC dose that attains half of E_{max}. This is based from the sigmoidal patterns observed by (Jamieson and Yamamoto 2000; Madretsma et al 1996).
(8)
where R_{max} is the maximum probability for cytokine release, D_{t} is the DNA damage at time t, H_{r} is the hill factor, M_{r} is the amount of DNA damage that achieves half of R_{max}. The findings from (Barr et al 1999; Young et al 1998b) were reliably described by this equation.
where R_{max} is the maximum probability for cytokine release, D_{t} is the DNA damage at time t, H_{r} is the hill factor, M_{r} is the amount of DNA damage that achieves half of R_{max}. The findings from (Barr et al 1999; Young et al 1998b) were reliably described by this equation.
Table 1. Summary of agent behavior and interactions.
Agent Behavior^{§} | Figures and Equations |
KERATINOCYTE | |
Assess DNA damage | Eqns 1 to 3 |
Detect/release cytokine | Eqn 8, Fig 1 edges 2, 5, 21 |
Differentiate into a nonviable cell | |
Undergo apoptosis | Eqn 6 |
Undergo DNA repair | Eqn 4 |
Acquire precancerous behavior | Eqn 5 |
MELANOCYTE | |
Maintain melanin | Fig 1 edge 1 |
LANGERHANS CELL | |
Detect Cytokine | Fig 1 edges 14, 15, 18 and 20 |
Present antigen to Th1 cell | Fig 1 edge 4 |
Activate Th1 cell | Fig 1 edge 8 |
Recruit CTL and Th1 cells | Fig 1 edges 6 and 7 |
Tc CELL | |
Assess DNA damage | Eqns 1 to 3 |
Eliminate autoantigenic keratinocyte | |
Undergo apoptosis | Eqn 6 |
Undergo DNA repair | Eqn 4 |
Move | |
Th1 CELL | |
Assess DNA damage | Eqns 1 to 3 |
Release cytokine | Fig 1 edge 13 |
Undergo apoptosis | Eqn 6 |
Undergo DNA repair | Eqn 4 |
Move | |
MACROPHAGE | |
Assess DNA damage | Eqns 1 to 3 |
Release cytokine | Fig 1 edge 12 |
Present antigen to Th1 cell | Fig 1 edge 3 |
Activate Th1 cells | Fig 1 edge 11 |
Undergo apoptosis | Eqn 6 |
Undergo DNA repair | Eqn 4 |
Move | |
NEUTROPHIL | |
Assess DNA damage | Eqns 1 to 3 |
Eliminate autoantigenic keratinocyte | |
Undergo apoptosis | Eqn 6 |
Undergo DNA repair | Eqn 4 |
Move | |
PORTAL AGENT | |
Detect cytokine | Fig 1 edges 16, 17, 19 and 22 |
Recruit macrophages and neutrophils | Fig 1 edges 9 and 10 |
Move |
The model was initialized such that there is already a population of immune cells that have infiltrated the modeled epidermis. By approximating the results from (Szabo et al 1998; Hoath and Leahy 2003; Bauer et al 2001; Fuentes-Duculan et al 2010), the corresponding relative numbers of the epidermal and immune cells for a psoriasis severity of SEI 8.0 are given in Table 2.
Table 2. Values of equation constants and model parameters.
Constants and Parameters | Value | References* |
A | 0.45 | de Winter et al 2001; Palmer et al 2006; Gange et al 1985; |
B | 5 | |
H | 20.65 | |
c | 13.33 | Jonason et al 1996 |
H_{a} (keratinocyte) | 1.468 | Bertrand-Vallery et al 2010 |
M_{a} (keratinocyte) | 750 | |
H_{e} | 0.830 | Jamieson and Yamamoto 2000; Madretsma et al 1996 |
M_{e} | 0.005 | |
H_{r}(CK2) | 1.0 | Arbitrary hill factor |
H_{r}(PK) | 3.17244 | Barr et al 1999; Young et al 1998b |
H_{a}(immune cell) | 1.192 | Teunissen et al 1993 |
M_{a}(immune cell) | 8319.93^{§} | Automated tweaking |
M_{r}(PK) | 124.33^{§} | Automated tweaking |
M_{r}(CK2) | 68.36^{§} | Automated tweaking |
A_{max} | 100 | Arbitrary damage limit |
E_{max} | 100 | Probability limit |
R_{max} | 100 | |
Minimum CK1 to trigger an effect | 49.43^{§} | Automated tweaking |
Minimum CK2 to suppress CK1 effect | ||
Minimum PK to trigger an effect | 87.24^{§} | Automated tweaking |
Standard GC dose | 0.002421^{§} | Automated tweaking |
Initial number of keratinocytes | 9375 | Szabo et al 1998; Hoath and Leahy 2003; Bauer et al 2001; Fuentes-Duculan et al 2010 |
Initial number of melanocytes | 625 | |
Initial number of Tc cells | 46 | |
Initial number of Th1 cells | 51 | |
Initial number of macrophages | 44 | |
Initial number of neutrophils | 51 | Arbitrary |
Initial number of portal agents | 30 |
The iterative interaction of the genetic algorithm and the agent-based model is summarized in Figure 2a. Both single-objective and multi-objective optimization experiments were done for this study using single- and multi-objective genetic algorithm (NSGA-II), respectively.
Figure 2. (a) Overview of the optimization process. (b) The Pareto front (black dots) among the candidate solutions for a minimization experiment.
In a single-objective optimization, only one criterion (i.e. objective) is used for the evaluation of the candidate decision variables during each iteration. For multi-objective optimization, the concept of Pareto dominance is used in which a candidate solution is said to dominate another solution if and only if it is equally good or better in all objectives and that there exists one objective where it is absolutely better. The sets of candidate solutions that are not dominated, known as the Pareto front (Figure 2b), are highly favored in the multi-objective optimization process.
The model output is in terms of the total count of immune cells over time along with the count of severely damaged (precancerous) keratinocytes. In order to produce realistic trends, it was necessary to first fit the model outputs to the clinical ones (i.e., as training or calibration data set). The findings from (Dawe et al 2002; Asawanonda et al 2005; Jamkhedkar et al 1996; Taneja et al 2003) are used as basis to compute the weekly difference of the clinical output and the model output which serves as the model error. This data set consists of a total of four plots showing the progression of psoriasis under varying treatment scenarios namely a) 0.70 MED initial UVR, 20% increment, 3X per week exposure (also referred to as the standard regimen) b) constant 4.0 MED UVR, 3X per week exposure c) topical application of glucocorticoid at standard dose (no UVR) d) no treatment.
Using a single-objective genetic algorithm (Deb et al 2002), the model error computed was then subjected to minimization by the automated tweaking of six model parameters identified in Table 2. More formally, the minimization problem can be stated as follows:
(9)
where
For fitted model validation, test runs involving clinical results not used during calibration were performed. Such test data set is composed of other results from (Asawanonda et al 2005; Dawe et al 2002) which similarly show the trends in psoriasis clearance upon exposure to specific dosage of UVR.
Lastly, the model was also subjected to qualitative validation by observing the formation of the so-called precancerous keratinocytes. For comparison, the results from (Jonason et al 1996) were used which roughly describe the qualitative aspects of the formation of precancerous keratinocytes upon exposure to varying intensities of sunlight.
Figure 3. Common treatment options often encountered in clinical practice.
After the model showed comparable trends with clinical results, it was then subjected to predictive experiments via multi-objective evolutionary optimization using NSGA-II (Deb et al 2002). The goal is to find the value of initial UVR dose, percentage increment, frequency of exposure along with the dose of adjunct topical GC which will satisfy the criteria for each case shown in Figure 3 namely aggressive UVB-only treatment without time limit (Case 1), mild UVB-only treatment lasting 8 weeks (Case 2) or 4 weeks (Case 3) and combined UVB and topical GC treatment lasting 8 weeks (Case 4) or 4 weeks (Case 5). The treatment options identified in the figure are clinically among the most common options chosen by patients for psoriasis therapy. To perform the optimization, treatment variables and objectives for optimization were set in the manner shown in Table 3.
Table 3. Treatment cases and their associated objectives and variables.
Category | Case 1 | Case 2 | Case 3 | Case 4 | Case 5 |
A) Objectives | |||||
Precancerous Cell Count | f_{1} | f_{1} | f_{1} | f_{1} | f_{1} |
Clearance Time | f_{2} | ||||
Relative UVR Dose | f_{2} | f_{2} | f_{2} | f_{2} | |
Relative GC Dose | f_{3} | f_{3} | |||
(auxiliary objective) | f_{3} | f_{3} | f_{3} | f_{4} | f_{4} |
B) Treatment Variables | |||||
UVR Initial Dose | ✓ | ✓ | ✓ | ✓ | ✓ |
Percentage Increment | ✓ | ✓ | ✓ | ✓ | ✓ |
Exposure Frequency | ✓ | ✓ | ✓ | ✓ | ✓ |
GC Dose | ✓ | ✓ |
Mathematically, the multi-objective optimization problem can be stated as follows:
(10)
where
Figure 4a shows the result of the automated model fitting using single objective GA. The lowest overall model error found by the genetic algorithm is 0.39 SEI score difference. Hence, the values for the six selected model parameters (i.e. the decision variables) that produced such small discrepancy were obtained and used for the next series of runs.
Figure 4. Trends for psoriasis clearance shown by the model compared with experimental data during (a) fitting or calibration and (b) validation or testing involving 10 trials. Qualitative assessment of model performance was also done by taking note of the cluster formation of precancerous keratinocytes (c and d).
The accuracy of the trends for psoriasis clearance shown by the model was further verified after comparing the model output with another batch of clinical results involving different sets of regimens which were not used during calibration. Figure 4b shows a relatively good fit from such test plots involving 10 runs.
Interesting patterns of cluster formation also emerged after the qualitative assessment of the clusters of developing precancerous keratinocytes. In general, as the UVR treatment becomes more aggressive, the number of large clusters formed increases. Figure 4c shows the average number of precancerous clusters reaching a particular size after 10 simulation runs at various UVR dosages. An example of a precancerous cell cluster in the model is shown in Figure 4d.
Lastly, predictions from the model showed some potential optimal phototherapy regimens that satisfy the treatment conditions and constraints set for various cases. The Pareto optimal solution set returned by NSGA-II contains numerous solutions for the values of the regimen variables and each of these solutions has its own priority objectives. However, clinically, the safest (i.e. have the least no. of precancerous cells) regimen values are more relevant. Hence, for three different trials involving 25 generations of NSGA-II each, the top three regimen sets with the minimum risk for skin cancer were chosen. The suggested regimens for all cases are shown in Figure 5a while the associated objective values for each case are plotted in Figure 5b to d. It can be seen that for Case 1, the suggested UV initial is twice that of the standard regimen while the suggested percent increment and frequency are very close to the standard values. The treatment regimen suggested for Case 2 is almost similar to that of Case 1. Interestingly, for Case 3, the suggested UV initial is five times that of the standard regimen and the values for percent increment and frequency are twice the standard values. For Case 4 and 5 which involve the combination of UV and topical GC, the UV treatment regimen they both suggested are similar and so much less than the standard values. From a clinical perspective, the UV treatment suggested for both cases are almost insignificant. Nevertheless, these two cases differ on the value of suggested topical GC dose. Case 5, having shorter treatment period, required twice the standard GC dose while Case 4 required a GC treatment close to the standard amount.
Figure 5. (a) Mean optimal values of the treatment variables (obtained from 3 NSGA-II runs) along with the resulting objective values for (b) Case 1, (c) Case 2 and 3, and (d) Case 4 and 5.
An agent-based model was constructed based on the recent knowledge of psoriasis pathogenesis as obtained from extensive literature search as well as numerous consultations with an expert dermatologist. The constructed model was able to demonstrate some of the trends and properties of the biological system being modeled. For instance, the single-objective genetic algorithm-guided model calibration satisfactorily fitted the clinical data with a discrepancy that would be very difficult to detect using a visual inspection-based clinical method. Comparison with the test data set ensured that even the trends not used during calibration are also replicable with less error hence minimizes the possibility of overfitting. With these tests, the ability of the model to predict the progression of psoriasis clearance using varying amounts of UVR exposures has been validated.
Another important biological pattern which was captured in the model is the pattern of precancerous keratinocyte cluster formation. The resulting precancerous cell cluster formation in the model somehow replicates the finding from (Jonason et al 1996) where larger clusters and higher number of precancerous cells are formed under aggressive sunlight exposures which also have UV wavelengths. This result also provides support for the hypothesis that UVB-induced apoptosis serves two contrasting roles on skin tumor development (Zhang et al 2005). According to their findings, apoptosis causes the elimination of highly damaged cells but in doing so, it leaves vacant spaces within the tissue which favors the clonal expansion of already existing mutant cells. Such dynamics between mutant cell elimination and clonal expansion by apoptosis can be further explored in the model.
To verify the plausibility of the suggested values for the phototherapy regimen set generated by multi-objective optimization, consultation with an expert dermatologist was done.
For Case 1, the resulting regimens are all slightly more aggressive than the commonly used standard regimen. This result seems plausible when the objectives are taken into account. The objective involving faster time for clearance may have caused the regimen values to be more aggressive than the standard regimen values since higher UVR dosages could result in faster elimination of the infiltrating immune cells as shown from (Asawanonda et al 2000). However, due to the safer-treatment objective, the values were prevented from becoming too high since that would cause higher incidence of precancerous cell formation.
For Case 2 and 3, a time constraint of 8 weeks and 4 weeks were given respectively. Also, the objective involving faster clearance was removed and replaced by an objective that minimizes the total dosage of UVR to be given. Hence, this means that the values found for these two cases are the mildest possible dosages that would ensure disease clearance within the specified period. It seems that having a very short treatment period (i.e. Case 3 involving a period of just four weeks) needs a very high UV treatment regimen. On the other hand, having a period of 8 weeks (i.e. Case 2) did not cause a significant increase in the values for UV treatment. In fact, the regimen values suggested for Case 2 are identical to that of Case 1 where no restriction of treatment period was set at all aside from the fact that faster clearance time was set as an objective. These results have an important impact when it comes to prescribing an appropriate dosage given the time preference of the patients. For instance, patients who opt to have a shorter time for therapy may need to undergo additional precautionary evaluation for their tolerance to endure higher UV dosages.
Lastly, for Case 4 and 5, the scenarios where UVR exposure combined with adjunct topical glucocorticoid were analyzed. In these experiments, the least possible amounts of both UVR and topical glucocorticoid to attain clearance were obtained for both 8-week and 4-week treatment session, respectively. The results, however, returned almost negligible UVR dosage. Hence, the regimen sets suggested are almost entirely glucocorticoid-based. Due to the objective that minimizes the glucocorticoid dose, we can say that the values for glucocorticoid dose given are the least possible values that would ensure clearance at the given time period. The required dosage for Case 4 is close to the standard amount while the required dosage for Case 5 is around twice as much due to the fact that a shorter treatment period was set for the latter.
The constructed genetic algorithm-driven agent-based model of epidermis undergoing psoriasis phototherapy can be an excellent tool for the exploration of two main biological processes, namely the progression of psoriasis treatment and the consequent initiation and promotion of precancerous cell cluster formation caused by UVR exposure. With respect to the former, it was shown that the model is accurate enough to provide dermatologically feasible trends which could be good starting grounds for further clinical testing. With respect to the latter, the observed pattern of precancerous cell formation caused by UVR and the contrasting effect of apoptosis provide interesting areas which could be further developed and studied. As the amount of dermatological data increases every year, more details can still be incorporated into the constructed model that would increase its validity in order to provide a bigger picture of the processes occurring during psoriasis treatment including skin cancer development. For instance, the increase in knowledge regarding the complex role of the cytokine and chemokine network may need to be incorporated in detail in the model to get a more accurate representation of the disease. Also, the recent findings from studies making use of modern methods and tools for disease severity assessment will greatly help in minimizing the subjective aspect involved in the study.
The first author was supported by DOST-SEI. The authors are indebted to Aletta T. Yñiquez and Marianne G. Camoying for helpful discussions and feedback.
None.
All the authors jointly conceived the study. A.C.C. and F.C.R. developed the skin model. A.C.C. with P.C.N. designed the simulation experiments while E.R.M. contributed extensively to the analysis of results. Both P.C.N and E.R.M. supervised the project. All authors discussed results and commented on the manuscript at all stages.
STRIDE Grants The next call for Grants will commence on February 27, 2015. To know more about the Grants program of STRIDE, you may visit our website at stride.org.ph.