A mathematical study for psoriasis transmission with immune-mediated time delays and optimal control strategies

TL;DR

This paper uses a mathematical model to study how immune delays affect psoriasis and how biologic treatments can help manage the disease.

Contribution

The study introduces a novel mathematical model with immune-mediated time delays and optimal control strategies for psoriasis.

Findings

Immune-mediated delays can cause instability and oscillations in psoriasis dynamics through Hopf bifurcations.

Optimal control strategies using TNF-α and IL-17 inhibitors can stabilize the immune response and reduce keratinocyte density.

Including time delays in the model provides a more realistic understanding of psoriasis progression and treatment.

Abstract

In psoriasis, dendritic cells activate T cells, which then release excessive pro-inflammatory cytokines, leading to abnormal growth of keratinocytes in the epidermis. At the same time, anti-inflammatory cytokines attempt to restore balance. In reality, these immune processes are not immediate; they involve biological time gaps due to signal processing, cell communication, and cytokine feedback. Such immune-related delays may play a key role in triggering unstable or oscillatory behavior observed in psoriasis flare-ups. In this study, we present and analyze a mathematical model of psoriasis that explicitly includes two intracellular immune-mediated time delays to demonstrate their biological significance in disease progression. The model captures the interactions among T cells, dendritic cells, keratinocytes, and local mature stem cells. It features two cytokine-mediated feedback loops between T cells and dendritic cells, while stem cells attempt to regulate the immune response through anti-inflammatory signaling. A key challenge is identifying the critical time delays that modulate these interactions. To address this, we introduce two different delays in different interaction terms of the model system. We test the hypothesis that these delays can critically influence the onset and persistence of psoriatic pathology mathematically. Using stability analysis of the interior equilibrium, we determine parametric relations, their ranges, and delay thresholds that give rise to Hopf bifurcations, thereby linking delays to disease and deriving conditions of instability. Our analysis demonstrates that both immune-mediated delays critically influence system stability, with threshold values of τ1* and τ2* inducing oscillations through Hopf bifurcations. Further, we apply optimal control strategies on the delayed system using the effects of two biologic agents: TNF-α and IL-17 inhibitors. Incorporation of optimal controls effectively stabilizes the immune response. Numerical simulations support these analytical findings and show that biologic interventions can effectively reduce keratinocyte density. Inclusion of immune-related delays, based on both analytical and numerical results, provides a more realistic understanding of psoriasis dynamics and helps optimize therapeutic approaches for psoriasis management.

Figures

Schematic representation of model populations and their interrelationships mediated by cytokine signaling.

The dynamics of the model populations, namely T cells, dendritic cells, keratinocytes, and mesenchymal stem cells, for the delayed system with τ1=τ2=1.2.

This figure illustrates the populations of three cell types T cells, dendritic cells, and keratinocytes, for both the non-delayed system (τ1=τ2=0, shown with blue trajectories) and the delayed system (τ1=1.2,τ2=1, shown with orange trajectories).The introduction of non-zero delays results in periodic solutions, with the corresponding limit cycles depicted in a 3D plot.

In sub-figure (a), the populations of T cells, dendritic cells, and keratinocytes are shown for the delay parameter τ2=0, with trajectories plotted for two distinct values of τ1 (0.8 and 1), each represented by a different color code.Sub-figure (b) depicts the Hopf bifurcation with respect to the delay parameter τ1 while keeping τ2=0.

In Sub-figure (a), the populations of T cells, dendritic cells, and keratinocytes are shown with a fixed delay of τ1=0.The trajectories are plotted for two distinct values of τ2 (0.85 and 1.2), each represented by a different color. Sub-figure (b) illustrates the Hopf bifurcation with respect to the delay parameter τ2, while the other delay parameter, τ1, is fixed at 0.

The stability region of the endemic equilibrium point is shown in the τ1−τ2 plane.The color coding denotes the max[Re(λ)], where λ is the eigenvalues corresponding to the delayed system.

Numerical solution of the optimal control problem with τ1=τ2=1.1 and remaining parameter values taken from Table 1.The red dotted trajectories represent the disease dynamics without optimal control, while the blue trajectories depict the population behavior under the influence of optimal control.

Numerical approximation of Pontryagin’s extremals and the corresponding optimal control profiles v1*(t) and v2*(t).

PRCC scatter plots for all system parameters against the keratinocyte population, with a sample size of N = 1000.The PRCC values are calculated at 100 days and all p values are taken as up to four decimal places. The x-axis shows the residuals obtained from regressing the rank-transformed values of the parameter of interest on the rank-transformed values of all other parameters. The y-axis displays the residuals from regressing the rank-transformed keratinocyte values on the rank-transformed values of all other parameters, excluding the one being analyzed.

Equations

Full text

1 Introduction

Psoriasis is a chronic autoimmune skin disease that affects around 2–3% of the global population [1]. It is mainly caused by abnormal immune regulation and excessive production of pro-inflammatory cytokines, which lead to rapid growth of keratinocytes and the formation of thick, scaly skin patches [2,3]. The disease not only reduces the physical quality of life but also creates psychological and social stress for patients [4]. Environmental factors such as smoking and alcohol consumption further increase the severity of psoriasis [5].

The interaction between T cells and dendritic cells (DCs) is considered the central mechanism in the development of psoriasis. Dysfunctional DCs activate T cells by releasing cytokines such as TNF-α and IL-23, while activated T cells further stimulate DCs through feedback signals including TNF-α and IL-17 [6,7]. This creates a self-sustaining loop that amplifies inflammation and promotes abnormal keratinocyte proliferation [8]. Although mesenchymal stromal cells (MSCs) attempt to regulate this imbalance by controlling cytokine release, they often lose their protective role during psoriatic lesions [9]. Thus, psoriasis progression can be seen as the outcome of multiple interacting immune processes, sustained by feedback mechanisms.

Several researchers have studied psoriasis using mathematical modeling to capture such complex biological interactions. Oza et al. [10] examined convergence properties in cytokine-driven systems, Roy et al. [11] analyzed cytokine signaling pathways and biologic control strategies. Kushary et al. [12] recently proposed a mathematical model for psoriasis. Later, they studied the fractional form of this model using three operators, namely Caputo, Caputo–Fabrizio, and Atangana–Baleanu in the Caputo sense [13].

Delays naturally arise in immune systems because cytokine signaling, cellular activation, infiltration, and subsequent responses are not instantaneous. Inspired by the work of Das et al. [14], who highlighted the role of delays in immune system modeling, we extend our previously proposed psoriasis model [12] by investigating the influence of immune-mediated delays. In the context of psoriasis, such delays can significantly alter disease dynamics by destabilizing the system or inducing oscillatory immune responses. To explore this, we formulate a delay-based mathematical model in which intracellular time delays are incorporated into the separate interaction terms of T cells and DCs mediated through cytokine signaling. We analyze the impact of these delays on disease progression, perform delay-induced Hopf bifurcation analysis, and design optimal control strategies on the delayed model using biologic inhibitors to regulate the system. Numerical simulations are further carried out to validate the theoretical findings. The novelty of this study lies in demonstrating how immune-mediated delays affect the stability of the psoriasis model and how biologic control strategies can be optimized within this delay-induced framework. We expect that the results of this work will provide useful insights for designing improved control strategies for psoriasis.

2 The model introducing immune delays

In this section, we present a mathematical model for psoriasis in the presence of immune-mediated delays. In our earlier work, Kushary et al. [12] formulated a four-dimensional deterministic model based on ordinary differential equations (ODEs). The detailed assumptions underlying that model are given in the cited reference. The system considered the main cellular populations involved in psoriasis, namely T cells, dendritic cells, keratinocytes, and mesenchymal stem cells, under appropriate biological assumptions. The present study extends that earlier model by incorporating time delays.

The interactions between T cells and dendritic cells, mediated by pro-inflammatory cytokines, are not instantaneous. Anti-inflammatory cytokines attempt to counterbalance these signals, which leads to a lag in immune activation and inflammatory cell infiltration. To capture this effect, we introduce two distinct intracellular immune response delays into the previous model. These delays represent the time gap in the interaction between T cells and dendritic cells, and their subsequent influence on keratinocyte proliferation.

Based on these considerations, the delayed mathematical model is formulated as follows:

$$[eqn]$$

subject to the initial conditions:

$$[eqn]$$

Here, T(t), D(t), K(t), and S(t) denote the population densities of T cells, dendritic cells, keratinocytes, and mesenchymal stem cells (MSCs), respectively, at time t. The interaction network among these cell populations is shown in the schematic diagram (Fig 1), and the model parameters are described in detail in Table 1. The parameter τ represents the time delay, measured in days.

Schematic representation of model populations and their interrelationships mediated by cytokine signaling.

Table 1: Parameter values used in the numerical simulations for the entire model system.Several parameter ranges were obtained from previous studies [12,24,25]. Due to the lack of sufficient primary data related to immune-mediated delay factors, some parameter values were chosen to ensure the biological plausibility of the model behavior.

Let $[eqn]$ be the Banach space of all continuous functions $[eqn]$ , equipped with the usual supremum norm defined by:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ with $[eqn]$ . To make the model feasible from the biological point of view, initial functions are assumed to be $[eqn]$ 0 for $[eqn]$ and $[eqn]$ for $[eqn]$ . Using the fundamental theory of functional differential equations [15], we can guarantee the uniqueness of the solutions of the system (1) with the initial conditions mentioned above.

3 Mathematical properties of the delayed model

3.1 Positivity and boundedness

Positivity and boundedness in a mathematical model are essential to establish well-posedness. It provides a reliable basis for the prediction of the system and gives meaningful insight into the outcome of the model. In this subsection, we analyze the positivity and boundedness of the delay-induced system (1) through the following theorem.

Theorem 1. Consider system (1) with initial history $[eqn]$ such that $[eqn]$ for $[eqn]$ . Then:

Every solution remains positive for all $[eqn]$ , i.e., $[eqn]$ . Hence, $[eqn]$ is a positively invariant region for (1).Every positive solution is bounded and ultimately enters the domain of attraction $[eqn]$ , where

$$[eqn]$$

with $[eqn]$ and $[eqn]$ .

Proof: Let $[eqn]$ and

$$[eqn]$$

Then system (1) can be written as

$$[eqn]$$

1. Let $[eqn]$ . From the form of F it follows that

$$[eqn]$$

By Lemma 2 of Yang et al. [16] and Theorem 1.1 of Bodnar et al. [17], any solution X(t) of (2) with $[eqn]$ satisfies $[eqn]$ for all $[eqn]$ . Therefore, $[eqn]$ is a positively invariant region for (1).

2. Adding the first three equations of (1) gives

$$[eqn]$$

where $[eqn]$ and $[eqn]$ . The quadratic term satisfies $[eqn]$ . Hence,

$$[eqn]$$

By the comparison principle [18], for t > 0,

$$[eqn]$$

so for large t,

$$[eqn]$$

For the fourth equation we have the explicit solution

$$[eqn]$$

Therefore, every positive solution ultimately enters and remains in

$$[eqn]$$

which satisfies $[eqn]$ .

3.2 Feasibility of positive equilibrium

To determine the interior equilibrium point E (T , D , K , S ) of the delay-induced system, we set all the equations of the system (1) to zero and solve for the state variables under the conditions T  > 0, D  > 0, K  > 0, and S  > 0. From the fourth equation of the system (1), we get $[eqn]$ Using the value of S , we introduce two constants, $[eqn]$ and $[eqn]$ , to simplify the computations. These constants are defined as follows

$$[eqn]$$

Subsequently, we determined the values of D  and K  in terms of T  using the aforementioned constants and by setting the second equation of the system (1) to zero. These values are given by

$$[eqn]$$

By setting the first equation of the system (1) to zero and substituting the obtained values of the state variables into this equation, we derive a cubic equation in terms of T .

$$[eqn]$$ $$[eqn]$$

Using the parameter values listed in Table 1, we observe that P1>0 is positive, while P2<0 and P3<0. According to Descartes’ rule of signs [19], the system (5) has exactly one positive root. By substituting the value of T  into (4), we can determine the interior equilibrium point E (T , D , K , S ).

3.3 Stability and Hopf-bifurcation analysis

The non-delayed system admits a unique positive interior equilibrium point, since the extinction of any population is biologically unrealistic. Without loss of generality, we denote this equilibrium as E (T , D , K , S ). The exact values of the state variables at E  were determined earlier in our previous work by Kushary et al. [12].

In this subsection, we study the local stability of the delay-induced system (1) in a neighborhood of the interior equilibrium point E .

Let $[eqn]$ and $[eqn]$ represent the perturbed variables near the equilibrium point E . Linearizing the delay-induced system (1) about E , we derive the following system

$$[eqn]$$

where $[eqn]$ and $[eqn]$ , and $[eqn]$ represent the 4 × 4 Jacobian matrices of the system evaluated at E  for the current state, the delayed state with delay $[eqn]$ , and the delayed state with delay $[eqn]$ , respectively. Their explicit forms are

$$[eqn]$$

Now, the characteristic equation of the system (1) is given by

$$[eqn]$$

where the coefficients are

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The interior steady state E  is said to be locally asymptotically stable (LAS) if all roots λ of the characteristic Eq (8) are negative or possess negative real parts. The following cases will now be analyzed in detail.

Case - I ( $[eqn]$ = $[eqn]$ = 0)

In this case, the characteristic Eq (8) simplifies to the following form

$$[eqn]$$

One of the eigenvalues is $[eqn]$ . Thus, the interior equilibrium E  is LAS in the absence of delays if following conditions are satisfied

$$[eqn]$$

Case - II ( $[eqn]$ and $[eqn]$ = 0)

For this case, the characteristic Eq (8) takes the form

$$[eqn]$$

The root $[eqn]$ is negative. Thus, we are only interested in analyzing the following equation

$$[eqn]$$

Since, the characteristic Eq (11) is transcendental in λ, the stability of the system (1) cannot be analyzed using the classical Routh-Hurwitz criteria [20]. The Eq (11) admits purely imaginary solutions as a necessary condition for a stability change at the interior equilibrium E . Let $[eqn]$ be a root of Eq (11). Substituting $[eqn]$ into Eq (11) and separating the real and imaginary components, we obtain

$$[eqn]$$

Squaring and adding the above two equations and substituting $[eqn]$ , we get the following equation

$$[eqn]$$

The Eq (11) does not admit purely imaginary roots if the Eq (13) satisfies the Routh-Hurwitz’s criteria [20]. The results are summarized in the following proposition.

Proposition 1. The interior equilibrium E  is LAS for all $[eqn]$ if the following conditions are hold

$$[eqn]$$

If R3<0, the Eq (13) has at least one positive root. Let $[eqn]$ denote the smallest positive root. In this case, the characteristic Eq (11) will have purely imaginary roots, says $[eqn]$ . We now determine the critical value of $[eqn]$ , denoted as $[eqn]$ that corresponds to the point at which the stability of the interior equilibrium E  changes. The critical value $[eqn]$ is determined from the system of Eq (12) to analyze the conditions under which the system transitions from stability to instability as $[eqn]$ increases, as follows

$$[eqn]$$

Let us assume that for $[eqn]$ , $[eqn]$ . Based on this assumption and summarizing the above discussions, we have derived the following theorem.

Theorem 2. If R3<0*, according to Butler lemma [21] the interior equilibrium E*  is LAS when $[eqn]$ and becomes unstable when $[eqn]$ , where

$$[eqn]$$

Furthermore, when $[eqn]$ , a Hopf bifurcation [22] occurs, meaning that a family of periodic solutions of the system (1) bifurcates as $[eqn]$ crosses the critical value $[eqn]$ . This bifurcation is ensured by the transversality condition

$$[eqn]$$

Proof: Now, considering λ as a function of $[eqn]$ , i.e., $[eqn]$ , and differentiating Eq (11) with respect to $[eqn]$ , we derive the following expression

$$[eqn]$$

From this above equation, we get

$$[eqn]$$

Hence, we conclude that the delayed system (1) undergoes a Hopf bifurcation at the critical value $[eqn]$ if the transversality condition $[eqn]$ i.e., when $[eqn]$ is satisfied.

Case - III ( $[eqn]$ and $[eqn]$ )

In this case, the characteristic Eq (8) reduces to

$$[eqn]$$

Following a similar approach as for the system described in (10), it can be shown that a critical value $[eqn]$ exists for the positive interior equilibrium point E , at which the stability transitions for the delayed system (1). Without repeating the detailed analytical calculations presented in Case - II, we provide the general expression that characterizes the transition from stability to instability as the value of $[eqn]$ increases.

$$[eqn]$$

Assuming that for $[eqn]$ , the sequence $[eqn]$ represents the critical values of $[eqn]$ , where $[eqn]$ . The minimum value in this sequence, $[eqn]$ , is identified as the critical point. The results are summarized in the following theorem without proving it.

Theorem 3. The interior equilibrium E  is LAS if $[eqn]$ and becomes unstable if $[eqn]$ , where

$$[eqn]$$

Moreover, the delayed system (1) undergoes a Hopf bifurcation at $[eqn]$ when transversality condition $[eqn]$ hold, where $[eqn]$ and $[eqn]$ .

Case - IV ( $[eqn]$ and $[eqn]$ )

The analysis of this case is extensive and complex, making it challenging to derive precise information regarding the nature of the eigenvalues and the conditions under which stability switches occur. Instead of performing detailed analytical calculations, we have numerically demonstrated the emergence of Hopf-bifurcating periodic solutions for the system (1) and provided the necessary biological interpretations. The result, stated without proof, is presented in the following theorem.

Theorem 4. The non-delayed system is asymptotically stable under the given conditions. However, there exists a critical threshold $[eqn]$ beyond which the stability of the system undergoes a change. Specifically

The steady state E  remains LAS when $[eqn]$ .The steady state E  becomes unstable when $[eqn]$ .A Hopf bifurcation occurs at $[eqn]$  +  $[eqn]$ , provided the transversality condition $[eqn]$ is satisfied.

4 Optimal control strategy on delayed system

The incorporation of optimal control approaches with delays has become a prominent area of research, particularly due to its applications in biological control systems. In the context of autoimmune responses in psoriasis, the cyclic interactions between T cells and dendritic cells are significantly amplified by the effects of two key cytokines: TNF-α and IL-17. This interaction leads to excessive growth of keratinocytes, which is a key characteristic of the disease. To mitigate this, our objective is to determine optimal control strategies for two biologic treatments that minimize the overall treatment cost of psoriasis. We address this by incorporating two control functions, $[eqn]$ and $[eqn]$ into the delayed system (1).

$[eqn]$ : Represents the effect of a TNF-α inhibitor, which can modulate the reaction rates between T cells and dendritic cells. $[eqn]$ : Represents the inhibition of the IL-17 axis, targeting its associated pathways.

Thus, the delayed system (1), augmented with optimal control functions, is described by the following set of equations over the specified time interval [0,t_f_].

$$[eqn]$$

with $[eqn]$ and $[eqn]$ .

4.1 Description of the objective functional

Our primary objective is to minimize the cost associated with each biologic treatment to reduce the overall treatment expense of psoriasis. To achieve this, we define an objective cost functional for the minimization problem, expressed as follows:

$$[eqn]$$

subject to the control system (21). The parameters L1 and L2 represent positive weight constants on the benefit of the cost of the TNF-α and IL-17 inhibitors, respectively. Our aim is to determine the optimal control $[eqn]$ that satisfies the following:

$$[eqn]$$

The Lebesgue measurable control set is defined on the interval [0,t_f], where tf_ represents the terminal time of control.

$$[eqn]$$

Pontryagin’s Minimum Principle with delays has been applied to establish the necessary conditions for solving this optimal control problem [23]. The optimal control-induced system (21) has been shown to admit non-negative, bounded solutions for a bounded Lebesgue-measurable control function and the non-negative initial conditions.

4.2 Optimal control’s characteristics

To characterize the optimal control, Pontryagin’s Minimum Principle reduces the problem (21)–(24) to a problem of minimizing the Hamiltonian H, defined as

$$[eqn]$$

where, $[eqn]$ ’s $[eqn]$ are the adjoint variables to be determined suitably. The adjoint system with transversality criteria $[eqn]$ for $[eqn]$ can be obtained by using Pontryagin’s Minimum Principle [23] as

$$[eqn]$$

The optimality of the considered system consists of the optimal states and the corresponding adjoint system takes the form

$$[eqn]$$

We determine the optimal controls $[eqn]$ ( $[eqn]$ ), using the optimality condition provided by

$$[eqn]$$

Taking partial differentiation of the Eq (25) with respect to $[eqn]$ and $[eqn]$ separately and applying the condition (27), we obtain

$$[eqn]$$ $$[eqn]$$

Now, using the conventional control’s boundedness criteria and following the characteristics of the control set $[eqn]$ that the admissible control takes the values such that $[eqn]$ , we can have

$$[eqn]$$

From this result, we can derive the following theorem.

Theorem 5. If the objective cost functional $[eqn]$ attains its minimum value for optimal controls $[eqn]$ and $[eqn]$ , moreover (T , D , K , S ) be the corresponding optimal state for the optimal control problem (21), then there exist adjoint functions $[eqn]$ satisfying the transversality conditions $[eqn]$ , for $[eqn]$ . Furthermore, the optimal control solutions are given by

$$[eqn]$$

5 Numerical simulations

In this section, we perform several numerical simulations to validate our analytical results for both the delayed system (1) and the optimal control system on the delay-induced model (21). The initial values for the model populations are chosen to satisfy the biological assumptions underlying this study. To demonstrate the behavior of the model, numerical simulations have been performed using the parameter values provided in Table 1.

5.1 Simulations of the stability of equilibrium

The initial values of all model populations for the simulation are chosen as follows

$$[eqn]$$

Initially, when both delay parameters are set to $[eqn]$ , representing the non-delayed system, the model populations, namely, T cells [T(t)], dendritic cells [D(t)], keratinocytes [K(t)], and mesenchymal stem cells [S(t)], exhibit asymptotically stable behavior. However, introducing time delays $[eqn]$ and $[eqn]$ (see Fig 2), all system populations, except mesenchymal stem cells [S(t)], show oscillatory dynamics. For this simulation, the intracellular delay parameters are set to $[eqn]$ , while the remaining parameter values are taken from Table 1.

The dynamics of the model populations, namely T cells, dendritic cells, keratinocytes, and mesenchymal stem cells, for the delayed system with τ1=τ2=1.2.

Since introducing non-zero delays into the system (1) does not alter the asymptotic behavior of the MSC population, our focus shifts to the dynamics of the other state populations. In (Fig 3), we illustrate the behavior of the first three cell types: T cells, dendritic cells, and keratinocytes, for both non-delayed system ( $[eqn]$ ), represented by blue trajectories and the delayed system ( $[eqn]$ ), represented by orange trajectories. Without delays, all system populations are asymptotically stable and converge to the endemic steady state. The introduction of non-zero delays exhibits stable periodic oscillations in the populations of these cell types. Moreover, the resulting limit cycles are depicted in a 3-D plot with respect to these populations, providing a comprehensive visualization of the periodic behavior. The simulation demonstrates the significant effect of intracellular delays on the autoimmune response for psoriasis transmission for the considered model system.

This figure illustrates the populations of three cell types T cells, dendritic cells, and keratinocytes, for both the non-delayed system (τ1=τ2=0, shown with blue trajectories) and the delayed system (τ1=1.2,τ2=1, shown with orange trajectories).The introduction of non-zero delays results in periodic solutions, with the corresponding limit cycles depicted in a 3D plot.

5.2 Simulations with different time delays

In this subsection, we examine numerically the influence of time delays on first three model populations of the system (1) and perform some simulations in order to depicted the delay induced Hopf bifurcating positions by considering each delay parameters separately at a time.

If the time delay parameter $[eqn]$ is fixed to zero ( $[eqn]$ ), and $[eqn]$ and $[eqn]$ respectively, we observe in Fig 4(a) that the oscillation increases as $[eqn]$ increases and periodic oscillation is seen at $[eqn]$ day. For $[eqn]$ , the cell densities initially oscillate but eventually achieve local asymptotic stability, confirming our analytical finding stated in Case-II in the Sect 3.3. The trajectories for $[eqn]$ and $[eqn]$ are illustrated using two distinct color codes. In this scenario, the stability of the endemic equilibrium E  depends on the value of the time delay $[eqn]$ , where $[eqn]$ is fixed to zero. Larger values of $[eqn]$ lead to regular oscillatory solutions, indicating that an increased time delay between the interactions of T cells with dendritic cells and keratinocytes results in sustained oscillations. Stability switch occurs through Hopf bifurcation which are depicted in the Fig 4(b) where we present the delay-induced Hopf-bifurcation, focusing on T cells, dendritic cells, and keratinocyte cells populations only. In this analysis, $[eqn]$ serves as the bifurcation parameter. The parameter $[eqn]$ is varied within the range of 0.7 to 1. The figure illustrates that, for $[eqn]$ , there exists a critical value $[eqn]$ at which the endemic equilibrium point E  transitions from stability to instability. This observation highlights the critical role of the delay parameter $[eqn]$ in determining the stability of the system (1). Biologically, this implies that an increase in the time lag within the interaction term of the T cell growth equation—mediated by dendritic cell mediated cytokines infiltrating and accelerates fluctuations in the model cell concentrations, excluding the MSC population. These fluctuations manifest as oscillatory dynamics and this behavior arises from the inflammatory response in the upper epidermal layer caused by the over-expression of keratinocytes, which, in turn, amplifies the progression and dissemination of the disease.

In sub-figure (a), the populations of T cells, dendritic cells, and keratinocytes are shown for the delay parameter τ2=0, with trajectories plotted for two distinct values of τ1 (0.8 and 1), each represented by a different color code.Sub-figure (b) depicts the Hopf bifurcation with respect to the delay parameter τ1 while keeping τ2=0.

In Fig 5(a), the time delay $[eqn]$ is fixed, and the trajectories are plotted for three distinct values of $[eqn]$ , color coded as follows: red for $[eqn]$ , and blue for $[eqn]$ . It is evident that as the value of $[eqn]$ increases, the oscillations in the system populations become more pronounced initially but gradually diminish over time, indicating stabilization of the system. This demonstrates that larger values of $[eqn]$ contribute to enhanced system stability. In Fig 5(b), we present the bifurcation diagram by varying the bifurcation parameter $[eqn]$ from 0.8 to 1.2 while keeping $[eqn]$ . The diagram reveals the occurrence of a Hopf bifurcation at a critical value $[eqn]$ . For $[eqn]$ , all system populations except the mesenchymal stem cell (MSC) exhibit stable behavior, whereas for $[eqn]$ , the stability is lost. This indicates that increasing the time delay in the interaction term of the dendritic cell growth equation (stimulated by T cell-derived cytokines) infiltrating the keratinocyte population induces more oscillatory behavior in the system, excluding the MSC population. This observation highlights the significant impact of the bifurcation parameter $[eqn]$ on system stability from a biological point of view.

In Sub-figure (a), the populations of T cells, dendritic cells, and keratinocytes are shown with a fixed delay of τ1=0.The trajectories are plotted for two distinct values of τ2 (0.85 and 1.2), each represented by a different color. Sub-figure (b) illustrates the Hopf bifurcation with respect to the delay parameter τ2, while the other delay parameter, τ1, is fixed at 0.

In Fig 6, the stability region of the endemic equilibrium point is depicted in the $[eqn]$ − $[eqn]$ plane. The color gradient represents the maximum real part of the eigenvalue, denoted as $[eqn]$ , where λ corresponds to the eigenvalue of the Jacobian matrix at the interior equilibrium point of the delayed system (1). The equilibrium point is unstable in the region where $[eqn]$ . It is observed that the critical values of $[eqn]$ and $[eqn]$ , denoted as $[eqn]$ and $[eqn]$ , increase in the figures. Consequently, the area of the stability region for the endemic equilibrium expands with increasing values of $[eqn]$ and $[eqn]$ . The values of all other parameters used in the analysis are provided in Table 1.

The stability region of the endemic equilibrium point is shown in the τ1−τ2 plane.The color coding denotes the max[Re(λ)], where λ is the eigenvalues corresponding to the delayed system.

5.3 Simulations of optimal control for the delayed model

In this subsection, we numerically solve the optimal control problem for the delayed system (21) and present the results obtained through simulation. The parameter values used for the simulations are taken from Table 1, while the initial values are specified in Eq (30). The weight constants in the objective functional are set as $[eqn]$ .

The optimal system (21) represents a two-point boundary value problem (BVP) with boundary conditions provided at two time points, t = 0 and t = tf. For the numerical simulations, we utilized the bvp4c solver in MATLAB, which is specifically designed for solving nonlinear two-point BVPs. For a detailed understanding of this numerical method, readers may refer to the work by Torres et al. [26].

To solve the system, we define the state variables T(t), D(t), K(t), S(t), the adjoint variables $[eqn]$ $[eqn]$ , and the control functions $[eqn]$ . A combination of forward-backward difference approximations is employed to solve the control-induced system along with the adjoint system. The solutions are then plotted, with all state variables shown over time in Fig 7. The red dotted trajectories represent the system population concentrations without control, while the blue trajectories correspond to the population behavior under optimal control.

Numerical solution of the optimal control problem with τ1=τ2=1.1 and remaining parameter values taken from Table 1.The red dotted trajectories represent the disease dynamics without optimal control, while the blue trajectories depict the population behavior under the influence of optimal control.

To determine the optimal control profiles, we apply the transversality conditions and enforce the boundedness of controls within the interval [0,1] as per Pontryagin’s principle. The control functions $[eqn]$ and $[eqn]$ for the delayed system (21) are computed and plotted over 30 days, as shown in Fig 8. From the results, we observe a similar pattern for both control profiles, though notable differences exist. The control $[eqn]$ (represented by trajectories with blue dots), associated with the TNF-α inhibitor drug, modulates the system within a short span of approximately 10 days. In contrast, the control $[eqn]$ (indicated by red dotted trajectories), corresponding to the IL-17 inhibitor, requires about 18 days to regulate the system effectively. This indicates that the combination of these biologic inhibitors optimally controls the delayed model system, mitigating the auto-immune response and keratinocyte over-expression. The inflammatory cytokine loop, regulated predominantly by T cell and dendritic cell mediated interactions, is thus brought under control, achieving an optimal therapeutic outcome.

Numerical approximation of Pontryagin’s extremals and the corresponding optimal control profiles v1(t) and v2*(t).*

In addition to these, Fig 9 presents the global sensitivity analysis of all parameters in our system using Latin Hypercube Sampling (LHS) with the Partial Rank Correlation Coefficients (PRCC) method, focusing on the keratinocyte population. We evaluate the PRCC and p-value for each parameter, where the p-value indicates the level of uncertainty associated with each parameter’s PRCC value. p-value <0.05 for each parameter signifies the statistical significance of the corresponding PRCC value. This scatter plot shows that the parameters $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ have a positive influence on disease progression, while the parameters δ, $[eqn]$ , $[eqn]$ , and $[eqn]$ have a negative influence. In this study, we consider the keratinocyte level as a marker for disease prediction. The overall sensitivity analysis helps in identifying the precise point of treatment for psoriatic disease.

PRCC scatter plots for all system parameters against the keratinocyte population, with a sample size of N = 1000.The PRCC values are calculated at 100 days and all p values are taken as up to four decimal places. The x-axis shows the residuals obtained from regressing the rank-transformed values of the parameter of interest on the rank-transformed values of all other parameters. The y-axis displays the residuals from regressing the rank-transformed keratinocyte values on the rank-transformed values of all other parameters, excluding the one being analyzed.

6 Discussion and conclusion

Autoimmune diseases like psoriasis involve complex immune mechanisms, where time delays play a crucial role in the onset and progression of symptoms. These delays arise from processes such as immune activation, cytokine production, and cellular recruitment. Specifically, the formation of synergistic loops by pro-inflammatory cytokines and the subsequent regulatory response by anti-inflammatory cytokines are inherently non-instantaneous. Discrete intracellular delays represent measurable periods between the triggering event and disease manifestation, varying based on patient characteristics, environmental factors, and the specific autoimmune condition. For psoriatic lesions, while symptoms may appear soon after a trigger, the full extent of the disease evolves over a longer time-frame.

In this study, we incorporated two intracellular immune response delays, $[eqn]$ and $[eqn]$ , into our mathematical model to represent the temporal dynamics of immune interactions in psoriasis. The first delay, $[eqn]$ , reflects the time required for the initial immune response, such as T cell activation and proliferation triggered by dendritic cells. The second delay, $[eqn]$ , accounts for subsequent processes, including cytokine production by activated T cells and the recruitment of additional immune cells to the site of inflammation.

The bifurcation analysis revealed the critical role of these delays in the system’s stability. When $[eqn]$ is fixed and $[eqn]$ is varied, we identified a critical threshold $[eqn]$ at which the system transitions from stability to instability. This suggests that the second delay, $[eqn]$ , significantly influences the overall immune response. Similarly, varying $[eqn]$ while keeping $[eqn]$ fixed revealed a critical value $[eqn]$ , beyond which the system loses stability. These findings underscore the importance of both delays in regulating the immune response and maintaining equilibrium in the system.

The incorporation of two distinct delays provides a realistic representation of the temporal dynamics of psoriasis-related immune responses. Our results highlight that both $[eqn]$ and $[eqn]$ critically affect the stability of the endemic equilibrium point E . Beyond certain threshold values, these delays can induce oscillatory behavior in the system, highlighting the destabilizing effects of prolonged immune response delays. This phenomenon has been identified through the stability region plot presented in the numerical simulation section.

Furthermore, we implemented optimal controls on the delayed model by incorporating the effects of two biologic inhibitors targeting the TNF-α and IL-17 pathways. These effects are incorporated into the model with associated cost functions to reduce their high therapeutic expense. Using Pontryagin’s Minimum Principle, we analytically and numerically derived the optimality conditions for the delayed system. The results indicate that the optimal control profiles for the biologics effectively regulate the dynamics of the system, except for the last population, by reducing oscillations and stabilizing the immune response. Notably, the TNF-α inhibitor demonstrated faster modulation, achieving control within approximately 10 days, while the IL-17 inhibitor required a longer duration of 18–20 days to achieve effective regulation.

The results of this study, particularly the analysis of delay effects, have important implications for clinical strategies in managing psoriasis. Our findings indicate that therapeutic interventions targeting the timing of immune responses can play a key role in disease control. Specifically, addressing cytokine signaling delays and stabilizing immune responses may help to reduce the chronic inflammation typically observed in psoriasis. The use of delay differential equations in the control-induced model enables clinicians and researchers to better understand and predict the timing of keratinocyte proliferation and immune cell interactions. In addition, the incorporation of optimal control provides a cost-effective approach to biologic therapy, which may help to lower the economic burden of long-term treatment. This study highlights the importance of including temporal dynamics in therapeutic models, thereby offering new insights for the design of advanced treatment strategies for psoriasis. In practice, this could support the development of optimized treatment protocols—such as adjusting the timing of immunosuppressive drug delivery—to minimize flare-ups and improve patient outcomes. Finally, the insights into delay-induced instability also support personalized medicine approaches, where treatment schedules can be adapted to the specific response profile of individual patients.