Dynamics of a network mediated by IL-36 and involved in the pathogenesis of psoriasis

The pathogenesis of the inflammatory, chronic, and common skin disease psoriasis involves immune cells, skin cells (keratinocytes), and the cytokines they secrete. Hyperproliferation and abnormal differentiation of keratinocytes are hallmarks of the disease. The roles of cytokines such as TNFα, IL-15, IL-17, and IL-23 in psoriasis have been studied through mathematical/computational models as well as experiments. However, the role of proinflammatory cytokine IL-36 in the onset and progression of psoriasis is still elusive. To explore the role of IL-36, we construct a network embodying indirect cell–cell interactions of a few immune and skin cells mediated by IL-36 based on existing knowledge. We also develop a mathematical model for the network and perform a global sensitivity analysis. Our results suggest that the model is most sensitive to a parameter that represents the level of cytokine IL-36. In addition, a steady-state analysis of the model suggests that an increase in the level of IL-36 could lead to the hyperproliferation of keratinocytes and, thus, psoriasis. Our analysis also highlights that the plaque formation and progression of psoriasis could occur through either a gradual or a switch-like increase in the keratinocyte population. We propose that the switch-like increase would be due to a bistable behavior of the network toward either a psoriatic or healthy state and could be used as a novel treatment strategy.

1 Introduction

Psoriasis is one of the most common skin diseases, and approximately 2%–3% of the world’s population is affected by it (Parisi et al., 2013). It is inflammatory and chronic in nature, distinguished by the hyperproliferation of keratinocytes (skin cells) and infiltration of immune cells in the psoriatic lesions. The lesion could be scaly or in the form of pustules, characterizing different types of psoriasis. The involvement of immune cells, such as T-lymphocytes, dendritic cells, macrophages, mast cells, and neutrophils, in the pathophysiology of psoriasis has been well-established (Lowes et al., 2014). The interaction among these immune and skin cells is believed to be mediated by cytokines. Recently, the roles of cytokines TNFα, IL-15, and the IL-17/IL-23 axis have been demonstrated by constructing networks of indirect cell–cell interactions among immune cells and keratinocytes (Oza et al., 2017; Pandey et al., 2023). Based on the steady-state behavior of the mathematical models of these networks, authors have demonstrated how an increase in the level of modeled cytokines could lead to the progression of psoriasis. Cytokines involved in psoriasis vulgaris have been more extensively explored than those in other types of psoriasis (Chen et al., 2021; Sachen et al., 2022), such as generalized pustular psoriasis (GPP) in which inflammatory cytokine IL-36 is believed to have a key role (Johnston et al., 2017). However, the underlying mechanism of its role in the pathogenesis of psoriasis vulgaris is still elusive. To address that, we have constructed a network for the indirect cell–cell interactions among immune and keratinocyte cells mediated by the cytokine IL-36. Our network is built on the published observations about the interactions of epidermal keratinocytes, T cells, macrophages, and dendritic cells (DC). Modulation of cells through IL-36 occurs upon binding of IL-36 to the receptors present on the surface of modeled cells in the network. Keratinocytes (K) are a predominant source of IL-36 cytokines, and in an autocrine manner, they enhance the expression of IL-36 (Carrier et al., 2011). A skin injury in the presence of psoriasis stimulates the release of IL-36 by neighboring keratinocytes in a paracrine manner by releasing cathelicidin LL-37 from dead keratinocytes (Li et al., 2014; Zhang et al., 2016). Due to that, T cells and dendritic cells are recruited at the injury site (Samotij et al., 2021). The receptors of IL-36 are predominantly present on naive CD4⁺ T cells, and one of the isoforms of the cytokine, IL36β, is also constitutively expressed in these cells (Vigne et al., 2011). Many studies have reported that the CD4⁺ T cells present in the psoriatic lesion hyperproliferate in response to an IL-36 stimulus (Günther and Sundberg, 2014). It has been observed that bone marrow-derived dendritic cells show Th1 responses as a result of IL-36 induction (Vigne et al., 2011). The DCs activated by IL-36 secrete IL-12, and a synergistic effect of both cytokines results in the polarization of naive CD4⁺T cells toward Th1 (Vigne et al., 2011). In addition, IL-36 promotes Th17 and Th22 lymphocyte polarization through a similar DC-mediated mechanism (Mercurio et al., 2018; Samotij et al., 2021). It is also observed that IL-36-mediated T-cell activation by dendritic cells leads to the hyperproliferation of keratinocytes (Foster et al., 2014).

Macrophages are also one of the sources of IL-36 in a psoriatic lesion (Dietrich et al., 2016). Macrophages treated with IL36γ stimulated the release of TNFα and IL-23. These cytokines have multiple roles, such as stimulating the differentiation and proliferation of T cells (Bridgewood et al., 2018), inducing cytokine release by keratinocytes (as illustrated by (D’erme et al., 2015)), and promoting the maturation of dendritic cells (Sachen et al., 2022). Human dermal macrophages have a high concentration of IL-36R on their surface, which is converted into a proinflammatory form (M1 phenotype) upon the action of IL-36 (Dietrich et al., 2016). Cytokine-stimulated macrophages increase the adherence of monocytes to endothelial cells and cause the polarization of lymphocytes, in this case, T cells, due to the upregulation of the IL-17/IL-23 axis (Bridgewood et al., 2018).

These indirect cell–cell interactions mediated by IL-36 led to the construction of a network, as illustrated in Figure 1. Sources and effect of IL-36 cytokines on the modelled cells presented in the network are summarized in Table 1. Here, we have developed a mathematical model for the constructed network mediated by IL-36 to explore the steady-state behavior of the network. We have denoted the IL-36 agonists (α, β, and γ) through a single IL-36 cytokine as they bind to the same signaling receptor (IL-36R) and are expressed in all the cell types in the model. The three separate cytokine agonists lead to the pathogenesis of psoriasis through hyperproliferation and, in some cases, maturation of the modeled immune cells. We have denoted them by a single cytokine that binds to the cell surface receptors IL-36R on the modeled immune cells to keep the model simple and remain trackable. We conducted a global sensitivity analysis to explore the impact of changing model parameters and to assess the sensitivity of the model outcome to each parameter. Our results demonstrate that an increase in IL-36 could lead to the hyperproliferation of keratinocytes, which is considered a hallmark of psoriasis. We also observed that the route of the plaque formation and or progression of psoriasis could be via a gradual or switch-like increase in the population of keratinocytes. This switch-like behavior could be used as a new treatment strategy, as argued in a few recent studies (Pandey et al., 2021; Pandey et al., 2023).

Figure 1

Figure 1. A schematic diagram explaining IL-36-mediated interactions among immune cells and keratinocytes involved in the pathogenesis of psoriasis. Here, L, K, D, and M denote T-lymphocytes, keratinocytes, dendritic cells, and macrophages, respectively. Different colored arrows represent different types of interactions among the cells involved.

Table 1

Table 1. IL-36 cytokines with their source and effect on different cell types. Here, S represents the origin of the cytokines, while + denotes a net increase in the cell population. These changes are outcomes of the indirect cell–cell interactions via modeled cytokines, that is, IL-36 agonists. The agonists represent the three receptor agonists that are members of the IL-36 cytokine family, that is, IL-36α, IL-36β, and IL-36γ. They initiate the IL-36 receptor signaling and activate various proinflammatory mediators like NFκB and MAPK.

2 Materials and methods

Based on the available published knowledge about the role of cytokine IL-36 in the pathogenesis of psoriasis and GPP, a schematic diagram is prepared to depict the indirect cell–cell interactions among skin cells (keratinocytes) and immune cells (dendritic, T cells, and macrophages) (Figure 1). Interactions depicted in the diagram elucidate the involvement of IL-36 in the proliferation, maturation, and differentiation of modeled cells in a psoriatic lesion. In a recent study, Pandey et al. (2021) developed a mathematical model and demonstrated the involvement of cytokines TNFα, IL-23/IL-17, and IL-15 in the pathogenesis of psoriasis (Pandey et al., 2021). The mathematical modeling framework they proposed is ideal for the present study. Therefore, we follow that framework for modeling the network mediated by IL-36 as follows.

2.1 Mathematical modeling framework

In our model, we consider the lesional skin of a psoriatic patient with a resident population of keratinocytes and immune cells (T cells, dendritic cells, and macrophages). We assume that all the modeled cells can be infiltrated into this lesion. This region of skin will be treated with different levels of IL-36 to investigate the role of that cytokine in the plaque formation and progression of psoriasis.

The dynamics of the modeled cell population are governed by following ordinary differential equations (ODEs) if the framework proposed elsewhere (Pandey et al., 2021) is applied:

$\frac{dx}{dt}={\alpha }_x+\beta \left(x\right)-{\gamma }_{x}x(1)$

Here, x represents the population of a modeled cell type that is, x = {K, D, T, M} where the populations of keratinocytes, dendritic cells, T cells, and macrophages are denoted by K, D, T, and M, respectively. The parameter α_x denotes a fixed rate of increase in the population of modeled cells x) through their migration toward the lesion or their differentiation/maturation in the lesion (independent of the cytokine IL-36). β denotes the net rate of change in the cell population that accounts for the cumulative effects of the cytokine in cell differentiation, maturation, and proliferation. The rate of decrease in the cell population through other processes, including apoptosis not mediated by IL-6, is denoted by the parameter γ_x. The function β is considered to be a function that saturates at β_x and depends on the population of all modeled cells or a few based on the interaction network shown in Figure 1. This dependence is modeled through the following two forms of β.

(i) when β increases due to the modulatory effect of the population of all cell types or a few

$\beta ={\beta }_x\frac{x^n}{k_x^n+x^n}(2)$

(ii) when β decreases because of the modulatory effect of the population of all cell types or a few

$\beta ={\beta }_x\frac{k_x^n}{k_x^n+x^n}.(3)$

In Eqs 2, 3, β_x represents the saturated value of function β due to the modulatory effect of modeled cells. Here, the parameter k_x represents the population of a modeled cell type that is required to show the half-maximal change of it or another cell’s population, depending on the interaction of the network shown in Figure 1. The slope of the function β is determined by the Hill-coefficient like parameter n. The rationale behind using the Hill-like function here is its wide acceptance as an empirical saturating function.

2.2 Indirect cell–cell interaction network mediated by cytokine IL-36

All modeled indirect cell–cell interactions driven by IL-36 have been depicted in Figure 1 in the form of a schematic network. The following system of ordinary differential equations will be obtained if the above-mentioned modeling framework is applied to model this interaction network.

$\begin{array}{cc}\hfill \frac{\mathrm{d}D}{\mathrm{d}t}& ={\alpha }_D+{\beta }_{IL36}\left(\frac{D^n}{k_D^n+D^n}\frac{T^n}{k_T^n+T^n}\frac{M^n}{k_M^n+M^n}\frac{K^n}{k_K^n+K^n}\right)-{\gamma }_{D}D\hfill \\ \hfill \frac{\mathrm{d}K}{\mathrm{d}t}& ={\alpha }_K+{\beta }_{IL36}\left(\frac{D^n}{k_D^n+D^n}\frac{T^n}{k_T^n+T^n}\frac{M^n}{k_M^n+M^n}\frac{K^n}{k_K^n+K^n}\right)-{\gamma }_{K}K\hfill \\ \hfill \frac{\mathrm{d}T}{\mathrm{d}t}& ={\alpha }_T+{\beta }_{IL36}\left(\frac{D^n}{k_D^n+D^n}\frac{T^n}{k_T^n+T^n}\frac{M^n}{k_M^n+M^n}\frac{K^n}{k_K^n+K^n}\right)-{\gamma }_{T}T\hfill \\ \hfill \frac{\mathrm{d}M}{\mathrm{d}t}& ={\alpha }_M-{\gamma }_{M}M\hfill \end{array},(4)$

where D, K, T, and M denote the population of dendritic cells, keratinocytes, T cells, and macrophages, respectively. The parameter β_IL36 represents the level of cytokine IL-36. The parameter k_D represents the population of dendritic cells in the lesion required to attain the half-maximal effect of dendritic cells (D) on D, K, and T (illustrated in the network 1). Similarly, k_K, k_T, and k_M denote the population of keratinocytes, T cells, and macrophages, respectively, which are required for the half-maximal effect of the K, T, and M populations on modeled cells. The parameters α_D, α_K, α_T, and α_M represent the rates of migration (in cells/day) of dendritic cells, keratinocytes, T cells, and macrophages, respectively, to the site of psoriatic lesion. Their rate of apoptosis (in day⁻¹) is represented by γ_D, γ_K, γ_T, and γ_M, respectively. The aforementioned set of ODEs indicates that the dynamics of the macrophage population are independent of the population of other cell types (D, K, and T). Hence, a quasi-steady state is assumed for macrophages, and the dynamics of the system of ordinary differential equations are governed by the following system of equations.

$\begin{array}{cc}\hfill \frac{\mathrm{d}D}{\mathrm{d}t}& ={\alpha }_D+{\beta }_{IL36}\left(\frac{D^n}{k_D^n+D^n}\frac{T^n}{k_T^n+T^n}\frac{M^{*n}}{k_M^n+M^{*n}}\frac{K^n}{k_K^n+K^n}\right)-{\gamma }_{D}D\hfill \\ \hfill \frac{\mathrm{d}K}{\mathrm{d}t}& ={\alpha }_K+{\beta }_{IL36}\left(\frac{D^n}{k_D^n+D^n}\frac{T^n}{k_T^n+T^n}\frac{M^{*n}}{k_M^n+M^{*n}}\frac{K^n}{k_K^n+K^n}\right)-{\gamma }_{K}K\hfill \\ \hfill \frac{\mathrm{d}T}{\mathrm{d}t}& ={\alpha }_T+{\beta }_{IL36}\left(\frac{D^n}{k_D^n+D^n}\frac{T^n}{k_T^n+T^n}\frac{M^{*n}}{k_M^n+M^{*n}}\frac{K^n}{k_K^n+K^n}\right)-{\gamma }_{T}T\hfill \end{array}(5)$

Here, $M^{*}=\frac{{\alpha }_M}{{\gamma }_M}$ denotes the macrophage population at a steady state. To further reduce the complexity of the system of ODEs, we make the system dimensionless by assuming $\bar{K}=\frac{K}{k_K}$, $\bar{T}=\frac{T}{k_T}$, $\bar{D}=\frac{D}{k_D}$, and $\bar{M}=\frac{M^{*}}{k_M}$. The dimensionless ODE system reads as

$\begin{array}{cc}\hfill \frac{\mathrm{d}\bar{D}}{\mathrm{d}t}& ={\alpha \prime }_D+{\beta \prime }_{IL36}\left(\frac{{\bar{D}}^n}{1+{\bar{D}}^n}\frac{{\bar{T}}^n}{1+{\bar{T}}^n}\frac{{\bar{M}}^n}{1+{\bar{M}}^n}\frac{{\bar{K}}^n}{1+{\bar{K}}^n}\right)-{\gamma }_D\bar{D}\hfill \\ \hfill \frac{\mathrm{d}\bar{K}}{\mathrm{d}t}& ={\alpha \prime }_K+{\beta \prime }_{IL36}\cdot e\cdot \left(\frac{{\bar{D}}^n}{1+{\bar{D}}^n}\frac{{\bar{T}}^n}{1+{\bar{T}}^n}\frac{{\bar{M}}^n}{1+{\bar{M}}^n}\frac{{\bar{K}}^n}{1+{\bar{K}}^n}\right)-{\gamma }_K\bar{K}\hfill \\ \hfill \frac{\mathrm{d}\bar{T}}{\mathrm{d}t}& ={\alpha \prime }_T+{\beta \prime }_{IL36}\cdot f\cdot \left(\frac{{\bar{D}}^n}{1+{\bar{D}}^n}\frac{{\bar{T}}^n}{1+{\bar{T}}^n}\frac{{\bar{M}}^n}{1+{\bar{M}}^n}\frac{{\bar{K}}^n}{1+{\bar{K}}^n}\right)-{\gamma }_T\bar{T}\hfill \end{array}(6)$

where ${\alpha \prime }_D=\frac{{\alpha }_D}{k_D}$, ${\beta \prime }_{IL36}=\frac{{\beta }_{IL36}}{k_D}$, ${\alpha \prime }_K=\frac{{\alpha }_K}{k_K}$, $e=\frac{k_D}{k_K}$, ${\alpha \prime }_T=\frac{{\alpha }_T}{k_T}$, and $f=\frac{k_D}{k_T}$.

2.3 Model parameters

The typical values of all the model parameters used for obtaining numerical solutions of the system of ODEs are summarized in Table 2. The values for the rate of apoptosis of T cells (γ_T), dendritic cells (γ_D), and keratinocytes (γ_K) are taken from a recent study (Pandey et al., 2021). The values of the remaining parameters were assumed subsequent to a systematic model exploration and calibration. Because the dependence of model results is explored through bifurcation analysis, the exact value of each parameter is not required.

Table 2

Table 2. Default values of the model parameters.

2.4 Sensitivity analysis

We performed a global sensitivity analysis using Sobol GSA Software (Kucherenko and Zaccheus, 2024) with MATLAB R2022b. Sobol’s method (Sobol, 1993) of sensitivity analysis quantifies the relative importance of input parameters by decomposing the total variance into contributions from individual parameters as well as combinations of parameters. Default sample values were generated using Sobol sequences according to parameter ranges, and their distributions are defined in Table 3. During this analysis, the software samples the parameter space simultaneously, considering all the parameters together, and reports two types of sensitivity index: first-order and total sensitivity. The first-order sensitivity indices are used to show the fractional contribution of a single parameter to the output variance, while higher-order sensitivity indices are calculated to measure the fractional contribution of parameter interactions along with the fractional contribution of individual model parameters to the output variance. In our results, we report the first-order sensitivity index to show the contribution of each model parameter to the total model output variance. Therefore, all model parameters are varied simultaneously in this approach, rather than one parameter at a time. The parameter ranges are chosen based on the observation that significant changes in the keratinocyte population were observed in the given range of parameters.

Table 3

Table 3. Parameters, their range of variation, and assumed distribution for global sensitivity analysis.

2.5 Bifurcation analysis

A bifurcation analysis is performed to investigate the steady-state behavior of the model and how it changes when parameter values are varied using an open-access tool, MATCONT (Govaerts et al., 2008), and MATLAB software (The MathWorks, Inc., 2022).

3 Results

We investigate the steady-state behavior of the network mediated by IL-36 (Figure 1) using a standard bifurcation analysis technique of dynamical systems theory. IL-36 is a proinflammatory cytokine, and its role in the progression of psoriasis has been observed (Mahil et al., 2017). Here, we study the effect of changes in the level of IL-36 cytokine on the population of keratinocytes for varying model parameters.

Our results suggest that an increase in the population of K and D would be observed for an increase in the level of the IL-36 (represented by ${\beta }_{IL-36}^{\prime }$) (Figure 2). A four-fold increase in the keratinocyte population is believed to be an indicator of a psoriatic condition, based on the observation by Valeyev et al. (2010). Therefore, a significant increase in the level of IL-36 would lead to a psoriatic state, and increased IL-36 levels have been observed in psoriasis (Mahil et al., 2017). In addition, a significant increase in the dendritic cell population has been observed in the pustule formation stage of GPP (Vyas et al., 2019). We observed that an increase in the value of the parameter e ($e=\frac{k_D}{k_K}$) would cause a sharp increase in the keratinocyte population, indicating the modulatory effect of keratinocytes on the other modeled cell types is stronger than the modulatory effect of dendritic cells, whereas a gradual increase in the dendritic cell population is observed for an increase in the values of parameter e (Figures 2A, B), which saturates for large values of e.

Figure 2

Figure 2. Steady-state response of the model to changes in different model parameters. (A) Change in the population of keratinocytes in response to varying model parameters. (B) Change in the population of dendritic cells in response to varying model parameters. Default parameters are: ${\alpha }_K^{\prime }$ = 0.1 day⁻¹, β_IL−36 = 4 days⁻¹, n = 1, $\bar{M}$ = 1, γ_K = 0.5 days⁻¹, γ_D = 0.2 days⁻¹, γ_T = 0.5 days⁻¹, ${\alpha }_D^{\prime }$ = 1 day⁻¹, e =2, ${\alpha }_T^{\prime }$ = 0.5 days⁻¹, f = 0.5. The initial condition is $\bar{K}=4,\bar{D}=2,\bar{T}=0.1$ The dashed line in the figure serves as a threshold for the keratinocyte population. Keratinocytes above this threshold are believed to indicate the psoriatic state (Valeyev et al., 2010).

In addition, the K and D populations would increase in the psoriatic lesion with an increase in the effective rate of migration of keratinocytes, dendritic cells, T cells, and macrophages (represented by ${\alpha }_K^{\prime }$, ${\alpha }_D^{\prime }$, ${\alpha }_T^{\prime }$, and ${\alpha }_{\bar{M}}^{\prime }$, respectively), as shown in Figure 2A.

Furthermore, an increase in the values of parameter f $(f=\frac{k_D}{k_T})$ would result in a gradual increase in the population of K and D, which saturates for further increase in f (Figures 2A, B). This suggests that the modulatory effect of the dendritic cell population on other modeled cells and itself is weak compared to that of T cells.

Our results also suggest an increase in the Hill-coefficient-like parameter n would lead to an initial increase in the keratinocyte and dendritic cell populations that saturate at large values of n.

For an increase in the steady-state population of macrophages$(\bar{M})$, the populations of keratinocyte and dendritic cells both increase gradually (Figure 2).

Through Figure 3, we demonstrate how the change in the population of keratinocytes upon varying the level of IL-36 (represented by β_IL−36) depends on other model parameters. Our model results suggest a sharp increase in the keratinocyte population for an increase in the effective rate of inward migration of keratinocytes $({\alpha }_K^{\prime })$, as depicted in Figure 3A. A switch-like increase would be observed for small values of ${\alpha }_K^{\prime }$ in Figure 3A. The bistable region bounded by the limit points would decrease with the increase in the ${\alpha }_K^{\prime }$. A sharp increase in the keratinocyte population would also be observed for increases in the effective rate of inward migration of dendritic cells $({\alpha }_D^{\prime })$ and the effective rate of inward migration of T cells $({\alpha }_T^{\prime })$ (Figure 3B and Supplementary Figure S1). In the case of a switch-like increase in the population of keratinocytes, the bistable region would decrease with an increase in ${\alpha }_D^{\prime }$ and ${\alpha }_T^{\prime }$ (Figure 3B and Supplementary Figure S1).

Figure 3

Figure 3. Keratinocyte population as a function of the level of IL-36 denoted by β′_IL−36. Change in keratinocyte population as a function of the level of IL-36 for different values of (A)${\alpha }_K^{\prime }$, (B) ${\alpha }_D^{\prime }$, (C) n,(D) f, (E) e, and (F) $\bar{M}$. Default parameter values are ${\alpha }_K^{\prime }$ = 0.1 day⁻¹, β_IL−36 = 4 days⁻¹, n = 2, $\bar{M}=1$, γ_K = 0.5 days⁻¹, γ_D = 0.2 days⁻¹, γ_T = 0.5 days⁻¹, ${\alpha }_D^{\prime }$ = 1 day⁻¹, e = 2, and ${\alpha }_T^{\prime }$ = 0.5 days⁻¹, and f = 0.5. The initial condition is $\bar{K}=4,\bar{D}=2,\bar{T}=0.1$.

An increase in the keratinocyte population would also be observed for an increase in the Hill-coefficient like parameter n (Figure 3C). For a small value of n, the population of keratinocytes would increase gradually. A switch-like increase is expected for a large value of n, and the bistable region would be larger for n = 3 than for n = 2.

Our results suggest gradual and switch-like increases in the population of keratinocytes for an increase in values of f (denoting the modulatory effect of keratinocytes on itself) and e (denoting the modulatory effect of dendritic cells on other modeled cell types, respectively (Figures 3D, E). Similar observations would be expected for an increase in the steady-state population of macrophages $(\bar{M})$ (Figure 3F).

3.1 Based on sensitivity analysis

Sensitivity indices for five parameters were obtained to assess their impact on the rate of change in the keratinocyte population (Figure 4) at different time points in the process of attaining the steady state. Our analysis suggests that the change in the parameter β_IL−36, which represents the level of IL-36 cytokines, was most influential over all time points for the change in keratinocyte cell population followed by the modulatory effects of keratinocytes on other cell populations and itself (e), the population of macrophages $(\bar{M})$, and Hill-coefficient like parameter (n), respectively. In addition, we observed that the change in keratinocyte population was unaffected by the change in modulatory effects of dendritic cells over the other cell population and itself (f).

Figure 4

Figure 4. Global sensitivity indices for keratinocyte population. (A) at the start (t = 0 days), (B) before the onset of steady state (t = 5 days), (C) at the onset of steady state (t = 10 days), and (D) post the onset of steady state (t = 20 days).

3.2 Based on bifurcation analysis

The switch-like increase in keratinocytes would be a robust behavior as we found large bistable regions in two-parameter space spanned by (n, ${\beta }_{IL-36}^{\prime }$) (Figure 5A). A typical temporal behavior of the network in bistable and monostable regions is shown in Figures 5B, C.

Figure 5

Figure 5. (A) Two-parameter bifurcation curve for n and ${\beta }_{IL-36}^{\prime }$. (B) Dynamics of the network for a model parameter combination that lies inside a bistable region, point b in Panel A of Figure 5. (C) Dynamics of the network for a model parameter combination that lies inside a monostable region at point a in Panel A. The values of parameters and initial condition are the same as in Figure 3.

4 Discussion

Psoriasis is characterized by the formation of a psoriatic lesion that is infiltrated by immune cell populations (Lowes et al., 2014; Dietrich et al., 2016). Several mathematical models have been developed to elucidate the role of cytokines in the progression of psoriasis (Roy et al., 2010; Valeyev et al., 2010; Roy and Datta, 2012). One of the models explored the roles of IL-15, IL-17/IL-23, and TNFα cytokines mediating indirect cell–cell interactions among keratinocytes, dendritic cells, and T cells in a psoriatic lesion (Pandey et al., 2021). These cytokines have been targeted in earlier studies that demonstrate their clinical potential.

Different types of psoriasis have been reported, and psoriasis vulgaris is the most prevalent type. GPP is another severe form of psoriasis. In GPP, relapsing sterile pustules are present in the affected individuals (Samotij et al., 2021). Recent studies attribute this relapse to T cells (skin resident memory TRM cell, and memory like γδT cells) whose proliferation and differentiation are caused by dendritic cells through a feed-forward loop involving inflammatory cytokines (reviewed in (Tiang and Lai, 2022)). The dermal dendritic cells colocalize with neutrophils that secrete elastase, which is believed to have a role in pustule formation (Skrzeczynska-Moncznik et al., 2012). Studies have shown that the elevated presence of human neutrophil elastase in GPP contributes to tissue damage, degradation of extracellular matrix, and disruption of the skin’s barrier function (Johnston et al., 2017). Hence, an increase in DC levels and their indirect interaction with elastase (Wittamer et al., 2005) could lead to worsening of psoriasis. Additionally, neutrophil elastase activates the cytokine IL-36γ, resulting in the exacerbation of psoriasis (Clancy et al., 2017).

Recent studies have unveiled IL-36 as a potential drug target treatment of GPP/psoriasis (Bachelez et al., 2019; Choon et al., 2021a; Choon et al., 2021b). This hypothesis is augmented by gene expression data that report a significant upregulation of the IL-36 cluster of cytokines in psoriasis patients (D’erme et al., 2015). To explore the role of proinflammatory cytokine IL-36 in the plaque formation and disease progression of psoriasis, we have built a network of indirect cell–cell interactions among key immune cells (dendritic cells, macrophages, and T-lymphocytes) and skin cells (keratinocytes) mediated by IL-36. Our results are consistent with the finding that IL-36 could lead to worsening of psoriasis. Furthermore, a sensitivity analysis for parameters (in Table 3) affecting the function of the network demonstrates that the cell populations are most sensitive to the level of IL-36. Our results suggest an increase in the population of keratinocytes; thus, the plaque formation and progression of psoriasis could occur either in a gradual or a switch-like manner. The switch-like-on set or progression of psoriasis would occur due to the bistable behavior of the network due to the indirect immune and keratinocyte cell interactions via IL-36. An exploration of the parameters controlling the bistable region suggests the switch-like abrupt changes in the keratinocyte population are a robust phenomenon. For instance, we found a sizeable bistable region in a two-parameter space spanned by the level of IL-36 $({\beta }_{IL-36}^{\prime })$ vs. the Hill-coefficient-like parameter n).

One limitation of our model is that it is formulated for only a single cytokine, whereas the role of several cytokines in psoriasis has been well-reported. Here, it is assumed that although several cytokines are involved in psoriasis, IL-36 has a dominant role, especially in GPP. Another limitation of the model is that all the IL-36 agonists (α, β, and γ) are represented as one single cytokine IL-36 instead of considering its isoforms separately. A future model could consider all the isoforms of the cytokine IL36, α, β, and γ, separately, and the role of different isoforms of IL-36 in the pathogenesis of psoriasis can be explored systematically. In addition, model results are yet to be experimentally verified. However, the present framework could also be used to study the role of a dominant cytokine in other inflammatory and autoimmune diseases.

This study highlights the network dynamics approach toward understanding the pathogenesis of psoriasis, which is a novel way to study cell–cell interactions in the context of human physiology. Several other studies have applied network dynamics to interpret organ interactions (Ivanov et al., 2021) and cortico-muscular interactions (Rizzo et al., 2020) and lay the groundwork in the field of network physiology and network medicine (Ivanov et al., 2016). Other models use fractional order differential equations to investigate the phenomenon of memory trace and study the influence of immune-boosting drugs on psoriasis (Özköse, 2024). This broadens the conventional approach of disease models framed by linear ODEs to include those featuring fractional order.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author contributions

SP: conceptualization, data curation, investigation, methodology, software, validation, writing–original draft, writing–review and editing. ST: data curation, formal analysis, investigation, methodology, resources, software, writing–original draft, writing–review and editing. SB: data curation, formal analysis, investigation, methodology, supervision, validation, writing–original draft, writing–review and editing. RM: conceptualization, formal analysis, project administration, supervision, writing–original draft, writing–review and editing. RP: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, validation, visualization, writing–original draft, writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. The funding numbers for grants acknowledged are as follows: IoE: R/Dev/D/IoE/Equipment/Seed Grant/ 2021-22/42671 UGC Startup: F.30-544/2021(BSR). The PMRF IDs of the scholars are as follows: ST: 0102940, SB: 0101840.

Acknowledgments

RP acknowledges the Department of Science and Technology, India, for the DST-INSPIRE Faculty Award (DST/INSPIRE/04/2015/001939) and Banaras Hindu University for the Institute of Eminence seed grant. RP also acknowledges the University Grant Commission, India, start-up grant awarded to him. SB and ST acknowledge the Prime Minister Research Fellowship awarded to them by the Ministry of Education Government of India.

Conflict of interest

Author RM was employed by Oligos Technologies Private Limited, India.

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

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnetp.2024.1363791/full#supplementary-material

References

Bachelez, H., Choon, S. E., Marrakchi, S., Burden, A. D., Tsai, T. F., Morita, A., et al. (2019). Inhibition of the interleukin-36 pathway for the treatment of generalized pustular psoriasis. N. Engl. J. Med. 380, 981–983. doi:10.1056/NEJMc1811317

PubMed Abstract | CrossRef Full Text | Google Scholar

Bridgewood, C., Fearnley, G. W., Berekmeri, A., Laws, P., Macleod, T., Ponnambalam, S., et al. (2018). IL-36γ is a strong inducer of IL-23 in psoriatic cells and activates angiogenesis. Front. Immunol. 9, 200. doi:10.3389/fimmu.2018.00200

PubMed Abstract | CrossRef Full Text | Google Scholar

Carrier, Y., Ma, H. L., Ramon, H. E., Napierata, L., Small, C., O’Toole, M., et al. (2011). Inter-regulation of Th17 cytokines and the IL-36 cytokines in vitro and in vivo: implications in psoriasis pathogenesis. J. Investigative Dermatology 131 (12), 2428–2437. doi:10.1038/jid.2011.234

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Du, Y., Han, Y., and Wei, Z. (2021). Correlation analysis between IL-35, IL-36 γ, CCL27 and psoriasis vulgaris. J. Dermatological Treat. 32 (6), 621–624. doi:10.1080/09546634.2019.1689226

CrossRef Full Text | Google Scholar

Choon, S. E., Lebwohl, M. G., Marrakchi, S., Burden, A. D., Tsai, T. F., Morita, A., et al. (2021a). Study protocol of the global Effisayil 1 Phase II, multicentre, randomised, double-blind, placebo-controlled trial of spesolimab in patients with generalized pustular psoriasis presenting with an acute flare. BMJ open. 11 (3), e043666. doi:10.1136/bmjopen-2020-043666

PubMed Abstract | CrossRef Full Text | Google Scholar

Choon, S. E., Marrakchi, S., Burden, A. D., Tsai, T. F., Morita, A., et al. (2021b). Trial of spesolimab for generalized pustular psoriasis. N. Engl. J. Med. 385, 2431–2440. doi:10.1056/nejmoa2111563

PubMed Abstract | CrossRef Full Text | Google Scholar

Clancy, D. M., Henry, C. M., Sullivan, G. P., and Martin, S. J. (2017). Neutrophil extracellular traps can serve as platforms for processing and activation of IL-1 family cytokines. FEBS J. 284 (11), 1712–1725. doi:10.1111/febs.14075

PubMed Abstract | CrossRef Full Text | Google Scholar

D’erme, A. M., Wilsmann-Theis, D., Wagenpfeil, J., Hölzel, M., Ferring-Schmitt, S., Sternberg, S., et al. (2015). IL-36γ (IL-1F9) is a biomarker for psoriasis skin lesions. J. Investigative Dermatology 135 (4), 1025–1032. doi:10.1038/jid.2014.532

CrossRef Full Text | Google Scholar

Dietrich, D., Martin, P., Flacher, V., Sun, Y., Jarrossay, D., Brembilla, N., et al. (2016). Interleukin-36 potently stimulates human m2 macrophages, langerhans cells and keratinocytes to produce pro-inflammatory cytokines. Cytokine 84, 88–98. doi:10.1016/j.cyto.2016.05.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Foster, A. M., Baliwag, J., Chen, C. S., Guzman, A. M., Stoll, S. W., Gudjonsson, J. E., et al. (2014). IL-36 promotes myeloid cell infiltration, activation, and inflammatory activity in skin. J. Immunol. 192 (12), 6053–6061. doi:10.4049/jimmunol.1301481

PubMed Abstract | CrossRef Full Text | Google Scholar

Govaerts, W., Kuznetsov, Y. A., Meijer, H., GaétanDhooge, E. A., and Sautois, B. (2008). New features of the software matcont for bifurcation analysis of dynamical systems. Math. Comput. Model. Dyn. Syst. 14 (2), 147–175. doi:10.1080/13873950701742754

CrossRef Full Text | Google Scholar

Günther, S., and Sundberg, E. J. (2014). Molecular determinants of agonist and antagonist signaling through the IL-36 receptor. J. Immunol. 193 (2), 921–930. doi:10.4049/jimmunol.1400538

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P.Ch., Liu, K. K. L., and Bartsch, R. P. (2016). Focus on the emerging new fields of network physiology and network medicine. New J. Phys. 18 (100201), 100201. doi:10.1088/1367-2630/18/10/100201

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P.Ch., Wang, J. W. J. L., and Zhang, X. (2021). “Signal processing in network physiology: quantifying network dynamics of organ interactions,” in 2020 28th European Signal Processing Conference (EUSIPCO), Amsterdam, Netherlands, 18-21 Jan. 2021, 945–949.

CrossRef Full Text | Google Scholar

Johnston, A., Xing, X., Wolterink, L., Barnes, D. H., Yin, Z., Reingold, L., et al. (2017). IL-1 and IL-36 are dominant cytokines in generalized pustular psoriasis. J. Allergy Clin. Immunol. 140 (1), 109–120. doi:10.1016/j.jaci.2016.08.056

PubMed Abstract | CrossRef Full Text | Google Scholar

Kucherenko, S., and Zaccheus, O. (2024). Sobolgsa software. Imperial College. Available at: https://www.imperial.ac.uk/process-systems-engineering/research/free-software/sobolgsa-software/.

Google Scholar

Li, N., Yamasaki, K., Saito, R., Fukushi-Takahashi, S., Shimada-Omori, R., Asano, M., et al. (2014). Alarmin function of cathelicidin antimicrobial peptide LL37 through IL-36γ induction in human epidermal keratinocytes. J. Immunol. 193 (10), 5140–5148. doi:10.4049/jimmunol.1302574

PubMed Abstract | CrossRef Full Text | Google Scholar

Lowes, M. A., Suárez-Fariñas, M., and Krueger, J. G. (2014). Immunology of psoriasis. Annu. Rev. Immunol. 32, 227–255. doi:10.1146/annurev-immunol-032713-120225

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahil, S. K., Catapano, M., Di Meglio, P., Dand, N., Ahlfors, H., Carr, I. M., et al. (2017). An analysis of IL-36 signature genes and individuals with IL1RL2 knockout mutations validates IL-36 as a psoriasis therapeutic target. Sci. Transl. Med. 9 (411), eaan2514. doi:10.1126/scitranslmed.aan2514

PubMed Abstract | CrossRef Full Text | Google Scholar

Mercurio, L., Morelli, M., Scarponi, C., Eisenmesser, E. Z., Doti, N., Pagnanelli, G., et al. (2018). IL-38 has an anti-inflammatory action in psoriasis and its expression correlates with disease severity and therapeutic response to anti-IL-17A treatment. Cell death Dis. 9 (11), 1104–1113. doi:10.1038/s41419-018-1143-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Oza, H. B., Pandey, R., Roper, D., Al-Nuaimi, Y., Spurgeon, S. K., and Goodfellow, M. (2017). Modelling and finite-time stability analysis of psoriasis pathogenesis. Int. J. Control 90 (8), 1664–1677. doi:10.1080/00207179.2016.1217566

CrossRef Full Text | Google Scholar

Özköse, F. (2024). Modeling of psoriasis by considering drug influence: a mathematical approach with memory trace. Comput. Biol. Med. 168, 107791. doi:10.1016/j.compbiomed.2023.107791

PubMed Abstract | CrossRef Full Text | Google Scholar

Pandey, R., Al-Nuaimi, Y., Mishra, R. K., Spurgeon, S. K., and Goodfellow, M. (2021). Role of subnetworks mediated by TNFα, IL-23/IL-17 and IL-15 in a network involved in the pathogenesis of psoriasis. Sci. Rep. 11 (1), 1–13.

PubMed Abstract | CrossRef Full Text | Google Scholar

Pandey, R., Jangid, A., Vinjamuri, R. G., and Ramaswamy, R. (2023). Modelling of indirect cell–cell interaction networks mediated by IFNγ/IL-4 cytokine involved in atopic dermatitis. J. Theor. Biol. 556, 111291. doi:10.1016/j.jtbi.2022.111291

PubMed Abstract | CrossRef Full Text | Google Scholar

Parisi, R., Symmons, D. P., Griffiths, C. E., and Ashcroft, D. M. (2013). Global epidemiology of psoriasis: a systematic review of incidence and prevalence. J. Investigative Dermatology 133 (2), 377–385. doi:10.1038/jid.2012.339

PubMed Abstract | CrossRef Full Text | Google Scholar

Rizzo, R., Zhang, X., Wang, J. W. J. L., Lombardi, F., and Ivanov, P.Ch. (2020). Network physiology of cortico–muscular interactions. Front. Physiology 11, 558070. doi:10.3389/fphys.2020.558070

CrossRef Full Text | Google Scholar

Roy, P. K., Bhadra, J., and Chattopadhyay, B. (2010). Mathematical modeling on immunopathogenesis in chronic plaque of psoriasis: a theoretical study. Proc. World Congr. Eng. 1, 550–555.

Google Scholar

Roy, P. K., and Datta, A. (2012). Negative feedback control may regulate cytokines effect during growth of keratinocytes in the chronic plaque of psoriasis: a mathematical study. Int. J. Appl. Math. 25 (2), 233–254.

Google Scholar

Sachen, K. L., Greving, C. N. A., and Towne, J. E. (2022). Role of IL-36 cytokines in psoriasis and other inflammatory skin conditions. Cytokine 156, 155897. doi:10.1016/j.cyto.2022.155897

PubMed Abstract | CrossRef Full Text | Google Scholar

Samotij, D., Szczęch, J., and Reich, A. (2021). Generalized pustular psoriasis: divergence of innate and adaptive immunity. Int. J. Mol. Sci. 22 (16), 9048. doi:10.3390/ijms22169048

PubMed Abstract | CrossRef Full Text | Google Scholar

Skrzeczynska-Moncznik, J., Wlodarczyk, A., Zabieglo, K., Kapinska-Mrowiecka, M., Marewicz, E., Dubin, A., et al. (2012). Secretory leukocyte proteinase inhibitor-competent dna deposits are potent stimulators of plasmacytoid dendritic cells: implication for psoriasis. J. Immunol. 189 (4), 1611–1617. doi:10.4049/jimmunol.1103293

PubMed Abstract | CrossRef Full Text | Google Scholar

Sobol, I. M. (1993). Sensitivity estimates for nonlinear mathematical models. Math. Model. Comput. Exp. (4), 407–414.

Google Scholar

The MathWorks, Inc. (2022). MATLAB version: 9.13.0 (R2022b). Available at: https://www.mathworks.com (Accessed January 01, 2023).

Google Scholar

Tiang, D., and Lai, Y. (2022). The relapse of psoriasis: mechanisms and mysteries. J. Investigative Dermatology Innovations 2, 100116. doi:10.1016/j.xjidi.2022.100116

CrossRef Full Text | Google Scholar

Valeyev, N. V., Hundhausen, C., Umezawa, Y., Kotov, N. V., Williams, G., Clop, A., et al. (2010). A systems model for immune cell interactions unravels the mechanism of inflammation in human skin. PLoS. Comput. Biol. 6 (12), 10010244–e1001112. doi:10.1371/journal.pcbi.1001024

PubMed Abstract | CrossRef Full Text | Google Scholar

Vigne, S., Palmer, G., Lamacchia, C., Martin, P., Talabot-Ayer, D., Rodriguez, E., et al. (2011). IL-36R ligands are potent regulators of dendritic and t cells. Blood, J. Am. Soc. Hematol. 118 (22), 5813–5823. doi:10.1182/blood-2011-05-356873

PubMed Abstract | CrossRef Full Text | Google Scholar

Vyas, N. S., Charifa, A., Desman, G. T., and McNiff, J. M. (2019). Distinguishing pustular psoriasis and acute generalized exanthematous pustulosis on the basis of plasmacytoid dendritic cells and mxa protein. J. Cutan. Pathology 46 (5), 317–326. doi:10.1111/cup.13430

PubMed Abstract | CrossRef Full Text | Google Scholar

Wittamer, V., Bondue, B., Guillabert, A., Vassart, G., Parmentier, M., and Communi, D. (2005). Neutrophil-mediated maturation of chemerin: a link between innate and adaptive immunity. J. Immunol. 175 (1), 487–493. doi:10.4049/jimmunol.175.1.487

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L. J., Sen, G. L., Ward, N. L., Johnston, A., Chun, K., Chen, Y., et al. (2016). Antimicrobial peptide LL37 and mavs signaling drive interferon-β production by epidermal keratinocytes during skin injury. Immunity 45 (1), 119–130. doi:10.1016/j.immuni.2016.06.021

PubMed Abstract | CrossRef Full Text | Google Scholar