Optimal Treatment Strategies for Control Model of Psoriasis

Downloaded 07/08/20 to 54.163.42.124. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Optimal Treatment Strategies for Control Model of Psoriasis Ellina Grigorieva Evgenii Khailov Abstract Psoriasis is an autoimmune disorder with the symptom of chronic inflammation of the skin. Psoriasis makes skin cells build up quickly on the skin surface which causes the skin to appear red, dry, and with scaly patches. Adequate treatment for psoriasis is very challenging and a total cure of the disease still does not exit. Mathematical models have long been effective means to predict cellular behaviors of skin regulation under normal or pathological circumstances. In this paper, a proposed control model of psoriasis is described on a given time interval by a nonlinear control system of three differential equations involving the concentrations of dendritic cells (tissues macrophages), T-lymphocytes, and keratinocytes with medication intake as a control function. An optimal control problem of minimizing the release of keratinocytes at the end of the time interval is stated and studied using the Pontryagin maximum principle. Different types of optimal control dependent on model parameters are obtained analytically. Possible applications to an optimal drug therapy are discussed. 1 Introduction One defining characteristic of the modern world is the increase of immune disorders and an absence of cure. There are some drugs that help fight the symptoms of the diseases but unfortunately, no presently known treatment defeats autoimmune illnesses such as penicillin helps to defeat some bacterial diseases. There are many autoimmune disorders, however, in this paper we will focus only on psoriasis. Inflammation is one of the weapons used by the immune system to fight an invader. For example, when you catch a virus or develop a bacterial infection, a type of immune cell called T-cells go into action. When T-cells recognize something as an invader (also called an antigen), T-cells begin an inflammatory attack against the invader. This attack is carried out by cytokines, which are proteins that help control the immune system’s inflammatory response. Cytokines trigger inflammation, causing the blood vessels to expand and send more immune cells to different parts of the body. In the case of psoriasis, this inflammation happens in the skin, leading to the red, itchy, and scaly patches known as plaques. The skin has two layers: one is the dermis, made Copyright © by SIAM Unauthorized reproduction of this article is prohibited 86 Paul Deignan of collagen-based matrix consisting of blood vessels, fibroblast cells, nerve endings, hair, and sweat glands; the other top layer is the epidermis, consisting of packed cells called keratinocytes. The generation of cytokines stimulates the proliferation of keratinocytes in the epidermal area of the skin and also the antigenic adhesion molecules in the dermal blood vessels. Genetically, psoriasis is viewed as a genetic heterogeneous disorder. The two causes of the disease development are the involvement of several genes and interactions with the environment [16]. Scientists continue to study the complex relationship between the immune system and psoriatic disease. Researchers are working to identify the antigens that trigger the autoimmune response in psoriasis, to better understand the role played by different kinds of immune cells in psoriatic disease, and develop new therapies that target cytokines and other parts of the immune system. The immune system is not only the key to understanding the causes of psoriatic disease – it may be the key to treating it too. In 1979, researchers discovered on accident that a drug called cyclosporine that suppresses the immune system also clears psoriasis ([2, 12]). This offered one of the first clues that psoriasis was actually an autoimmune disease. Since then, many effective treatments directed toward the immune system have been developed for the symptoms of psoriasis. In order to understand why certain drugs are beneficial in treating psoriasis, mathematical models can be very useful. Epidermal dynamics are usually modeled in two ways. The first, the deterministic models define analytical solutions to stationary cell populations, such as the model by [22], which investigate cell transition by a constitutive equation. The second, which includes agent-based models that involve keratinocytes under certain biological and physical rules. Both models simulate the multi-layers of the skin regulated by cell differentiation, proliferation, migration, and death in which certain biochemical factors change the behaviors of the epidermal structure [25]. There has been research performed by mathematicians jointly with biologists to investigate the psoriasis pathogenesis dynamics and treatments effects. Some models of psoriasis are given in [11, 14, 16, 17, 18, 20, 22, 24]. However, no optimal control problem has been stated and solved for Downloaded 07/08/20 to 54.163.42.124. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php those models. To our knowledge, the only optimal control problems for psoriasis were solved in the works of Roy ([1, 3, 19]) where the authors numerically found an optimal treatment that would minimize the cumulative concentration of keratinocytes and the cost of treatment. While the results of these works are very interesting, the presence of a square of the control under the integral of the objective function makes the optimal control problem numerically simplistic but prohibits the attainment of the optimal analytical solution. Before solving an optimal control problem numerically, it would be beneficial to investigate the problem analytically in order to reveal some features and properties of its optimal solutions. Therefore, we will consider the control model proposed by Roy [19] with a different objective function. Here, M is a positive constant depending on the parameters σ, ω, γ1 , γ2 , T of the system (2.1) and its initial conditions l0 , m0 , k0 . The boundedness, positiveness and continuation of solutions for system (2.1) are established by the following lemma. Lemma 2.1. Let the inclusion (l0 , m0 , k0 ) ∈ Θ hold. Then for any admissible control u(t) the corresponding solutions l(t), m(t), k(t) for system (2.1) are defined on the entire interval [0, T ] and satisfy the inclusion : (l(t), m(t), k(t)) ∈ Θ, t ∈ (0, T ]. (2.3) Relationship (2.3) implies that the region Θ is a positive invariant set for system (2.1). We note that the local and global stability analysis 2 Control Model of equilibria for system (2.1) in cases γ1 = γ2 = 0 and On a given time interval [0, T ] we consider a nonlinear γ1 , γ2 ̸= 0, u(t) = 1 are presented in [16, 19]. control model of psoriasis skin disease defined by the Using results from [1, 16, 19, 20], parameters γ1 , γ2 , following system: λ, µ, ν of system (2.1) satisfy the following constraints:  ˙ = σ − δl(t)m(t) − γ1 u(t)l(t)k(t) − µl(t), l(t)    γ1 ̸= γ2 , λ > µ, λ > ν. (2.4) ṁ(t) = ω − βl(t)m(t) − νm(t), (2.1)  k̇(t) = ηl(t)m(t) + γ2 u(t)l(t)k(t) − λk(t),   3 Optimal Control Problem l(0) = l0 , m(0) = m0 , k(0) = k0 ; l0 , m0 , k0 > 0. For system (2.1) on the set of admissible controls Here m(t), l(t), k(t) are the concentrations of dendritic Ω(T ) we consider the problem of minimization of the cells (tissue macrophages), T-lymphocytes, and ker- concentration of keratinocytes at the terminal time T : atinocytes at a specific time t, respectively. The param{ } eters are: µ the per capita removal of T-lymphocytes, ν min J(u) = k(T ) . (3.5) the per capita removal of dendritic cells, λ is the decay u(·)∈Ω(T ) rate of keratinocytes; γ1 is the rate of activation of keratinocytes by T-lymphocytes, γ2 is the rate of growth We note that the optimal control problem (2.1),(3.5) difof keratinocytes due to T-lymphocytes; δ is the rate of fers from problems that are typically considered in the activation of T-lymphocytes by dendritic cells, β is the literature on the control of psoriasis models ([1, 3, 19]) in rate of activation of dendritic cells by T-lymphocytes, that the functional of (3.5) does not include an integral η = β + δ is the proportion at which commonly stim- of the square of the control u(t), which is responsible for ulated T-lymphocytes and dendritic cells join the ex- the total cost of the drug dosage. In psoriasis therapy, pansion of epidermal keratinocytes. Further, σ, ω are in most cases, either a skin cream or an oral medicaconstant rates of influx of T-lymphocytes and dendritic tion, such as cyclosporine, are used. Both prescribed cells, respectively. Moreover, u(t) is the control, which medications have regular daily dosage and are not as represents medication intake in order to restrict the in- harmful for patients as the drugs used in chemotherapy teraction between T-lymphocytes and epidermal ker- for cancer treatment ([23]). Therefore the total cost of atinocytes. It satisfies the inequalities: psoriasis treatment in the meaning of “harm” to a patient and that usually mathematically is described by 0 < umin ≤ u(t) ≤ 1. (2.2) an integral of the square of the control, can be ignored. For system (2.1) the set of all admissible controls Moreover, using the terminal functional in (3.5) instead Ω(T ) is formed by all possible Lebesgue measurable of the corresponding integral functional (see [1, 19]), functions u(t), which for almost all t ∈ [0, T ] satisfy simplifies the subsequent analytical arguments. The existence in problem (2.1),(3.5) of the opticonstraints (2.2). mal control u∗ (t) and the corresponding optimal soluNow, we introduce a region { } tions l∗ (t), m∗ (t), k∗ (t) for system (2.1), follows from Θ = (l, m, k) : l > 0, m > 0, k > 0, l + m + k < M . Lemma 2.1 and Theorem 4 ([13], Chapter 4). Copyright © by SIAM Unauthorized reproduction of this article is prohibited 87 Then, using equations and initial conditions of system (4.6) we obtain for the functions L(t), G(t), In order to analyze problem (2.1),(3.5) we apply the ∗ (t) the Cauchy problem: ψ 1 Pontryagin maximum principle ([15]). We define the  Hamiltonian L̇(t) = a(t)L(t) + γ1 G(t), t ∈ [0, T ],     Ġ(t) = −b(t)L(t)  H(l, m, k, u, ψ1 , ψ2 , ψ3 ) = (σ − δlm − γ1 ulk − µl)ψ1  { − c(t)G(t)  −1 +m (t) α(λ − ν)m2∗ (t) +(ω − βlm − νm)ψ2 + (ηlm + γ2 ulk − λk)ψ3 , ∗ (4.8) +λ(λ − µ)m∗ (t) − ω(λ − µ)} ψ1∗ (t),     where ψ1 , ψ2 , ψ3 are adjoint variables. By the Pon ψ̇1∗ (t) = −e(t)L(t) − G(t) + λψ1∗ (t),   tryagin maximum principle, for the control u∗ (t) and L(T ) = −γ2 , G(T ) = 0, ψ1∗ (T ) = 0. corresponding solutions l∗ (t), m∗ (t), k∗ (t) for sysNext, we will study properties of the switching tem (2.1) there necessary exists the vector-function function L(t), which together with corresponding to it ∗ ∗ ∗ ψ∗ (t) = (ψ1 (t), ψ2 (t), ψ3 (t)) such that: functions G(t), ψ1∗ (t) satisfies system (4.8). (i) ψ∗ (t) is the nontrivial solution of the adjoint system: By the first initial condition of system (4.8) and the  ∗ ∗ ∗ continuity of the switching function L(t), the following (t)) (t) − γ ψ (t) = −u (t)k (t)(γ ψ ψ̇  1 ∗ ∗ 2 1 3 1   ∗ ∗ ∗ lemma can be stated.  −m (t)(ηψ (t) − βψ (t) − δψ (t))  ∗ 3 2 1   ∗  +µψ (t),  1 Lemma 4.1. There exists such a value t0 ∈ [0, T ) that   ∗ ψ̇2 (t) = −l∗ (t)(ηψ3∗ (t) − βψ2∗ (t) − δψ1∗ (t)) the inequality L(t) < 0 holds for all t ∈ (t0 , T ]. (4.6) +νψ2∗ (t),    Corollary 4.1. Formula (4.7) and Lemma 4.1 yield  ψ̇ ∗ (t) = −u (t)l (t)(γ ψ ∗ (t) − γ ψ ∗ (t))  ∗ ∗ 2 3 1 1  3  the following relationship for the optimal control u∗ (t) :  ∗  (t), +λψ  3  ∗ ψ1 (T ) = 0, ψ2∗ (T ) = 0, ψ3∗ (T ) = −1; u∗ (t) = umin , t ∈ (t0 , T ]. (4.9) Downloaded 07/08/20 to 54.163.42.124. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 4 Pontryagin Maximum Principle Now, we will focus on the case when parameters of system (2.1) satisfy the inequalities: (ii) the control u∗ (t) maximizes the Hamiltonian H(l∗ (t), m∗ (t), k∗ (t), u, ψ1∗ (t), ψ2∗ (t), ψ3∗ (t)) α < 0, λ2 + 4αω(λ − ν)(λ − µ)−1 < 0. with respect to u ∈ [umin , 1] for almost all t ∈ [0, T ], and therefore the following relationship holds:  , if L(t) > 0,  1 ∀u ∈ [umin , 1] , if L(t) = 0, u∗ (t) = (4.7)  umin , if L(t) < 0, where, by Lemma 2.1, the function L(t) = γ2 ψ3∗ (t) − γ1 ψ1∗ (t) (4.10) Let us consider system (4.8) and its second differential equation. We evaluate the discriminant, D of the quadratic expression inside the braces: ( ) D = (λ − µ)2 λ2 + 4αω(λ − ν)(λ − µ)−1 . By inequalities (2.4) and (4.10), it is easy to see that this discriminant is negative. Therefore, by Lemma 2.1, the function 2 d(t) = m−1 ∗ (t)(α(λ−ν)m∗ (t)+λ(λ−µ)m∗ (t)−ω(λ−µ)) is the switching function, which defines the type of the takes only negative values for all t ∈ [0, T ]. optimal control u∗ (t) according to formula (4.7); By introducing one more auxiliary function P (t) = (iii) the Hamiltonian ψ1∗ (t), we finally rewrite system (4.8) as  H(l∗ (t), m∗ (t), k∗ (t), u∗ (t), ψ1∗ (t), ψ2∗ (t), ψ3∗ (t)) L̇(t) = a(t)L(t) + γ1 G(t), t ∈ [0, T ],    Ġ(t) = −b(t)L(t) − c(t)G(t) + d(t)P (t), is constant on the given interval [0, T ]. (4.11)  Ṗ (t) = −e(t)L(t) − G(t) + λP (t),  Now, let us define the constant α = γ2−1 (ηγ1 − δγ2 )  L(T ) = −γ2 , G(T ) = 0, P (T ) = 0. and following functions: Analyzing this system, the validity of the following a(t) = u∗ (t)(γ1(k∗ (t) − γ2 l∗ (t)) + ηγ1 γ2−1 m∗)(t) + λ, lemma can be shown. 2 c(t) = m−1 ∗ (t) αm∗ (t) + (λ − µ)m∗ (t) − ω , Lemma 4.2. The switching function L(t) cannot be e(t) = u∗ (t)k∗ (t) + ηγ2−1 m∗ (t), −1 zero on any subinterval of the interval [0, T ]. b(t) = (αm∗ (t) + (λ − µ))e(t) − βηγ2 m∗ (t)l∗ (t). Corollary 4.2. Formula (4.7) and Lemma 4.2 imply that the optimal control u∗ (t) is bang-bang function taking values {umin ; 1} for almost all t ∈ [0, T ]. Also, we introduce the auxiliary function: G(t) = −m∗ (t)(βψ2∗ (t) − αψ1∗ (t)) + (λ − µ)ψ1∗ (t). Copyright © by SIAM Unauthorized reproduction of this article is prohibited 88 Downloaded 07/08/20 to 54.163.42.124. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 5 Estimation of the Number of Zeros Let us estimate the number of zeros of the switching function L(t). For this, we use system (4.11), which is a linear non-autonomous homogeneous system of differential equations. To analyze this system, we apply the ideas presented in [6, 7]. Let us rewrite system (4.11) in a matrix form: (5.12) Now, let us add to system (5.14) the initial conditions:    L(t) 0 d(t)  , r(t) =  G(t)  . P (t) λ h21 (0) = h021 , h31 (0) = h031 , h32 (0) = h032 , (5.16) ṙ(t) = A(t)r(t), where  a(t) A(t) =  −b(t) −e(t) γ1 −c(t) −1 Now, we reduce the matrix triangular form:  ∗ Π(t) =  0 0 of this system to the upper γ1 ∗ 0  0 d(t)  . ∗ To do this, let us make the following nonsingular transformation of the vector r(t): ρ(t) = H(t)r(t), where  1 H(t) =  h21 (t) h31 (t) In turn, transformed system (5.13) has the from:   e ė = (a(t) − γ1 h21 (t)) L(t) e + γ1 G(t), L(t)     ė e G(t) = (γ1 h21 (t) − c(t) − d(t)h32 (t)) G(t) (5.15) e  +d(t)P (t),     ė P (t) = (λ + d(t)h32 (t)) Pe (t).   0 0  1 0  , ρ(t) =  h32 (t) 1  e L(t)  e G(t) . Pe(t) Here, h21 (t), h31 (t), h32 (t) are the functions, which we will choose later. It is easy to see that under such transformation the function L(t) does not change, i.e. e = L(t). The other functions can be written as L(t) e = h21 (t)L(t) + G(t), G(t) e P (t) = h31 (t)L(t) + h32 (t)G(t) + P (t). In terms of the new variables given by the components of vector ρ(t), system (5.12) has the following form: ) ( (5.13) ρ̇(t) = Ḣ(t) + H(t)A(t) H −1 (t)ρ(t). where h021 , h031 , h032 are some constants. Hence, we obtain the Cauchy problem, and as a corollary of the corresponding Existence and Uniqueness Theorem ([9]), its solutions will be defined, generally speaking, locally, in some neighborhood of value t = 0. As it follows from [6, 7], it is important to find restrictions for the coefficients of the equations of system (5.14) at which there exist solutions of Cauchy problem (5.14),(5.16) defined on the entire interval [0, T ]. In order to simplify the analysis of the quadratic system of the third order (5.14),(5.16) let us reduce its order by one, i.e. we will split this system into two more simple quadratic subsystems of the second order. Equations for the functions h31 (t), h32 (t) together with the corresponding initial conditions from (5.16) form the first quadratic subsystem:  ḣ31 (t) = e(t) − (a(t) − λ)h31 (t) + b(t)h32 (t)     +d(t)h31 (t)h32 (t),  ḣ32 (t) = 1 − γ1 h31 (t) + (c(t) + λ)h32 (t) (5.17)    +d(t)h232 (t),   h31 (0) = h031 , h32 (0) = h032 . Then, in the differential equation for the function h21 (t) of system (5.14) we introduce a new variable h20 (t) using formula h20 (t) = h32 (t)h21 (t) − h31 (t). Substituting it into the first equation of system (5.14) and then obtaining the differential equation for the function h20 (t) and also its initial condition h20 (0), we find the second quadratic subsystem:    ḣ20 (t) = −e(t) − (a(t) − λ)h20 (t) + h21 (t)   +γ1 h20 (t)h21 (t),  (5.18) ḣ21 (t) = b(t) − (a(t) + c(t))h21 (t)  2   −d(t)h (t) + γ h (t), 20 1 21   h20 (0) = h020 = h032 h021 − h031 , h21 (0) = h021 . We will choose functions h21 (t), h31 (t), h32 (t) in a such way that relationship Π(t)H(t) = Ḣ(t) + H(t)A(t) is true. Then, we find non-autonomous nonhomogeneous quadratic system of differential equations to which these functions satisfy:  ḣ21 (t) = b(t) − (a(t) + c(t))h21 (t) + d(t)h31 (t)     −d(t)h32 (t)h21 (t) + γ1 h221 (t),    Note that the quadratic system of the third orḣ31 (t) = e(t) − (a(t) − λ)h31 (t) + b(t)h32 (t) (5.14) der (5.14),(5.16) is equivalent to the two quadratic sub+d(t)h (t)h (t),  31 32    systems of the second order (5.17) and (5.18). Hence,  ḣ (t) = 1 − γ1 h31 (t) + (c(t) + λ)h32 (t)   32 further we will investigate these quadratic subsystems +d(t)h232 (t). Copyright © by SIAM Unauthorized reproduction of this article is prohibited 89 Downloaded 07/08/20 to 54.163.42.124. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php in order to find the restrictions for the coefficients of the equations of these subsystems at which there simultaneously exist their solutions defined on the entire interval [0, T ]. By [9, 13], we will consider that the solutions h31 (t), h32 (t) for system (5.17) and solutions h20 (t), h21 (t) for system (5.18) are defined on the corresponding maximum possible intervals ∆3 = (t30 , t31 ) and ∆2 = (t20 , t21 ) of the existence of such solutions. Without loss of generality, we assume that t31 ≤ T and t21 ≤ T . Now, at systems (5.17) and (5.18) let us make the change of variables: Next, with nonnegative initial conditions: e h031 ≥ 0, e h032 ≥ 0, e h020 ≥ 0, e h021 ≥ 0 let us apply to systems (5.19) and (5.20) the necessary and sufficient condition of non-negativity of the solutions e h31 (t), e h32 (t) for all t ∈ [0, t31 ) and e h20 (t), e h21 (t) 2 for all t ∈ [0, t1 ) ([10]), which implies the execution of the so-called quasi-positivity condition: Φ31 (t, 0, e h32 ) ≥ 0 for all e h32 ≥ 0, Φ32 (t, e h31 , 0) ≥ 0 for all e h31 ≥ 0 for system (5.19) and e h31 = −h31 − K31 , e h32 = h32 − K32 , e h20 = −h20 − K20 , e h21 = −h21 − K21 , Φ20 (t, 0, e h21 ) ≥ 0 for all e h21 ≥ 0, e e where K20 , K21 , K31 , K32 are the constants, which we Φ21 (t, h20 , 0) ≥ 0 for all h20 ≥ 0 choose later. In the new variables e h31 , e h32 and e h20 , e h21 for system (5.20). systems (5.17) and (5.18) can be written as As a result, we find the system of inequalities:  ė e e  h (t) = Φ (t, h (t), h (t))   31 31 31 32   −γ1−1 < K31 < 0, K32 > 0,   e e   = d(t) h (t) h (t) 31 32   ( )   K d(t) − b(t) > 0,      31  h31 (t) + K32 d(t) − (a(t) − λ) e 2  K32 d(t) + K32 (c(t) + λ)   (5.21)  e  + (K d(t) − b(t)) h (t) +(K31 γ1 + 1) ≥ 0,   31 32 (       K (K d(t) − b(t))   + K32 (K31 d(t) − b(t))    32 31  )  −(K31 (a(t) − λ) + e(t)) ≥ 0,    −(K31 (a(t) − λ) + e(t)) ,   (5.19) related to system (5.19), and the system of inequalities: ė h32 (t) = Φ32 (t, e h31 (t), e h32 (t))      h31 (t) ) =(d(t)e h232 (t) + γ1 e 0 < K20 < γ1−1 , K21 > 0,         e K d(t) + b(t) < 0,   + 2K32 d(t) + (c(t) + λ) h32 (t)     20 ( 2  K21 γ1 + K21 (a(t) + c(t))   2 (5.22)  + K d(t) + K (c(t) + λ) 32  32 +(K20 d(t) + b(t)) ≤ 0,   )       K (K γ − 1)   +(K31 γ1 + 1) ,     21 20 1  +(K20 (a(t) − λ) − e(t)) ≤ 0,  e 0 0 e  h31 (0) = h31 = −h31 − K31 ,    e 0 0 h32 (0) = e h32 = h32 − K32 . corresponding to system (5.20). These systems give the  required restrictions on the functions a(t), b(t), c(t), ė  h20 (t) = Φ20 (t, e h20 (t), e h21 (t))  d(t), e(t) and the parameters γ1 , λ.     =(−γ1 e h20 (t)e h21 (t) In subsequent arguments, we will assume that the  )    e following condition holds.  h (t) − K γ + (a(t) − λ) 20 21 1     e  −(K  Assumption 1. At each value of t ∈ [0, T ] the sets of ( 20 γ1 − 1)h21 (t)     − K (K γ − 1) solutions for systems of inequalities (5.21) and (5.22)  21 20 1  )  ∗ ∗  are not empty. There exist the constants K31 , K32   +(K20 (a(t) − λ) − e(t)) ,  ∗ ∗  and K20 , K21 belonging to the corresponding sets of (5.20) solutions of these systems at once for all t ∈ [0, T ]. ė e e h (t) = Φ (t, h (t), h (t)) 21 21 20 21     =(−γ1 e h221 (t) − d(t)e h20 (t) )   A cursory analysis of systems (5.21) and (5.22) leads    h21 (t) − 2K21 γ1 + (a(t) + c(t)) e  to the conclusion about the validity of this assumption.  (    2 Hence, the subject of our further investigations is the  − K21 γ1 + K21 (a(t) + c(t))   )  search of relationships as between parameters of sys   +(K20 d(t) + b(t)) ,  tem (2.1) as between its initial conditions under which    0 0 e  e − K , = −h h (0) = h an accurate analysis of these systems confirms the va 20 20 20 20   e lidity of Assumption 1. h21 (0) = e h021 = −h021 − K21 . Copyright © by SIAM Unauthorized reproduction of this article is prohibited 90 Downloaded 07/08/20 to 54.163.42.124. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Now, we will consider that all previous arguments values of the functions a(t), b(t), c(t), d(t), e(t) on the interval [0, T ]. Here the sign ⊤ means transpose. These were made for the constants: inequalities are obtained by the constraints: ∗ ∗ ∗ ∗ K31 = K31 , K32 = K32 , K20 = K20 , K21 = K21 . e h231 (t) + e h232 (t) ≥ 1, e h220 (t) + e h221 (t) ≥ 1, e e Next, for nonnegative solutions h31 (t), h32 (t), which, in turn, occur due to the relationships: e h20 (t), e h21 (t) let us check the condition of Theorem 1.6 lim ||(e h31 (t), e h32 (t))⊤ || = +∞, from [10] about the continuability of these solutions on t→t31 −0 the given interval [0, T ]. To do this, we introduce the following Lyapunov functions: lim ||(e h20 (t), e h21 (t))⊤ || = +∞, t→t21 −0 ( ) h232 , h231 + e V3 (e h31 , e h32 ) = 0.5 e following from Lemma ([4], §14, Chapter 4). We note that in (5.23) and (5.24) the non-positivity of the first ( ) h221 , h220 + e V2 (e h20 , e h21 ) = 0.5 e terms in the scalar products plays a significant role for their evaluation. The inequalities obtained in (5.23) and then evaluate the appropriate scalar products of the and (5.24) imply the required continuability of the gradients of these Lyapunov functions and the corresolutions e h31 (t), e h32 (t) for system (5.19) and solutions sponding right-hand sides of systems (5.19) and (5.20): e e h (t), h (t) for system (5.20) on the given interval 20 21 ( [0, T ]. Returning back, to the functions h31 (t), h32 (t), e e ∇V3 (h31 (t), h32 (t)), ) h20 (t), h21 (t) satisfying systems (5.17) and (5.18), which (Φ31 (t, e h31 (t), e h32 (t)), Φ32 (t, e h31 (t), e h32 (t)))⊤ are also obviously defined on the entire interval [0, T ], ) ( we conclude that system (5.15) is defined on this entire 2 2 e e e = d(t)h32 (t) h31 (t) + h32 (t) ( ) interval as well. ∗ h2 (t) + K32 d(t) − (a(t) − λ) e ( ) 31 Remark 1. Conducted arguments show that the sets: ∗ + (K31 d(t) − b(t)) + γ1 e h31 (t)e h32 (t) { } ( ) ∗ ∗ , h ≥ K (h , h ) : h ≤ −K , 32 31 32 31 32 31 ∗ (5.23) h232 (t) + 2K32 d(t) + (c(t) + λ) e ( { } ∗ ∗ ∗ ∗ d(t) − b(t)) (K31 + K32 (h20 , h21 ) : h20 ≤ −K20 , h21 ≤ −K21 ) ∗ h31 (t) −(K31 (a(t) − λ) + e(t)) e are invariant sets for systems (5.17) and (5.18), respec( tively. ∗ 2 ∗ + K32 d(t) + K32 (c(t) + λ) ) By the first initial condition of system (4.11) and ∗ h32 (t) γ1 + 1) e +(K31 the generalized Rolle’s Theorem ([5]) applied one by ≤ Q3 V3 (e h31 (t), e h32 (t)), one to the each equation of system (5.15), the following ( lemma can be stated. ∇V2 (e h20 (t), e h21 (t)), ) Lemma 5.1. The switching function L(t) has on the (Φ20 (t, e h20 (t), e h21 (t)), Φ21 (t, e h20 (t), e h21 (t)))⊤ interval [0, T ) at most two distinct zeros. ) ( 2 2 e e e = −γ1 h21 (t) h20 (t) + h21 (t) Corollary 5.1. Formula (4.7), relationship (4.9) and ( ) 2 ∗ Lemma 5.1 imply that the optimal control u∗ (t) depende − K21 γ1 + (a(t) − λ) h20 (t) ) ( ing on the relationships as between parameters of sys∗ γ1 − 1) e h20 (t)e h21 (t) − d(t) + (K20 tem (2.1) as between its initial conditions has one of ( ) 2 ∗ the following three types : e (5.24) − 2K21 γ1 + (a(t) + c(t)) h21 (t) ( ∗ ∗ u0∗ (t) = umin , t ∈ [0, T ], (5.25) − K21 (K20 γ1 − 1) ) { ∗ umax , if 0 ≤ t ≤ θ∗ , h20 (t) (a(t) − λ) − e(t)) e +(K20 u1∗ (t) = (5.26) ( u if θ∗ < t ≤ T, min , ∗ 2 ∗ − K21 γ1 + K21 (a(t) + c(t))  )  umin , if 0 ≤ t ≤ θ1∗ , ∗ e 2 +(K20 d(t) + b(t)) h21 (t) umax , if θ1∗ < t ≤ θ2∗ , u∗ (t) = (5.27)  e e umin , if θ2∗ < t ≤ T, ≤ Q2 V2 (h20 (t), h21 (t)), ∗ ∗ where Q2 , Q3 are the constants depending on the where θ∗ ∈ [0, T ) and θ1 , θ1 ∈ [0, T ) are the moments of 1 2 parameters γ1 , λ and the estimates of maximal absolute switching of the corresponding controls u∗ (t), u∗ (t). Copyright © by SIAM Unauthorized reproduction of this article is prohibited 91 Downloaded 07/08/20 to 54.163.42.124. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 6 Solution to the Optimal Control Problem Now we outline a method of solving optimal control problem (2.1),(3.5). Let us define the set: { } Λθ = θ = (θ1 , θ2 ) : 0 ≤ θ1 ≤ θ2 ≤ T . For arbitrary point θ ∈ Λθ as   umin , umax , uθ (t) =  umin , we define the control uθ (t) if 0 ≤ t ≤ θ1 , if θ1 < t ≤ θ2 , if θ2 < t ≤ T. It is easy to see that the control uθ (t), defined in this way, includes the types of the optimal control u∗ (t) given by formulas (5.25)–(5.27). Substituting the control uθ (t) into system (2.1), and then integrating this system over interval [0, T ], we find the functions lθ (t), mθ (t), kθ (t), which correspond to this control. In particular, kθ (T ) is the value of the functional from (3.5) corresponding to this control, or, in other words, to the point θ ∈ Λθ . The above discussed method defines a function of two variables: F (θ) = J(uθ ), θ ∈ Λθ . Thus, problem (2.1),(3.5) is now reduced to a problem of constrained minimization: { } min F (θ) . (6.28) θ∈Λθ Methods for numerical solution of such problem are well-known ([26]). Problem (6.28) is considerably simpler than optimal control problem (2.1),(3.5) and can be solved numerically by using MAPLE. The results of the corresponding numerical calculations and their analysis will be presented in our report. The choice of the terminal time T depends on a patient and the type of his or her treatment (skin cream or an oral drug). For example, cyclosporine can provide rapid relief from symptoms. One may see some improvement in symptoms after two weeks of treatment, particularly with stronger doses. However, it may take from three to four months to reach optimal results. By [2, 3, 19, 20], this time typically varies from 30 to 120 days. The values of other parameters of the system (2.1) and its initial conditions are presented in [19]. 7 Discussion Let us compare formulas (5.25)–(5.27) for optimal control u∗ (t) with the results presented in [19]. To do this, we note that the convenient for analytical analysis, auxiliary control u(t) in system (2.1) is related to the corresponding physical control v(t) in the same system Copyright © by SIAM Unauthorized reproduction of this article is prohibited 92 from [19] by the formula: u(t) = 1 − v(t). Therefore, where the auxiliary optimal control u∗ (t) has a maximum value of 1, the appropriate physical optimal control v∗ (t) takes a minimum value of 0, and vice versa. Thus, physical optimal control v∗ (t) either according to formula (5.25) takes the maximum value on the entire interval [0, T ] that corresponds to the psoriasis treatment with the greatest intensity throughout a given time interval (emergency drug therapy); or according to formula (5.26) has one switching from the minimum value to the maximum value that describes the situation when, first there is the period of the psoriasis treatment with the lower intensity (adjustment period), and then, at the time t = θ∗ the switching occurs to the period of the treatment with the greatest intensity (drug therapy with adjustment period). Finally, according to formula (5.27) the physical optimal control v∗ (t) has two switchings: first there is the switching from the maximum value to the minimum value, and then again to the maximum value. This corresponds to the situation when, first the psoriasis treatment is provided with the greatest intensity, then, at the time t = θ1∗ the switching occurs, and the treatment is performed with the lower intensity (rest period), then, at the time t = θ2∗ the treatment is returned to the greatest intensity (drug therapy with rest period). Thus, three types of psoriasis optimal strategies can be distinguished: the emergency drug therapy, the drug therapy with adjustment period, and the drug therapy with rest period. These therapies are consistent both with the actual data of clinical studies ([8, 12, 21]) and the results of the numerical calculations from [19]. 8 Conclusions On the given time interval we considered a control model of psoriasis skin disease described by a nonlinear system of three differential equations. Phase variables of this system denote the concentrations of dendritic cells, Tlymphocytes and keratinocytes. The model contained a bounded control representing medication intake. This control was introduced into the model in order to investigate optimal treatment strategies for minimizing the concentration of keratinocytes at the end of the predefined period of time. We used such values of the model parameters for which the corresponding optimal control was bang-bang. A new approach for estimating the number of zeros of the corresponding switching function, different from the one that was used in our previous papers, was applied. Having the maximal possible number of switchings, we finally reformulated the original optimal control problem as a problem of constrained minimization, which then may be solved numerically using standard software. Downloaded 07/08/20 to 54.163.42.124. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php References [1] X. Cao, A. Datta, F. Al Basir, and P.K. Roy, Fractional-order model of the disease psoriasis: a control based mathematical approach, J. Syst. Sci. Complex, 29 (2016), pp. 1565–1584. [2] M.D. Colombo, N. Cassano, G. Bellia, and G.A. Vena, Cyclosporine regimens in plaque psoriasis: an overview with special emphasis on dose, duration, and old and new treatment approaches, The Scientific World Journal, 2013 (2013), Article ID 805705, pp. 1–11. [3] A. Datta and P.K. Roy, T-cell proliferation on immunopathogenic mechanism of psorias: a control based theoretical approach, Control Cybern., 42 (2013), pp. 365–386. [4] B.P. Demidovich, Lectures on Stability Theory, Nauka, Moscow, 1967. [5] A.V. Dmitruk, A generalized estimate on the number of zeros for solutions of a class of linear differential equations, SIAM J. Control Optim., 30 (1992), pp. 1087– 1091. [6] E.V. Grigorieva, E.N. Khailov, and A. Korobeinikov, Parametrization of the attainable set for a nonlinear control model of a biochemical process, Math. Biosci. Eng., 10 (2013), pp. 1067–1094. [7] E.V. Grigorieva and E.N. Khailov, Optimal intervention strategies for a SEIR control model of Ebola epidemics, Mathematics, 3 (2015), pp. 961–983. [8] J.E. Gudjonsson, A. Johnston, H. Sigmundsdottir, and H. Valdimarsson, Immunopathogenic mechanisms in psoriasis, Clin. Exp. Immunol., 135 (2004), pp. 1–8. [9] P. Hartman, Ordinary Differential Equations, John Wiley & Sons, New York, 1964. [10] M.A. Krasnosel’skii, The Operator of Translation along the Trajectories of Differential Equations, AMS, Providence, RI, 1968. [11] J.R. Lancaster, Simulation on the diffusion and reaction of endogenously produced nitric oxide, Proceedings of the National Academy of Science USA, 91 (1994), pp. 8137–8141. [12] M. Lebwohl, Psoriasis, Lancet, 361 (2003), pp. 1197– 1204. [13] E.B. Lee and L. Marcus, Foundations of Optimal Control Theory, John Wiley & Sons, New York, 1967. [14] G. Niels and N. Karsten, Simulating psoriasis by altering transit amplifying cells, Oxford Journals, Science and Mathematics. Bioinformatics, 23 (2007), pp. 1309– 1312. [15] L.S. Pontryagin, V.G. Boltyanskii, R.V. Gamkrelidze, and E.F. Mishchenko, Mathematical Theory of Optimal Processes, John Wiley & Sons, New York, 1962. [16] P.K. Roy, J. Bradra, and B. Chattopadhyay, Mathematical modeling on immunopathogenesis in chronic plaque of psoriasis: a theoretical study, Lecture Notes in Engineering and Computer Science, 1 (2010), pp. 550–555. [17] P.K. Roy and A. Datta, Negative feedback control may regulate cytokines effect during growth of keratinocytes Copyright © by SIAM Unauthorized reproduction of this article is prohibited 93 [18] [19] [20] [21] [22] [23] [24] [25] [26] in the chronic plaque of psoriasis: a mathematical study, International Journal of Applied Mathematics, 25 (2012), pp. 233–254. P.K. Roy, A. Datta, and A.N. Chatterjee, Saturation effects on immunopathogenic mechanism of psoriasis: a theoretical approach, Acta Analysis Functionalis Applicata, 13 (2013), pp. 310–318. P.K. Roy and A. Datta, Impact of cytokine release in psoriasis: a control based mathematical approach, Journal of Nonlinear Evolution Equations and Applications, 3 (2013), pp. 23–42. P.K. Roy and A. Datta, Impact of perfect drug adherence on immunopathogenic mechanism for dynamical system of psoriasis, Biomath, 2 (2013), Article ID 121201, pp. 1–6. R. Sabat, S. Philipp, C. Höflich, S. Kreutzer, E. Wallace, K. Asadullah, H.D. Volk, W. Sterry, and K. Wolk, Immunopathogenesis of psoriasis, Exp. Dermatol., 16 (2007), pp. 779–798. N.J. Savill, R. Weller, and J.A. Sherratt, Mathematical modelling of nitric oxide regulation of Rete Peg formation in psoriasis, J. Theor. Biol., 214 (2002), pp. 1–16. H. Schättler and U. Ledzewicz, Optimal Control for Mathematical Models of Cancer Therapies: An Application of Geometric Methods, Springer, New YorkHeidelberg-Dordrecht-London, 2015. J.A. Sherratt, R. Weller, and N.J. Savill, Modelling blood flow regulation by nitric oxide in psoriatic plaques, B. Math. Biol., 64 (2002), pp. 623–641. M. Trehan and C.R. Taylor, Medium-dose 308-nm excimer laser for the treatment of psoriasis, J. Amer. Acad. Dermatol., 47 (2002), pp. 701–708. F.P. Vasil’ev, Optimization Methods, Factorial Press, Moscow, 2002.