## Abstract

We present a computational model to study the spatio-temporal dynamics of epidermis homoeostasis under normal and pathological conditions. The model consists of a population kinetics model of the central transition pathway of keratinocyte proliferation, differentiation and loss and an agent-based model that propagates cell movements and generates the stratified epidermis. The model recapitulates observed homoeostatic cell density distribution, the epidermal turnover time and the multilayered tissue structure. We extend the model to study the onset, recurrence and phototherapy-induced remission of psoriasis. The model considers psoriasis as a parallel homoeostasis of normal and psoriatic keratinocytes originated from a shared stem cell (SC) niche environment and predicts two homoeostatic modes of psoriasis: a disease mode and a quiescent mode. Interconversion between the two modes can be controlled by interactions between psoriatic SCs and the immune system and by normal and psoriatic SCs competing for growth niches. The prediction of a quiescent state potentially explains the efficacy of multi-episode UVB irradiation therapy and recurrence of psoriasis plaques, which can further guide designs of therapeutics that specifically target the immune system and/or the keratinocytes.

## 1. Introduction

The epidermis, the outermost layer of skin, provides the human body a physiological barrier to the environment and protects the body from water loss, pathogenic infection and physical injury. The epidermis organizes into a stratified structure of keratinocytes at several differentiated stages [1], which constitute 95% cell population in the epidermis [2]. Like other regenerative tissues, the epidermis constantly renews itself to replace desquamated and apoptotic keratinocytes, repair tissue damage and establish homoeostasis. The renewal is orchestrated by a cascade of cellular processes including proliferation, differentiation, migration, apoptosis and desquamation [3–5]. A keratinocyte transits spatially from the stratum basale to the stratum corneum during its lifespan and meanwhile experiences multi-stage biochemical and morphological changes. Many endogenous and exogenous factors (e.g. Ca^{2+} concentration, cytokines, UV irradiation, etc.) affect epidermal dynamics and the homoeostasis of the epidermis by regulating one or more cellular processes.

Mathematical and computational models have long been useful tools to predict cellular behaviours of epidermis renewal under normal or pathological conditions. Previous models for epidermal dynamics have usually adopted two approaches. One approach includes deterministic models that derived analytical solutions to stationary cell populations. For example, the model by Savill [6] described proliferation of stem cells (SCs) and transit-amplifying (TA) cells and differentiation to post-mitotic cells, which predicted the influences of apoptosis, cell-cycle time and transit time on cell populations. Gandolfi *et al*. [7] proposed a spatio-temporal model to investigate the evolution of the epidermis, which described cell motion by a constitutive equation. The other approach includes agent-based models that treat individual keratinocytes as computing entities operating under specific physical and biological rules. Such models can simulate the multi-layer epidermal structure organized by cell proliferation, differentiation, death and migration, in which non-specific intracellular and extracellular biochemical factors affected cell proliferation and differentiation, while physical adhesive and repulsive forces governed cell motion [8–12].

In this paper, we present a hybrid model to combine advantages of the above two approaches. The model uses a mean-field cell population kinetics together with an agent-based model for cell migration. The model computes population dynamics of epidermal renewal without having to compute the cell movements simultaneously, allowing fast and analytical evaluations of modelling hypotheses and results. The population kinetics model describes cellular processes including cell division, differentiation, apoptosis and desquamation. Either deterministic or stochastic simulation can be used to generate the population dynamics of keratinocytes. Cell migration is described by a two-dimensional agent-based model that tracks the cell movement driven by cell–cell interactions. Simulation of cell migration can be integrated with the stochastic population dynamics to visualize tissue stratification and establishment of homoeostasis. A properly parametrized model reproduces experimentally observed epidermis growth, differentiation and desquamation dynamics, homoeostatic density distribution over different types of keratinocytes and the epidermis turnover times of different cell compartments.

To investigate an important pathological condition of the skin, we study the onset and recurrence of psoriasis plaques and their recovery under phototherapy using UVB irradiation. Psoriasis is a complex epidermal disorder characterized by keratinocyte hyperproliferation and abnormal differentiation due to intricate interactions with the immune system. The disease affects 2–4% of the general population and currently has no cure [13,14]. Specifically, we hypothesize a novel mechanism of SC–immune system interaction and predict the chronic disorder as a bimodal switch between a disease phenotypic and a quiescent (seemingly normal) state. The model hypothesizes a parallel epidermal homoeostasis simultaneously maintained by the normal and the psoriatic keratinocytes. The psoriatic homoeostasis is caused by permanent perturbations in cell division, apoptosis and differentiation, derived from defective SCs and their interactions with a weakened immune system. For treatment, the model predicts that to achieve a controlled remission, the effective treatment of UVB irradiation must reduce the high-density psoriatic epidermis below a threshold level by activating apoptosis in SCs, potentially explaining the chronicity and recurrences of the disorder and providing a guide to design feasible therapeutics.

## 2. The model

The model consists of (i) a kinetics model, which tracks the temporal evolution of the cell population of keratinocytes at several differentiation stages, and (ii) a migration model, which describes the motion of individual cells to generate the stratified structure of the epidermis. One can choose the kinetic model as a standalone module to compute the mean-field population dynamics or can integrate the two models to visualize the renewal kinetics and stratification of the tissue.

### 2.1. Model of cell population kinetics

Figure 1 illustrates the *central transition pathway* of epidermis renewal. The pathway considers the population dynamics of six categories of keratinocytes stratified from the stratum basale to the stratum corneum, including progenitors, SCs and TA cells in the proliferative compartment, and differentiated cells in the non-proliferative compartment with growth-arrested (GA) cells, spinous (SP) cells, granular cells (GC) and corneocytes (CC). The above classification is primarily based on known histological structure and physiological function of the human epidermis [15]. GA cells are precursors for non-proliferating cells. SP cells and GC are fully differentiated keratinocytes, and the non-nucleated CC represent the end-stage differentiation and eventually desquamate. The central transition pathway incorporates three main cellular processes: (i) proliferation of SCs and TA cells, (ii) differentiation including several inter-category cell conversions, and (iii) cell loss including apoptosis of nucleated cells and desquamation of the CC.

How SCs maintain the epidermal homoeostasis remains unresolved [16]. However, recent long-term *in vivo* lineage tracing studies on mouse tail tissue by genetic labelling revealed remarkable details about clonal dynamics of epidermal progenitors, suggesting the existence of either single (SC alone) or two progenitors (SC and committed progenitor) [17–19]. As illustrated in figure 2, our model considers a slow cycling SC population together with a faster proliferating committed progenitors (or, TA cells). A SC divides in one of three modes [20–22]: (i) self-proliferation, by which a SC divides into two daughter SCs, (ii) asymmetric division, by which a SC divides into a SC and a TA cell, or (iii) symmetric division, by which a SC divides into two TA cells. Considering a finite availability of SC niches [20,23], we assume a logistic growth of SCs to limit the SC density by a maximal growth capacity, which ensures the system to reach a well-defined steady state (see [24] for a more general model that guarantees a steady state). Previous models [17,25] required a precise balance between SC self-proliferation and symmetric division and were intolerable to arbitrary perturbations such as population random drift caused by intrinsic stochasticity in the three-mode SC division. Similarly, TA cells also divide in one of the three modes of self-proliferation, symmetric and asymmetric division into the GA cells [17,26,27]. In addition, a TA cell may resume the SC state and a GA cell may resume a TA cell state by backconversions [28].

The rate of progenitor division is often characterized by the cell-cycle time and by the subpopulation of cells that are actively dividing (also known as the ‘growth fraction’). Our mean-field model does not distinguish proliferative propensity in individual cells and therefore parametrizes cell divisions with empirical rate constants that integrate influences of the cell cycle and the growth fraction. Environmental changes regulate the proliferation rate of SCs. For example, the need to repair tissue damage promotes SC proliferation [19,21,29]. A recent study of hair follicles showed that TA cells may signal to SCs to regulate proliferation [30]. To incorporate this feedback mechanism, we assume an empirical dependence of SC division rate constants, *γ*_{1}, *k*_{1a} and *k*_{1s}, on the density of TA cells and define2.1where the subscript ‘h’ indicates a homoeostatic rate constant (see table 3 for numerical values), and *ω* ≡ *r _{x}*

_{,max}/

*r*

_{x}_{,h}is the ratio of the maximum division rate to the homoeostatic rate and is assumed identical for all division processes.

*ω*reflects the maximum increase in the growth fraction and/or decrease in the cell-cycle time when SC proliferation accelerates. At homoeostasis of the normal epidermis, the reported growth fraction varied from 20 to 70% [31,32]. A study in mice epidermis found a more than 10-fold decrease in the cell-cycle time from 5–7 days to 11 h following tissue abrasion [33], whereas no significant change in cell-cycle time was found in psoriasis [34,35]. The exponent

*n*models the sensitivity to the deviation of TA cell density from the homoeostasis. In the limit of a much reduced TA cell density (), SCs divide at the maximum rate,

*r*

_{x}_{,max}, whereas SCs divide at a minimal rate when TA cells overpopulate (). The total SC proliferating rate at homoeostasis is set about 10

^{−2}per day, aligned with four to six division events per year [19].

The non-proliferative compartment describes a cascade of differentiations from GA cells to CC. The model also considers an apoptosis process for all nucleated keratinocytes. Early studies suggested that apoptosis was only significant for SCs and TA cells in the proliferative compartment [36]. Recent experiment [37] also showed that apoptosis was evident in the differentiated keratinocytes. The extent of apoptosis is commonly characterized by the *apoptotic index*, which is typically quantified using the transferase-mediated uridine nick end labelling assay [36]. For a cell type *i*, the apoptosis index equals the probability that a cell undergoes apoptosis, defined as the ratio of apoptosis rate to the total outflux including rates of apoptosis and transition to the downstream cell category2.2For SCs and TA cells, *k*_{1} = *k*_{1s} and *k*_{2} = *k*_{2s} are symmetric division rate constants. Experiments often reported global apoptosis indices that did not differentiate apoptotic activities in proliferative and non-proliferative compartments. Apoptotic heterogeneities of cells within a compartment are even less known. Here, each cell category is assumed an identical apoptosis index that is used to calculate the apoptosis rate constants .

Governing system equations (ODEs) are listed in table 1 and computation of the model is described in the electronic supplementary material, figure S1.

### 2.2. Model of cell migration

The agent-based migration model describes movement of all keratinocytes in a two-dimensional (cross-sectional) epidermis volume. Keratinocytes once derived from SCs move outward from the stratum basale to the outermost stratum corneum, to compose a stratified epidermis. The model describes cell mechanics that propels cell movement. An individual cell is subject to three forces: (i) a viscous force due to cells moving in the surrounding environment; (ii) a repulsive force due to cell–cell compression, and (iii) an adhesive force due to interactions among adhesive molecules on cell membranes. Considering the sluggish keratinocyte motion (on a scale of µm h^{−1}) in a fluidic environment with a low Reynolds number (i.e. viscosity dominates inertia) [38], the model neglects the acceleration due to inertia. The model also considers keratinocytes as non-chemical tactic cells that do not move by self-propulsion. The force balance for the *i*th cell is2.3where vector *x _{i}* is the cell-centre coordinate, and

*μ*is the viscosity coefficient. The first term in equation (2.3) is the viscosity of the epidermis. and are repulsive and adhesive forces between neighbouring cells, which are sums of forces derived from all individual pairwise contacts,2.4where and

*Ω*(

*i*) denote sets of cells overlapping and neighbouring with the

*i*th cell, and and are force vectors produced onto cell

*i*by interaction (repulsion or adhesion) between cells

*i*and

*j*. By symmetry, and . The model treats individual cells as rigid-body agents and uses the extent of virtual cell overlap to determine the repulsive force (figure 3

*a*). Adhesion between two cells is a function of their spatial distance [39,40] (figure 3

*b*). Computation of and is given in the electronic supplementary material.

Previous agent-based models [9,11,41] treated keratinocytes as identically sized circles or spheres. However, cells progressively adopt varied shapes and sizes at different stages of differentiation. Cells in an outer layer generally have a more flattened cell body and larger surface area, compared with cells in layers underneath. Tissue location in the body can also influence cell geometry. For example, compared with the mean basal cell diameter of 6–8 µm at the forearm and hand [42], at unexposed sites in adult tissue [3], the average cell diameter in the proliferative compartment is about 10 µm, whereas the average differentiated cell diameter is about 16 µm. For simplicity, we use ellipsoids to model geometric heterogeneity in cell morphology and size. Cell types are distinguished by the mean major-to-minor axis ratio and the mean nominal size.

The basement membrane of the epidermis is undulant with rete ridges extending downward between the dermal papillae (figure 3*c*). In the adult human epidermis, the average rete ridge height is about 40 µm, and about six rete ridges along 1 mm cross-sectional tissue length were observed [43], which changes with age. The basement membrane is modelled by periodically repeating Gaussian functions (see the electronic supplementary material). Parameters of cell sizes and epidermis thickness are listed in table 2.

## 3. Results

The model recapitulates two important measures of epidermis homoeostasis: cell counts in different layers and epidermal turnover times in different compartments.

### 3.1. Homoeostatic cell density distribution

Mean homoeostatic cell densities can be analytically calculated from the ordinary differential equations in table 1 as follows (see table 3 for definitions of parameters):3.13.23.33.43.53.6and the total cell density is given as3.7The density of each cell category is proportional to the epidermis capacity of SCs, and the ratio between densities of any pair of cell types is a constant. Therefore, the observed cell density distribution can be used to identify kinetic parameters. The homoeostatic TA-cell density is3.8Kinetic parameters must satisfy two necessary conditions to establish a physiologically proper steady state3.9and3.10Given the small backconversion rate (*k*_{−1} and *k*_{−2} assumed at 10^{−6} per day, 3–5 orders of magnitude smaller than the SC and TA cell proliferating rate constants, table 3), the above conditions simplify to3.11The first condition ensures the establishment of a viable SC population, whereas the second condition prevents an unchecked growth of the epidermis because the proliferation of TA cells is not limited by a maximum capacity in the model. These conditions extend a previous treatment that required a precise balance between self-proliferation and symmetric division in progenitor cells when apoptosis was neglected [17,25]. Apoptosis is a rare event compared with progenitor self-proliferation and symmetric division (e.g. and ; see table 3) and plays a secondary role in regulating epidermal homoeostasis.

A recent study of progenitors in mice interfollicular epidermis by Mascré *et al.* [19] suggested that in homoeostasis SCs proliferate by four to six division events per year, 10–20 times slower than proliferation and differentiation (about 1.2 events per week) of committed progenitors (TA cells in our model). This information is reflected such that *k*_{1a,h} + 2*k*_{1s,h} is an order of magnitude smaller than *γ*_{2} + 2*k*_{2s}. However, considering that TA cells have a much larger population than SCs (*p*_{ta}/*p*_{sc} ≈ 5 [32]), rates of SC asymmetric and symmetric divisions into TA cells (*k*_{1a,h} + 2*k*_{1s,h}) must be about fivefold greater than the difference between TA-cell self-proliferation and symmetric division, i.e. *k*_{2s} + *β*_{2} − *γ*_{2} based on equation (3.2), assuming backconversion rates (*k*_{−1} and *k*_{−2}) are negligible. Therefore, the model predicts that mild perturbations on TA-cell dynamics including proliferation and apoptosis may have substantial effects on the homoeostatic keratinocyte population.

The model recapitulates the dynamics and the homoeostatic cell density distribution and the epidermal turnover time, in both deterministic and stochastic simulations (figure 4*a*). Starting from an initial SC population, the cell density rapidly increases to reach about 80% of the homoeostatic density in the transient phase of the first 100 days and significantly slows down when it approaches homoeostasis. The simulation reached homoeostasis in about 2000 days to a total density at 101 684 cells mm^{−2}. The transient dynamics is more than an order of magnitude slower than that in the model by Grabe & Neuber [9,61], where the authors reported a near 2000 h transient dynamics from an initial SC population to homoeostasis. This discrepancy is caused by a slow SC dynamics of four to six cell events per year reported recently [19], which also agrees with an early-reported cell cycle of 100–200 h in a 20% growth fraction in an animal model [62]. Therefore, the SC dynamics (equation (2.1)) becomes a rate-limiting step to the dynamic establishment of homoeostasis, notably when the system approaches homoeostasis.

Figure 4*b* shows a histogram of the simulated steady-state cell density distribution over different cell types. SCs, TA cells and GA cells in the stratum basale consist of 22% total cell population, in which proliferating cells account for 13.2% of the total cell population, consistent with the experimental data [32,48,49]. Studies by Bergstresser *et al*. [50] also suggested that 30% of nucleated keratinocytes were in the basal layer where SCs occupied 10% of the population [32]. SP and GC are the majority, consisting of 51.3% of the total keratinocytes with the population of SP cells nearly two times that of GC, in agreement with a previous prediction by Grabe & Neuber [9,61]. CC consist of 26.6% of the total cell population. Bauer *et al.* [48] measured a mean nucleated keratinocytes density of 75 346 cells mm^{−2}, whereas that of the non-nucleated CC was estimated to be about 18 000 cells mm^{−2}, consisting of 19.3% total cell population [49].

Figure 4*c* shows snapshots of temporal evolution of two-dimensional epidermal stratification from an initial group of SCs distributed along the basement membrane to homoeostasis. We compute the cell population dynamics and the cell migration within an area of 1 mm in length by 10 µm in width. The cell density is defined over a surface area number of keratinocytes per mm^{2} without explicitly considering the epidermis height. The choice of 10 µm (approximately mean cell size) is to visualize a two-dimensional single layer of keratinocytes. The simulated tissue histology shows that the thickness of the nucleated epidermis is about 60 µm, aligning with observations ranging from 38 to 77 µm with a mean of 60 µm in adult tissue, with little variation across age groups [42,43,48]. A movie of the epidermis renewal process of normal tissue is available at http://www.picb.ac.cn/stab/epidermal.html.

### 3.2. The epidermal turnover time

Another common measure of epidermis homoeostasis is the epidermal turnover time *τ*. At the tissue level, *τ* is interpreted as a time required for replacing the entire epidermis with new keratinocytes. The epidermal turnover time varies significantly with age groups, tissue locations and cell densities. In earlier studies, *τ* has been reported approximately six to seven weeks in the volar forearm with a nucleated cell density of 44 000 mm^{−2} [63,64]. Based on a more recent count of nucleated cell density of 75 346 mm^{−2} on breast skin [48], Hoath & Leahy [49] suggested that *τ* should be calculated as 59.3 days. Renewal of a specific layer of keratinocytes takes a shorter turnover time. The turnover times of the stratum basale and the differentiated compartment were reported as being approximately 22 and 12 days [64], whereas the stratum corneum has a turnover time that varies from 14 days [51,64] to about 20 days in young adult [65]. Following the convention [64], we calculate *τ* as the ratio of total cell density to the rate of cell loss including desquamation and apoptosis, or alternatively, the ratio of total cell density to the rate of cell birth by SC and TA cell divisions because cell birth and death rates are balanced at homoeostasis3.12where the rates of cell growth and cell loss are3.13and3.14An analytical solution to *τ* can be obtained by substituting equations (3.1)–(3.6) into the above equation. Apoptosis events represent rare ramifications from the central transition pathway in normal epidermis renewal. Rate constant of an apoptotic process *β*_{i} (10^{−4}–10^{−5} d^{−1}) is much smaller than the differentiation and desquamation rate constants *k _{i}* and

*α*(about 10

^{−1}d

^{−1}, see table 3). We can approximate

*r*

_{loss}by the rate of desquamation. By also neglecting the backconversions, we have3.15where3.163.173.18where

*τ*is the sum of contributions by sub-compartment turnover times as keratinocytes migrate from the proliferative compartment (

*τ*

_{prolif}) through the differentiated compartment (

*τ*

_{diff}) and then through the stratum corneum (

*τ*

_{corn}). We note that

*τ*is independent of SC self-proliferation

*γ*

_{1}and the SC capacity . These two parameters determine the steady-state SC density

*p*

_{sc}. Both cell density and cell growth rate are proportional to

*p*

_{sc}, which masks the effect of SC self-proliferation dynamics on the turnover times. The turnover time of a single cell category is the inverse of the rate constant for the transit to the next cell category. Therefore,

*τ*can be dominated by a rate-limiting step in the cell growth and differentiation cascade. From the calculation by equation (3.12) using parameters in table 3,

*τ*is about 52.5 days, partitioned into 7, 31.5 and 14 days in the proliferative compartment, differentiated compartment and the stratum corneum, respectively, within the wide range of reported adult-tissue measurements between 39 and 75 days [3,51,63].

### 3.3. Pathogenesis of psoriasis

To investigate an important epidermis disorder, we extend the above model to study the onset and recurrence of psoriasis and its management. Psoriasis, an immune system-mediated chronic skin condition, is characterized by an overproduction of keratinocytes accompanied by inflammation [66], resembling features found in many autoimmune diseases. Here, we focus on studying the most prominent type of the disorder, psoriasis vulgaris, a plaque-formed scaly silvery patch, affecting the majority of psoriasis patients.

Current studies are inconclusive about whether psoriasis (i) arises from genetically hyperproliferative progenitor keratinocytes, or (ii) is alternatively induced by a faulty immune system, particularly the over-response of dendritic cells and T lymphocytes to unresolved self-antigens, which in turn produces excessive cytokines (such as TNF*α*) to promote keratinocyte proliferation, or (iii) more likely an intricate interplay between keratinocytes and the immune system [66,67].

To examine the pathogenesis of psoriasis, we emphasize an interplay between keratinocytes overproduction and responses by the immune system. Our fundamental hypothesis is that a psoriatic epidermis assumes two groups of keratinocytes: normal and psoriatic, maintaining a parallel homoeostasis. The psoriatic keratinocytes are derived from a hyperproliferative SC population co-residing with the normal SCs at the basement membrane, which compete for limited SC niches (see figure 5 for illustration). This hypothesis can be justified by the existence of intrinsically hyperproliferative SCs or SCs that are more responsive to growth stimulants originated from an activated immune system. We modify the above model to describe the competition for niches between the normal and psoriatic SCs. Dynamics of psoriatic TA and non-proliferative cells are governed by rate equations as in table 1 for the normal keratinocytes, however, with different parameter values to generate known phenotypes in psoriatic plaques, with an exception of GC, which are missing in psoriasis. The dynamics of psoriatic SCs in our model is similar to the model of a spruce budworm outbreak by Ludwig *et al*. [68], which models the population growth of a single species under predation.3.193.203.21The psoriatic tissue activates the immune system to combat disease cells. We assume that repertoires for normal *p*_{sc} and psoriatic SCs are both limited by the available niche environment and psoriatic SCs can acquire a larger growth capacity. Parameter *λ* (more than 1) accounts for the fold increase in the growth capacity accessible to SCs. The density of normal SCs remains limited by , which is invaded by a fraction (1/*λ*) of psoriatic SCs. Equation (3.21) models the immune activities triggered by psoriatic SCs. An activated immune system induces apoptosis of psoriatic SCs. The activity of the immune system (the killing rate, ) is directly regulated by the psoriatic SC density, under the assumption that the immune system is activated in a faster timescale than the tissue growth. The immune response is significantly activated when exceeds a threshold parametrized by *K _{a}* and is saturated at the maximum rate

*K*at when the psoriatic SC population overwhelms that of cytotoxic T cells. This approach hypothesizes that the immune system combats disease SCs, but does not exclude the commonly believed role played by the immune system of inducing keratinocyte overproduction, even though the model does not explicitly couple the immune system to SC proliferation.

_{p}Proliferation rate constants *γ*_{1}, *k*_{1s} and *k*_{1a} for the normal SCs are regulated by the total TA cells, , similar to equation (2.1)3.22In comparison, we assume that psoriatic SCs are not subject to regulation by the TA cell population and proliferate with rates *ρ*_{sc}-fold higher than the homoeostatic rate constants (*γ*_{1,h}, *k*_{1a,h} and *k*_{1s,h}) of normal SCs. We assume that the immune response substantially switches on when the psoriatic SC population reaches 10% of the SC population in the normal tissue. This assumption is used to parametrize the steepest change in the removal rate in response to , which sets the half-activation density at . The model does not consider immune responses against cells derived from psoriatic SCs by observing that reduction in the SC population results in a subsequent reduction in the derived keratinocyte population.

Despite its mechanistic uncertainties, psoriasis plaques have well-defined tissue-level phenotypes, making it a good candidate for study by a predictive model. Depending on its severity, a plaque exhibits a two to five times increase in the total cell density [52,53,69] with a relatively higher growth in the proliferative compartment compared with the non-proliferative compartment [55,56]. A disordered tissue usually loses the granular layer due to abnormal differentiation and contains a subset of nucleated CC. A psoriatic plaque also has a turnover time several fold faster [51,70]. More specifically, studies of the cell kinetics of psoriasis found (i) a significant increase of cell-cycle marker Ki-67 in psoriatic tissue without much change in cell-cycle time [35,71], suggesting a substantial increase in growth fraction; (ii) the transit time of keratinocytes through the differentiated compartment is shortened to 48 h from 240 to 330 h, five to seven times faster than in normal tissue [54]; (iii) the transit time through the corneum is also shortened from 14 to 2 days [51]. These factors together result in a decrease in the epidermis turnover time (equation (3.15)). In addition, the cell apoptotic index decreases nearly fourfold from 0.12 to 0.035% [36], making a further contribution to keratinocyte overproduction. Table 3 lists parameter values based on the above observations, where coefficients *ρ*_{ta}, *ρ*_{tr} and *ρ*_{de} are fold changes over rate constants in normal kinetics of TA-cell proliferation (*γ*_{2}, *k*_{2a} and *k*_{2s}), transit in the non-proliferative compartment (*k*_{3} and *k*_{4}) and desquamation (*α*), respectively. Similar to variations in normal tissues, the severity and phenotype of psoriasis vary widely across individuals and disease subtypes and therefore for any specific condition or study the model should be parametrized accordingly.

As the main result, the model predicts two interconverting homoeostatic modes of psoriatic tissue (see appendix A for analytical details): (i) a *disease* state, which generates psoriasis phenotypes of keratinocytes overproduction and a shortened epidermal turnover time when psoriatic SCs outcompete normal SCs for available niches and overwhelm the immune system; and (ii) a *quiescent* state, which predicts a coexistence of a small number of psoriatic cells with a dominant population of normal keratinocytes when the immune system keeps the psoriatic SC population low. The quiescent mode with remanent psoriatic SCs is a symptomless state for a psoriasis-susceptible tissue and may relapse into the disease state given favourable conditions. On the other hand, the disease state may be reverted back to the quiescent state by a properly designed treatment.

The model is parametrized (table 3) to produce a specific psoriasis phenotype, in which the total cell density reaches 217 000 mm^{−2} at homoeostasis, more than two times the normal epidermis (101 600 mm^{−2}), with the pathological keratinocytes consisting of 99.5% of the total population. The relative growth of the proliferative compartment over the non-proliferative compartment is about three versus two times the normal densities (table 4). The epidermal turnover time is shortened more than five times from 52.5 to 9.8 days.

The psoriatic epidermis may retreat to the ‘quiescent’ state and achieve a remission provided that the psoriatic SC density can be managed below a threshold value (see appendix A). Narrow-band 311 nm controlled UVB irradiation is known as an effective treatment for managing psoriasis [72]. For example, a recent study by Weatherhead *et al*. [53] applied sequential episodes of 0.75–3 minimal erythemal dose UVB irradiation to achieve plaque remission by inducing strong apoptosis of proliferative cells. Figure 6*a* shows a model simulation of psoriasis remission after a simulated sequence of seven-episode UVB irradiations. Each UVB irradiation episode is simulated by increasing apoptosis rate constants 60 000-fold indiscriminately for normal and psoriatic SCs and TA cells for 48 h followed by a 8-h resting interval before starting the next episode. The entire treatment lasts 16 days. Model simulations unveil an intriguing interplay between the dynamics of psoriatic and normal cells. Upon the initiation of irradiation, the total cell population first declines due to the UVB-induced apoptosis in SCs and TA cells. Each episode of UVB irradiation induced apoptosis in about 22% of the SCs in the pre-episode population (see figure 6*a*, inset). The population of TA cells declined more substantially due to the combined effects of increased apoptosis and reduced SC symmetric and asymmetric divisions. A mild rebound of cell densities happens in each resting interval because of a continuing hyperproliferation of psoriatic SCs and TA cells. At the end of the treatment, the total keratinocyte density dramatically drops more than 95% from 217 000 to 12 100 mm^{−2}, with normal and psoriatic SCs, respectively, reduced to 260 and 960 mm^{−2}. Post-treatment the SC population continues to decline to the ‘quiescent’ steady state due to a relatively stronger immune response (figure 6*b*; electronic supplementary material, figure S3*c* and related text), which later brings down the psoriatic cell density to a minimum (less than 0.5% of the total cell density, table 4). The GC become visible during the recovery. The last phase indicates a recovery of keratinocytes derived from the normal SCs that reclaim their niche repertoire by a slower kinetics (figure 6*c*). The downstream differentiated cells follow similar dynamics of remission. Simulated dynamics are similar if the model considers UVB-induced apoptosis in all nucleated cells (results not shown).

A treatment with less UVB irradiation episodes and/or inadequate intensity may fail to clear a psoriatic plaque, which eventually returns to the disease state once the treatment stops due to an insufficient loss in psoriatic SCs. Figure 6*b*,*c* shows that after a six-episode UVB irradiation treatment, the psoriatic SC and TA cell populations bounce back to the disease state after terminating the treatment. The total cell density drops to 18 600 mm^{−2} at the end of the sixth episode with the normal and psoriatic SC densities as 247 and 1263 mm^{−2}, respectively. Interestingly, the end-treatment normal SC count is slightly less than that after a seven-episode treatment, implying that the increased normal SC proliferation rate due to loss of TA cells well below the healthy level, offsets the cell loss caused by apoptosis. Shortly after the treatment has stopped, both psoriatic and normal SCs and TA cells started to increase slowly, but later the density of psoriatic cells (figure 6*b*) rapidly expands and outcompetes the normal cells (figure 6*c*) that retreat from a maximum to the steady state at a lower level. During the entire course of the treatment, the cytotoxic rate remains below the growth rate of psoriatic SCs, giving no chance for the immune system to effectively reduce the psoriatic SC population (see electronic supplementary material, figure S3*d*).

Histograms in figure 6*d*,*e* show that a well-designed treatment can manage the psoriatic tissue to the quiescent state that is almost phenotypically indistinguishable from the healthy tissue regarding cell density distribution and the turnover time. Histologically, psoriasis causes thickening in the stratum corneum and the differentiated layer as well as an expanded proliferating compartment with more protruding rete ridges [73]. A dynamic model of rete ridge remodelling, simulation snapshots of the homoeostatic psoriatic epidermis, tissue under treatment and recovered tissue can be found in the electronic supplementary material, figures S3 and S4). Movies of simulations (normal and UVB-treated psoriatic tissues) are provided at http://www.picb.ac.cn/stab/epidermal.html.

## 4. Discussion

We presented a hybrid model that simulates and visualizes spatio-temporal dynamics of epidermal homoeostasis. The model represents an efficient approach that separates the computation of cell kinetics from that of an agent-based cell migration. Compared to previous agent-based models [8,9,11,12], our population kinetics model describes cell proliferation, differentiation and cell death as empirical rate processes and can be simulated by integrating the governing ordinary differential equations (table 1) or the master equations by a kinetic Monte Carlo algorithm. The cell population kinetics model can be combined with a two-dimensional cell migration model to visualize dynamics of epidermis renewal and stratification. The model reproduces observed characteristics of the normal epidermis. Model analysis and simulations show that balancing cell production and cell loss in each sub-compartment is critical to establishing and maintaining proper epidermis homoeostasis (figure 4).

The current model has some addressible limitations: (i) the cell population kinetics does not explicitly incorporate specific intracellular and extracellular factors that affect the dynamics and steady state of the epidermis homoeostasis. However, physiological and physical factors including age and UV irradiation as well as many commonly investigated signalling molecules can be coarsely coupled to model parameters such as proliferation and differentiation rate constants and morphology of keratinocytes and the epidermis. (ii) We neglected the effects of backconversions from TA cells to SCs and from differentiated cells to TA cells by assuming their minimal impact. These processes could be worth a close examination as suggested in a recent theoretical study [74] that showed rare backward transitions may cause catastrophic outcomes such as cancerous growth. (iii) Technically, as intensive modelling and computation is made possible by high-performance hardware [75,76], our two-dimensional cross-sectional model can be extended to simulate a more realistic three-dimensional epidermis even though we expect that the qualitative results obtained from the two-dimensional model remain valid in a three-dimensional model.

As an important application, a non-trivial extension to the above model allows us to investigate the pathogenesis of psoriasis, an immune-mediated skin disorder. Genetic origins of psoriasis have been recently explored by an increasing number of genome-wide association studies that identified a multitude of psoriasis susceptibility loci [77–80]. Many psoriasis-associated loci are connected to genes in the immune system (e.g. MHC class I molecules) and proteins expressed in keratinocytes, suggesting a complex nature of the disease [81]. However, the mechanistic epidermis–immune system interactions implicated by these loci are yet to be resolved.

In this study, we propose an alternative hypothesis of interactions between the immune system and keratinocytes, in which the disordered epidermis maintains a parallel homoeostasis of both normal and psoriatic keratinocytes, derived from respective SC populations. We examine this hypothesis in an extended model and demonstrate that treatment by UVB irradiation with consecutive episodes can potentially manage the disease and achieve a remission of the psoriatic phenotype.

Psoriasis has recently been studied by agent-based models [53,61]. The model by Grabe & Neuber [61] was able to generate the psoriasis phenotype of an increased cell density and a shortened epidermal turnover time, by adjusting the fractional time of TA cell proliferation. This parameter characterizes the amount of time for proliferation during a constant life time of a TA cell, which in our model is embedded in the TA cell self-proliferation rate constant *γ*_{2}. Increasing *γ*_{2} and keeping symmetric division rate constant *k*_{2s} unchanged (equivalent to keeping a constant TA cell life time) does increase TA cell population (equation (3.2)). To obtain a relative growth of the proliferating compartment, rate constants for cell differentiation must have relatively higher increases than *γ*_{2}, which was achieved by modulating the Ca^{2+} gradient in the Grabe and Neuber model. The model however did not propose possible management that can target the hypothesized mechanism. For example, it is not obvious how apoptosis induced by episodes of UVB irradiation can attain a remission via recovering the normal TA cell proliferating time. Weatherhead *et al*. [53] developed a model to simulate UVB-induced apoptosis in SCs and TA cells, which was able to demonstrate psoriasis remission after a few episodes of UVB irradiation. The model assumed a constant pool of SCs that derived TA cells by asymmetric division and the UVB-induced apoptotic hyperproliferative SCs were replaced with normal SCs by symmetric divisions. This model assumption consequently led to a permanent reduction in cell density after each UVB irradiation treatment and therefore did not explain relapses of psoriatic phenotypes once an ineffective treatment ends or recurrence of the disorder.

Interactions between the immune system and keratinocytes considered in our model can serve as a conceptual basis for interpreting the pathogenesis of psoriasis. Especially, we showed that the act of the immune system cytotoxicity against psoriatic keratinocytes plays a pivotal role in the onset, remission and recurrence of the disease phenotype. The most common paradigm considers that faulty immune responses triggered by unknown self-antigens or pathogens (introduced by injuries or trauma, known as Koebner phenomenon) produce cytokines and growth factors that promote keratinocyte hyperproliferation, immature differentiation and skin inflammation, establishing the psoriasis phenotype [14,70]. Drugs that inhibit T cell activity and cytokine productions do improve psoriatic conditions [82]. By contrast, our model demonstrates that a psoriasis lesion develops when the immune system is genuinely weakened or locally overwhelmed by a large population of psoriatic keratinocytes, staging up a chronic condition. This prediction may explain high occurrence and increasing severity of psoriasis in HIV-infected, especially late-stage AIDS patients with substantially compromised immune systems when CD4^{+} and naive CD8^{+} T cell counts substantially decrease [83,84]. Furthermore, the onset age of psoriasis has been known to have two separate populations, type I (early onset in patients younger than 40 with a peak at age 20) and type II (late onset after age 40 with a peak at about 60) [85]. Our model speculates that a vigorous immune system at a younger age can stimulate a strong hyperproliferation (a large *ρ*_{sc}) in a psoriatic epidermis and causes manifestation of plaque phenotypes. One the other hand, a weak immune system (a small *k _{p}*) at an older age can also have an equivalent effect.

The model prediction provides an alternative (a less explained) perspective of the pathogenesis of psoriasis, suggesting that psoriasis is a parallel epidermal homoeostasis due to heterogeneity in SC clones and that the immune system as a double-edged sword plays two essential however opposing roles: (i) cytotoxic (CD8^{+}) T-cells recruited to the epidermis induce apoptosis in psoriatic SCs, which is implicated by studies that demonstrated CD8^{+} T cell (especially, CD45RO+ memory subtype) count and cytotoxic proteins including perforin and granyme B substantially increase in psoriatic lesions [83]; and (ii) immune activities in the meantime promote progenitor keratinocytes proliferation by producing a multitude of cytokines and growth factors (TNF, IFN-*α*, IFN-*γ*, IL-17, IL-22, IL-23, etc.). Balance between the two acts determines the outcome of the disease. Our model predicts that the psoriasis-susceptible tissue is a bimodal system and can switch between a non-symptom state and a phenotypical state, potentially explaining the recurring nature of the disorder and suggesting the feasibility of disease management by an effective treatment. This perspective and the model could be in general applicable to other autoimmune diseases in regenerative tissues.

We showed in UVB phototherapy simulation that psoriasis plaque clearance can be attained by inducing strong apoptosis in keratinocytes. The model only considered UVB-induced apoptosis in keratinocytes. UVB irradiation may promote proliferation and differentiation in skin cells [86,87] or may alter immune responses. These effects could be incorporated and be examined in the model. For example, as suggested in equation (3.20), rebalance in self-proliferation *γ*_{1,h} and symmetric division *k*_{1s,h} and change in immune activities may affect the population of psoriatic SCs and thus the therapeutic outcome. Experiments did not demonstrate whether the UVB irradiation causes non-apoptotic cell death that was not reflected by an apoptosis marker [53]. Regardless of the actual mechanism, induced cell deaths will result in a shift of the SC density from the disease state to the quiescent state across the phase boundary, an important parameter that determines the design of a phototherapy regimen, including irradiation dosage in each treatment episode, the number of episodes and the time interval between consecutive episodes.

The model also suggests the limitation of UVB irradiation or similar treatments that attempt to achieve plaque remission by killing keratinocytes below a critical threshold. First, psoriasis plaques may recur under a temporally weakened immune system or a transient burst of cell proliferation caused by conditions such as wound healing. Second, the psoriatic severity may be worsened under a genuinely weak immune system (with a low *K _{p}* and/or a high threshold

*K*), in which the disease phenotype persists and cannot be adequately reverted by simple reduction in psoriatic SCs because the system does not possess a quiescent state (region I in figure 7

_{a}*b*). In the latter case, other treatment options such as cytokine-targeting biologic and small-molecule drugs should be considered to at least shift the disorder to the bimodal region where phototherapy becomes effective. However, we note that phototherapy may be effective through alternatively mechanisms other than induction of apoptosis [88,89], by which it alone may switch a plaque from persistent disorder (region I) to bimodal (region II) or symptomless (region III) as in figure 7

*b*.

Further insight from the model is that the bimodal psoriasis exhibits hysteresis, by which a mode of the epidermis, psoriatic or quiescent, once reached under a favourable condition, will tend to be relatively stable. For example, a symptomless epidermis as illustrated in figure 7*b* may switch to the disease state when *K*_{p} is reduced beyond its lower threshold because of weakening in the immune system. However, restrengthening the immune system to revert the disease state back to the quiescent mode requires elevating *K*_{p} in the model beyond the upper threshold, implying that the disease is resistant to mild natural or induced perturbations.

One definitive experimental test of our model hypothesis should aim at simultaneously tracing clones of psoriatic and normal SC lineages. Because of recent advancement in techniques, studies of epidermal SCs become central to the understanding of epidermal homoeostasis. Especially, emerging powerful *in vivo* lineage tracing technology allows dynamic monitoring of SCs and progenitors clone formation and differentiation and can achieve multicolour simultaneous tracing of separate clones [27]. The technique has already been pioneered to study progenitor fates in normal skin tissue [17] and benign epidermal tumours [90] and can be potentially applicable to identify the coexistence of normal and psoriatic progenitors in psoriasis. One potential challenge is that in either the quiescent or disease state, one type of SCs is present as a small population (about a 1 : 20 ratio, according to table 4), making it difficult for differential labelling. Alternatively, an experiment may tend to identify and isolate the two competing effects by the immune system on the epidermis: (i) apoptosis in epidermal progenitors by T-cell cytotoxic activities and (ii) hyperproliferation due to elevated cytokines and growth factors.

## Funding statement

This work was supported by LVMH Recherche through Sprim Inc. and institutional fund to J.Y.

## Acknowledgements

We thank Weiren Cui, Eric Perrier, Michaël Shleifer, Gallic Beauchef and Delphine de Queral for helpful discussions.

## Appendix A. Analytical details of the psoriasis model

**A.1. Psoriasis as a bimodal switch**

To illustrate the principle and simplify the analysis, we only present a two-dimensional model that describes dynamic interactions between normal and psoriatic SCs. The complete model that produced the results (figure 6) is more complex and of a higher dimension because of the feedback regulation by the TA cell population on the proliferation of normal SC (equation (3.22)), requiring analysis of the dynamics of normal and psoriatic TA cells together. Again, for convenience in the analysis, we neglect the influences of backconversions and apoptosis. This simplification allows us to show insights of the behaviours of the model without much mathematical complication. We work with normalized quantitiesA 1Ratios between the normal SC proliferation rate constants remain constant despite their dynamic modulations by the TA cell population (equation (2.1)), i.e. *k*_{1s}/*γ*_{1} = *k*_{1s,h}/*γ*_{1,h}. Ignoring the trivial steady states of zero density (*x* = 0 and/or *y* = 0), based on equations (3.19)–(3.21) we have the steady-state equations for SC densities asA 2andA 3whereA 4andA 5Equation (A 2) can solve for one or three steady-state psoriatic SC densities, only depending on choices of two parameters, *w* and *u*, which is illustrated in figure 7*a*.

To study plaque remission and relapse, we are interested in the region where the system has two stable steady states, ‘quiescent’ and ‘disease’, and one unstable ‘transition’ steady state. This scenario is illustrated in figure 7*a* with three intersections between the line *w*(1 − *x*/*u*) and the non-parametric curve *x*/(1 + *x*^{2}). Determination of the stability of these fixed points requires analysis similar to the spruce budworm model by Ludwig *et al*. [68]. Figure 7*a* shows that the ‘quiescent’ state is more sensitive to *w*, whereas the ‘disease’ state is determined by both *u* and *w*. This observation provides a guideline to parametrize the model to attain appropriate psoriatic phenotypes by identifying disease-related parameters *ρ*_{sc}, *λ* and *k _{p}*. The ratio

*ρ*

_{sc}/

*k*emerges as the key parameter that characterizes the interplay between proliferation of keratinocytes and activity of the immune system (magnitude

_{p}*K*and threshold

_{p}*K*), playing a critical role in the interpretation of the pathogenesis of psoriasis.

_{a}Bifurcation diagrams in figure 7*b* of steady-state SC densities shows three classes of behaviours as *ρ*_{sc}/*k _{p}* varies. A small population of psoriatic SCs survive in the ‘quiescent’ state (

*ρ*

_{sc}/

*k*< 159 day) under a vigorous cytotoxic activity (a large

_{p}*k*, i.e. a low threshold

_{p}*K*and/or a high magnitude

_{a}*K*) and/or a moderate SC proliferation (a small

_{p}*ρ*

_{sc}). By contrast, a weak immune system (a small

*k*, i.e. a high threshold

_{p}*K*and/or a low magnitude

_{a}*K*) cannot adequately counterbalance the psoriatic hyperproliferation (a large

_{p}*ρ*

_{sc}) and thus the epidermis assumes a persistent ‘disease’ state (

*ρ*

_{sc}/

*k*> 440 day). A bimodal system with an intermediate

_{p}*ρ*

_{sc}/

*k*has the potential to switch between the quiescent and disease modes when conditions change. Figure 7

_{p}*c*shows a two-dimensional phase plane of normal and psoriatic SC densities. Starting from initial SC densities, located inside region I or II, a temporal trajectory will be attracted to the quiescent or the disease state, respectively. Those starting on the boundary of the two regions will in theory converge to the transition state, which is unsustainable because slight perturbations will dislocate a trajectory into either region I or II.

**A.2. Cell density and turnover time**

Below, we provide equations for homoeostatic cell densities and the turnover time in the psoriatic model. As shorthands, we defineA 6The ratio *ξ* = *g*/*g*_{h} is a function of the steady-state TA cell population () by equation (3.22). The total steady-state cell density is given asA 7whereA 8andA 9Note the lack of the GC layer in the psoriatic cell population. The proliferative and differentiated cell populations areA 10andA 11CC are not included in the differentiated compartment. By comparison, the cell densities in the healthy tissue areA 12andA 13We can calculate the growth of the proliferative compartment relative to that of the non-proliferative compartment asA relative increase in the proliferative population requires the above ratio to be greater than 1. In psoriasis, *ξ* < 1, reflecting a decreased proliferation rate in the cohabitating normal SC population due to a much elevated TA cell population. It is unclear whether SCs and TA cells have differential increases in proliferation rates. For the lack of information, we assume *ρ*_{sc}/*ρ*_{ta} ≈ 1. Therefore, *ρ*_{tr} > *ρ*_{sc} becomes a sufficient (but not a necessary) condition to produce a relative increase in the proliferative compartment, requiring that the fold increase of transit rate in the differentiated compartment is higher than the fold increase of division rates in the proliferative compartment. This result agrees with an earlier prediction from Heenen *et al.* [55], showing that a relative growth in the proliferative compartment requires keratinocyte hyperproliferation and in the meantime an increased transit time in the differentiated compartment. A proportional increase in both TA-cell proliferation rate and transit rate (*ρ*_{sc} = *ρ*_{tr}) does not cause relatively differential growth of the two compartments unless difference exists between increases of proliferating rates of SCs and TA cells, *ρ*_{sc} ≠ *ρ*_{ta}. Fast migration of differentiated keratinocytes also offsets overgrowth in the proliferative compartment. An important insight from the model is that increasing keratinocyte population in psoriasis is not merely caused by elevated proliferation rates but also requires the increase in growth capacity in the progenitor repertoire. Observed morphological remodelling of the psoriatic epidermis suggests increased growth niches for SCs [73].

Similarly to equation (3.15), the epidermal turnover time can be approximated as the ratio of the total cell density to the rate of desquamationA 15which can be analysed for the two psoriasis modes (quiescent or disease). At the quiescent mode, *ξ* ≈ 1 and , the turnover time is approximated byA 16converging back to *τ*_{tot}. At the disease mode, *ξ* < 1 and . The turnover time is close toA 17

Both cases do not account for the relatively much smaller cell loss by the cytotoxic effect from the immune system. Compared to the turnover time of the healthy tissue *τ*_{tot} (equation (3.15)), of the psoriatic tissue is clearly shortened at each stage from the proliferative compartment (the first term) to the differentiated compartment (the second term) and to the corneum (the third term). The turnover time of the granular layer is ignorable because of the minimal normal cell population.

- Received September 25, 2014.
- Accepted December 10, 2014.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.