Fuzzy Modeling and Control of HIV Infection

Abstract

The present study proposes a fuzzy mathematical model of HIV infection consisting of a linear fuzzy differential equations (FDEs) system describing the ambiguous immune cells level and the viral load which are due to the intrinsic fuzziness of the immune system’s strength in HIV-infected patients. The immune cells in question are considered CD4+ T-cells and cytotoxic T-lymphocytes (CTLs). The dynamic behavior of the immune cells level and the viral load within the three groups of patients with weak, moderate, and strong immune systems are analyzed and compared. Moreover, the approximate explicit solutions of the proposed model are derived using a fitting-based method. In particular, a fuzzy control function indicating the drug dosage is incorporated into the proposed model and a fuzzy optimal control problem (FOCP) minimizing both the viral load and the drug costs is constructed. An optimality condition is achieved as a fuzzy boundary value problem (FBVP). In addition, the optimal fuzzy control function is completely characterized and a numerical solution for the optimality system is computed.

1. Introduction

Usage of fuzzy differential equations is a natural way to model dynamical systems under uncertainty [1]. For example these equations are used to modeling the cell growth and dynamic of population [2], dry friction [3], tumor growth [4], and the phenomenon of nuclear disintegration [5] under uncertainty. In [6], transition from HIV to AIDS (the acquired immunodeficiency syndrome) is described through a mathematical model with fuzzy transference rate correlated with the viral load and CD4+ T-cells level by rule bases. Moreover, in [7], the authors have proposed a methodology combining a macroscopic HIV-positive population model, which is a differential equation system whose transference rate from asymptomatic to symptomatic population is found through a fuzzy rule-based system, with an individual microscopic model to study the evolution of positive HIV population for manifestation of AIDS. In [8], a fuzzy delay differential equation is proposed to model HIV infection, assuming that there exists delay between the infection of a CD4+ T-cell by the virus and the production of new virus particles. In this model, the delay and the clearance rate of HIV particles are fuzzy numbers where correlation between them is restated by rule bases. However, it should be noted that the whole parameters of a model such as the production and clearance rate of viruses and immune cells can be the source of uncertainty. In real world, there are various HIV-infected patients with different strengths of immune system causing uncertainty as to the immune cells level and the viral load during the different stages of the disease. A number of mathematical models have been formulated to describe various aspects of the interaction between HIV and the immune cells. The basic and simple model of HIV infection that contains three state variables: healthy CD4+ T-cells, infected CD4+ T-cells, and viruses, is presented by Perelson et al. [9], and more complicated models containing other parts of the immune system such as the cytotoxic T-lymphocyte and the macrophages are presented in [10] and references therein. None of these models can mirror the mentioned uncertainties proposing a mathematical model with fuzzy parameters which could reflect such ambiguities would be desirable.

One of the earliest suggestions to define the concept of differentiability for fuzzy mappings and, in consequence, to study fuzzy differential equations is the Hukuhara derivative [11]. Nevertheless, the solution of fuzzy differential equation interpreted by Hukuhara derivative became fuzzier as time goes by [12]. Hence, the fuzzy solution behaves quite differently from the crisp solution. In order to overcome this difficulty, a more general definition of derivative for fuzzy-number-valued functions, which is called the strongly generalized differentiability, has been introduced and many papers have been published in this field (see, e.g., [1, 2, 5, 13–17] and references therein). In [13] a generalization of the Hukuhara differentiability to the case of interval-valued functions is introduced and the local existence and uniqueness of solutions for the interval differential equations are obtained under this type of differentiability.

First-order linear fuzzy differential equations are one of the simplest fuzzy differential equations which may appear in many applications. However, the form of such an equation is very simple, it raises many problems since under different fuzzy differential equation concepts, the behavior of the solutions is different (depending on the interpretation used). This type of equations has been studied by many researchers. For example, the general form of the solutions for the first-order fuzzy differential equations with crisp coefficients under the generalized differentiability concept is presented in [14]. Moreover, an operator method is proposed for solving a class of first- and second-order linear fuzzy differential equations under the assumption of strongly generalized differentiability which is constructed based on their equivalent integral forms [15]. The existence and uniqueness of the solutions has been demonstrated for a first-order linear fuzzy differential equation with impulses subject to boundary value conditions, and the explicit solutions are obtained by calculating the solutions on each level set [16]. The generalized Euler approximation method is applied to solve numerically fuzzy differential equations under generalized differentiability [17]. A number of works in these fields have dealt with the linear fuzzy differential dynamical systems. For instance, in [18], a complex number representation of the α-level sets of the linear first-order fuzzy differential dynamical systems where the initial condition is described by a vector of fuzzy numbers is presented and the solutions are obtained under such representation. Using this approach, a method is proposed to find the solutions of a class of linear differential dynamical systems with fuzzy matrices [19] and the proposed method is extended to provide the solutions of linear matrix differential dynamical systems with fuzzy matrices [20].

Classical control system is described by a differential equation. However, uncertainty is inherit in most dynamic systems. The concept of fuzzy optimal control was presented by Komolov et al. [21] in 1979. Since then, many researchers have studied this type of problems. In [22], fuzzy differential equations are generalized to be fuzzy set control differential equations (FSCDEs) and the problem of stability and controllability of FSCDE are presented. Furthermore, some properties of the fuzzy solution for the linear FSCDE as well as the necessary and sufficient optimality conditions for a linear fuzzy time optimal control problem are obtained in [23].

In this paper, we model the uncertain behaviors of CD4+ T-cells and CTLs level and the viral load in different patients by a system of linear fuzzy differential equations and analyze the optimal control regarding minimizing both the viral load and drug costs.

Following a preliminary introduction, in Section 3, a system of linear differential equations with fuzzy parameters describing the ambiguous behaviors of CD4+ T- cells and CTLs level and the HIV viral load in patients with a weak, moderate and strong immune system is introduced. Moreover, a method for finding explicit solutions to the proposed model is introduced in this section. Some authors have used mathematical models for HIV infection in conjunction with control theory to achieve appropriate goals. Although the proposed model is simple, it can be used to investigate the effects of antiretroviral therapy in preventing the HIV progression. Section 4 is devoted to the latter topic. The last section deals with the conclusion.

2. Preliminaries

In this section, we give some definitions and introduce the necessary notations which will be used throughout the paper. See, for example, [12].

Definition 1. (fuzzy set, α-level set, and fuzzy number) A fuzzy set 𝑢 in 𝑅^𝑛 is defined as a set of all pairs (𝑥,𝜇_𝑢(𝑥))∈𝑅^𝑛×[0,1] for some function 𝜇_𝑢∶𝑅^𝑛→[0,1], which is called the membership function of 𝑢, and 𝜇_𝑢(𝑥) is interpreted as the membership grade of a element 𝑥 in the fuzzy set 𝑢. We define [𝑢]^𝛼={𝑥∈𝑅∣𝜇_𝑢(𝑥)≥𝛼} the α-level set of 𝑢, with 0<𝛼≤1. For 𝛼=0, the support of 𝑢 is defined as [𝑢]⁰=supp(𝑢)=cl{𝑥∈𝑅∣𝜇_𝑢(𝑥)>0}, where cl denotes the closure of a subset. A fuzzy set 𝑢 in 𝑅 is called a fuzzy number in 𝑅 if(i)𝜇_𝑢 is upper semicontinuous on 𝑅,(ii)𝑢 is a convex fuzzy set, that is, 𝜇_𝑢(𝜆𝑥+(1−𝜆)𝑦)≥min{𝜇_𝑢(𝑥),𝜇_𝑢(𝑦)}, for all 𝑥,𝑦∈𝑅, 𝜆∈[0,1],(iii)𝜇_𝑢 is normal, that is, there exists a unique 𝑥₀∈𝑅 such that 𝜇_𝑢(𝑥₀)=1,(iv)the support of 𝑢 is compact.

Example 2. The normal level of CD4+ T-cells in blood, that is: “close to 1000 cells/μL,” can be represented as a fuzzy set 𝑢 with a membership function defined as 𝜇_𝑢(𝑥)=exp(−𝛽(𝑥−1000)²) where  β is a positive real number. It is easy to see that [𝑢]^𝛼=[1000−−𝛽⁻¹ln𝛼,1000+−𝛽⁻¹ln𝛼]. The membership function 𝜇_𝑢(⋅) with 𝛽=0.0005 and [𝑢]^𝛼 are shown in Figure 1. Note that [𝑢]⁰=(−∞,+∞); hence, 𝑢 is not a fuzzy number.
The set of all fuzzy numbers in 𝑅 is denoted by 𝐹(𝑅). From (i)–(iv), it follows that if 𝑢 belongs to F(R) then the α -level set [𝑢]^𝛼 is a nonempty compact interval for all 𝛼∈[0,1]. The notation [𝑢]^𝛼=[𝑢_𝛼,𝑢_𝛼] denotes explicitly the α-level set of 𝑢. Triangular fuzzy numbers are one of the most commonly used fuzzy numbers. The membership function of a triangular fuzzy number 𝑢 is completely characterized with the peak (or center) m, left width 𝜎≥0, and right width 𝛽≥0 and has the following form: 𝜇_𝑢(𝑥)=1−𝑚−𝑥𝜎,if𝑚−𝜎≤𝑥<𝑚,1−𝑥−𝑚𝛽,if𝑚≤𝑥≤𝑚+𝛽,0,otherwise, (1) and we use the notation 𝑢=(𝑚,𝜎,𝛽). A triangular fuzzy number 𝑢=(𝑚,𝜎,𝛽) is named a symmetric triangular fuzzy number if its left width and right width are equal and we denote it by 𝑢=(𝑚,𝜎), for brevity. Some examples of triangular and symmetric triangular fuzzy numbers are depicted in Figures 2 and 3, respectively. We consider a crisp number 𝑎∈𝑅 as a symmetric triangular fuzzy number 𝑎=(𝑎,0).
Let 𝑓 be a real-valued mapping on 𝑅^𝑛. Assume 𝑢_𝑗, 𝑗=1,…,𝑛 are fuzzy numbers in 𝑅. Using the extension principle, we can define 𝑌=𝑓(𝑢₁,…,𝑢_𝑛) as a fuzzy set in 𝑅 such that 𝜇_𝑌𝑦=supmin𝜇𝑢₁𝑥₁,…,𝜇𝑢_𝑛𝑥_𝑛∣𝑥₁,…,𝑥_𝑛∈𝑓⁻¹𝑦,if𝑓⁻¹𝑦≠𝜙,0,if𝑓⁻¹𝑦=𝜙. (2)

Figure 1 

The membership function of the fuzzy set of real numbers “close to 1000” and its α-level set.

Figure 2 

The fuzzy numbers 𝑢=(4,4,1) (- - -), 𝑣=(3,1,4) (—), max{𝑢,𝑣} (∘∘), and min{𝑢,𝑣} (**).

Figure 3 

The fuzzy numbers 𝑢=(−1,1) (—), 𝑣=(2,2/3) (- - -) , 𝑢Θ𝐻𝑣 (⋯), and 𝑢⊗𝑍𝑣 (–·–).

Example 3. Let 𝑓(𝑥)=𝜆𝑥 be a linear function. Suppose 𝑢∈𝐹(𝑅), and let 𝑌=𝑓(𝑢). Then, using the extension principle, we obtain 𝜇_𝑌(𝑦)=sup{𝜇_𝑢(𝑥)∣𝑦=𝜆𝑥}=𝜇_𝑢(𝑦/𝜆). Especially, if 𝜆=−1, then we write 𝑓(𝑢)=Θ𝑢 and we have 𝜇_Θ𝑢(𝑦)=𝜇_𝑢(−𝑦).

Theorem 4. Let 𝑓 be a real-valued mapping on 𝑅^𝑛, and let 𝑢_𝑗, 𝑗=1,…,𝑛, be fuzzy numbers in 𝑅. Then, [𝑓(𝑢₁,…,𝑢_𝑛)]^𝛼=𝑓([𝑢₁]^𝛼,…,[𝑢_𝑛]^𝛼) where 𝑓(𝑢₁,…,𝑢_𝑛) is defined by the extension principle and 𝑓([𝑢₁]^𝛼,…,[𝑢_𝑛]^𝛼)={𝑓(𝑢₁,…,𝑢_𝑛)∣𝑢₁∈[𝑢₁]^𝛼,…,𝑢_𝑛∈[𝑢_𝑛]^𝛼}.

Proof (see [12]). The above theorem defines arithmetic operations of fuzzy numbers in terms of their α-level sets by [𝑢⊕𝑣]^𝛼=[𝑢_𝛼+𝑣_𝛼,𝑢_𝛼+𝑣_𝛼], [𝑢Θ𝑣]^𝛼=[𝑢_𝛼−𝑣_𝛼,𝑢_𝛼−𝑣_𝛼], and [𝑢⊗𝑣]^𝛼=[min{𝑢_𝛼𝑣_𝛼,𝑢_𝛼𝑣_𝛼,𝑢_𝛼𝑣_𝛼,𝑢_𝛼𝑣_𝛼}, max{𝑢_𝛼𝑣_𝛼,𝑢_𝛼𝑣_𝛼,𝑢_𝛼𝑣_𝛼,𝑢_𝛼𝑣_𝛼}], where ⊕, Θ, and ⊗ denote the addition, minus, and multiplication operators on 𝐹(𝑅), respectively. Moreover, [Θ𝑢]^𝛼=[−𝑢_𝛼,−𝑢_𝛼].

Definition 5 (fuzzy max, fuzzy min, fuzzy inequality, and the weighted center of gravity). Let 𝑢,𝑣∈𝐹(𝑅), and 𝑓(𝑥,𝑦)=max{𝑥,𝑦}. Then, the maximum of 𝑢 and 𝑣 is defined by 𝑓(𝑢,𝑣) and applying the extension principle. Similarly, setting 𝑓(𝑥,𝑦)=min{𝑥,𝑦}, then 𝑓(𝑢,𝑣) defines the minimum of 𝑢 and 𝑣. We denote the maximum and the minimum of 𝑢 and 𝑣 by max{𝑢,𝑣} and min{𝑢,𝑣}, respectively. The notation “≤” will be used for the inequality relation between fuzzy numbers 𝑢 and 𝑣 and is defined as 𝑢≤𝑣⇔miñ{𝑢,𝑣}=𝑢. Moreover, the quantity of fuzzy number 𝑢 can be given by its weighted center of gravity (WCOG) defined as 𝑢=∫¹₀𝛼(𝑢_𝛼+𝑢_𝛼)𝑑𝛼, where the weights are the membership degrees.

Example 6. Two triangular fuzzy numbers 𝑢=(4,4,1), 𝑣=(3,1,4) and their maximum and minimum are shown in Figure 2. Moreover, if 𝑢=(𝑚,𝜎,𝛽), then the WCOG of 𝑢 is 𝑢=∫¹₀𝛼[2𝑚+(1−𝛼)(𝜎−𝛽)]𝑑𝛼=𝑚+1/6(𝜎−𝛽).
Applying Theorem 4, we get [miñ(𝑢,𝑣)]^𝛼=min{[𝑢]^𝛼,[𝑣]^𝛼}=[min{𝑢_𝛼,𝑣_𝛼},min{𝑢_𝛼,𝑣_𝛼}]; hence, 𝑢≤𝑣⇔𝑢_𝛼≤𝑣_𝛼, 𝑢_𝛼≤𝑣_𝛼, for all 𝛼∈[0,1].

Definition 7 (H-difference and Z-product). Let 𝑢,𝑣∈𝐹(𝑅). If there exists 𝑧∈𝐹(𝑅) such that 𝑢=𝑣⊕𝑧, then 𝑧 is called the H-difference of 𝑢 and 𝑣 and it is denoted by 𝑢Θ_𝐻𝑣. Moreover, if there exist 𝑧∈𝐹(𝑅) such that 𝑧⊗𝑣⁻¹=𝑢, then we call it the Z-product of 𝑢 and 𝑣 and we denote it by 𝑧=𝑢⊗_𝑍𝑣.
Note that the Z-product of fuzzy numbers is a new concept which is introduced in this paper for the first time. It is easy to see that [𝑢Θ_𝐻𝑣]^𝛼=[𝑢_𝛼−𝑣_𝛼,𝑢_𝛼−𝑣_𝛼] and [𝑢⊗_𝑍𝑣]^𝛼=[𝑢_𝛼𝑣_𝛼,𝑢_𝛼𝑣_𝛼], if 𝑢≤(0,0) and (0,0)≤𝑣.

Example 8. Let 𝑢=(−1,1) and 𝑣=(2,2/3). A straightforward calculation shows [𝑢Θ_𝐻𝑣]^𝛼=[−3−1/3(1−𝛼),−3+1/3(1−𝛼)]; hence, 𝑢Θ_𝐻𝑣=(−3,1/3). Furthermore, we obtain [𝑢⊗_𝑍𝑣]^𝛼=[−8/3+2/3𝛼²,−8/3𝛼+2/3𝛼²] which is valid α-level set of a fuzzy number (not a triangular fuzzy number). The fuzzy numbers 𝑢, 𝑣 and their H-difference and Z-product are shown in Figure 3, for the sake of clarity.

Definition 9 (the strongly generalized differentiability). Let 𝑥∶(𝑎,𝑏)→𝐹(𝑅) be a fuzzy function and 𝑡₀∈(𝑎,𝑏). We say that. 𝑥 is differentiable at 𝑡₀ if it exists an element 𝑥̇(𝑡₀)∈𝐹(𝑅) such that, for all ℎ>0 sufficiently near to 0,(i)there are 𝑥(𝑡₀+ℎ)Θ_𝐻𝑥(𝑡₀) and 𝑥(𝑡₀)Θ_𝐻𝑥(𝑡₀−ℎ) and the limits lim_ℎ→0+(𝑥(𝑡₀+ℎ)Θ_𝐻𝑥(𝑡₀)/ℎ)=lim_ℎ→0+(𝑥(𝑡₀)Θ_𝐻𝑥(𝑡₀−ℎ)/ℎ)=𝑥̇(𝑡₀) or(ii)there are 𝑥(𝑡₀−ℎ)Θ_𝐻𝑥(𝑡₀) and 𝑥(𝑡₀)Θ_𝐻𝑥(𝑡₀+ℎ) and the limits lim_ℎ→0+(𝑥(𝑡₀)Θ_𝐻𝑥(𝑡₀+ℎ)/−ℎ)=lim_ℎ→0+(𝑥(𝑡₀−ℎ)Θ_𝐻𝑥(𝑡₀)/−ℎ)=𝑥̇(𝑡₀), where the limits are taken in the metric 𝐷 defined as 𝐷(𝑢,𝑣)=sup_𝛼∈[0,1]max{|𝑢_𝛼−𝑣_𝛼|,|𝑢_𝛼−𝑣_𝛼|}, for all 𝑢,𝑣∈𝐹(𝑅).

Theorem 10. If 𝑥(𝑡) is differentiable in the first form (i), then 𝑥_𝛼(𝑡) and 𝑥_𝛼(𝑡) are differentiable functions and [𝑥̇(𝑡)]^𝛼=[𝑥̇_𝛼(𝑡),𝑥̇_𝛼(𝑡)]. If 𝑥(𝑡) is differentiable in the second form (ii), then 𝑥_𝛼(𝑡) and 𝑥_𝛼(𝑡) are differentiable functions and [𝑥̇(𝑡)]^𝛼=[𝑥̇_𝛼(𝑡),𝑥̇_𝛼(𝑡)].

Proof (see [17]). Consider the following fuzzy initial value problem: 𝑥̇(𝑡)=𝑓(𝑡,𝑥(𝑡)),𝑥(0)=𝑥₀, (3) where 𝑓∶𝐼×𝐹(𝑅)→𝐹(𝑅) is a fuzzy function and 𝑥₀∈𝐹(𝑅). Let [𝑓(𝑡,𝑥)]^𝛼=[𝑓_𝛼(𝑡,𝑥_𝛼,𝑥_𝛼),𝑓_𝛼(𝑡,𝑥_𝛼,𝑥_𝛼)]. From Theorem 10, if we consider 𝑥(𝑡) by using the derivative in the first form (i), then the solution of problem (3) is obtained by solving the following system of ordinary differential equations: 𝑥̇_𝛼(𝑡)=𝑓_𝛼𝑡,𝑥_𝛼(𝑡),𝑥_𝛼(𝑡),𝑥_𝛼(0)=𝑥_0𝛼,𝑥̇_𝛼(𝑡)=𝑓_𝛼𝑡,𝑥_𝛼(𝑡),𝑥_𝛼(𝑡),𝑥_𝛼(0)=𝑥_0𝛼, (4) and ensuring that [𝑥_𝛼(𝑡),𝑥_𝛼(𝑡)] and [𝑥̇_𝛼(𝑡),𝑥̇_𝛼(𝑡)] are valid α-level sets. Moreover, if we consider 𝑥(𝑡) by using the derivative in the second form (ii), then the solution of problem (3) is obtained by solving the following system of ordinary differential equations, 𝑥̇‌_𝛼(𝑡)=𝑓_𝛼𝑡,𝑥_𝛼(𝑡),𝑥_𝛼(𝑡),𝑥_𝛼(0)=𝑥_0𝛼,𝑥̇_𝛼(𝑡)=𝑓_𝛼𝑡,𝑥_𝛼(𝑡),𝑥_𝛼(𝑡),𝑥_𝛼(0)=𝑥_0𝛼, (5) and ensuring that [𝑥_𝛼(𝑡),𝑥_𝛼(𝑡)]and [𝑥̇_𝛼(𝑡),𝑥̇_𝛼(𝑡)] are valid α-level sets.
The integral of fuzzy function 𝑥(𝑡) using the Riemann integral concept can be defined as follows.

Definition 11. The integral of a fuzzy mapping 𝑥∶[𝑎,𝑏]→𝐹(𝑅) is defined levelwise by ^𝑏_𝑎𝑥(𝑡)𝑑𝑡^𝛼=^𝑏_𝑎𝑥_𝛼(𝑡)𝑑𝑡,^𝑏_𝑎𝑥_𝛼(𝑡)𝑑𝑡. (6) Note that if 𝑥∶[𝑎,𝑏]→𝐹(𝑅) is continuous in the metric D, then it is integrable, that is, ∫^𝑏_𝑎𝑥(𝑡)𝑑𝑡∈𝐹(𝑅).

Example 12. Define the fuzzy mapping 𝑥∶[0,0.5]→𝐹(𝑅) by 𝑥(𝑡)=(𝑡,sin𝜋𝑡,cos𝜋𝑡). Then, [∫^𝑏_𝑎𝑥(𝑡)𝑑𝑡]^𝛼=[∫₀^0.5(𝑡−(1−𝛼)sin𝜋𝑡)𝑑𝑡,∫₀^0.5(𝑡+(1−𝛼)cos𝜋𝑡)𝑑𝑡]=[0.125−(1−𝛼)𝜋⁻¹,0.125+(1−𝛼)𝜋⁻¹]; hence, ∫^𝑏_𝑎𝑥(𝑡)𝑑𝑡=(0.125,𝜋⁻¹).

3. Linear Fuzzy Model of HIV Infection

HIV infection can be characterized as a disease of the immune system, with progressive depletion of defensive cells, resulting in immunosuppression and susceptibility to opportunistic infections. CD4+ T-cells, CTLs, and the virus particles play important roles in HIV infection. CD4+ T-cells are a fundamental component of the human immune response system. These cells can be considered “messengers” or the command centers of the immune system, and they signal other immune cells that an invader is to be fought. The immune response cells, or cytotoxic lymphocytes (CTLs), are the cells that respond to this message and set out to eliminate infection by killing infected cells. HIV can infect a number of cells in the body however, its main target is the CD4+ T-cells. HIV enters these cells by a complex process and begins to replicate, then the new virus particles are released by bursting the infected cells. CD4+ T-cells are generated from sources within the body and are lost either by having finite life span or by bursting during the proliferation of HIV, which leads to a drop in the number of these cells, after infection and an accelerated decrease during the later stages of the disease that signals the onset AIDS. In accordance with experimental findings, too high a level of HIV impairs establishment of a lasting CTL response. This is a delicate task, since CD4+ T-cell population, which plays an essential role in stimulation of immune response, depletes dramatically with raising the HIV load. The rate of CD4+ T-cells depletion varies greatly from patient to patient, depending on the strength or weakness of the immune system. More precisely, a stronger immune system leads to a lower rate of CD4+ T-cells depletion and vice versa. We have a similar argument about the proliferation rate of HIV particles. Therefore, the levels of the immune cells as well as the HIV viral load during the different stages of the disease can be considered as fuzzy quantities. According to these descriptions, the interaction of HIV with the immune system can be modeled by a system of linear differential equations with fuzzy parameters as follows:̇̇̇𝑥=𝜆Θ𝜎⊗𝑥Θ𝑐̃⊗𝑣,𝑣=𝑘⊗𝑣Θ𝑎⊗𝑧,𝑧=ℎ⊗𝑥Θ𝜏⊗𝑣, (7) where the fuzzy functions 𝑥(𝑡), 𝑧(𝑡), and 𝑣(𝑡) indicate the level of CD4+ T-cells, CTLs, and the HIV viral load at time t, respectively. Most of the terms in the model have straightforward interpretations as follows.

The first equation in (7) represents the dynamics of the concentration of CD4+ T-cells. The CD4+ T-cells are produced from a source, such as the thymus, at a constant rate 𝜆. Here, we have assumed that CD4+ T-cells have a finite life-span and die at a rate 𝜎 per cell. Therefore, the number of these cells, which are lost due to natural death, is represented through the loss term 𝜎⊗𝑥 in the first equation. Moreover, the CD4+ T-cell population is lost through infection by a virus particle at a rate of 𝑐̃, and so the term 𝑐̃⊗𝑣 models the rate that free viruses destroy CD4+ T-cells. The second equation in (7) depicts the rate of change in the virus population. An HIV particle uses a host cell to replicate itself and thus proliferates with a growth rate 𝑘. Thus, the total amount of produced viruses is given by the term 𝑘⊗𝑣. Infected cells are killed by CTLs, and hence viruses are lost through an immune response. Assuming that a CTL eliminates the virus particles at a rate 𝑎, the number of virus particles eliminated by the immune response is given by the term 𝑎⊗𝑧. The third equation in (7) describes the dynamics of CTLs during HIV infection. A CD4+ T-cell stimulates CTLs to proliferate at a rate ℎ. Therefore, CD4+ T-cells effect on proliferation of CTLs is expressed by the term ℎ⊗𝑥. The term 𝜏⊗𝑣 takes into account loss of CTLs due to increasing the HIV viral load where 𝜏 is the rate at which the virus-induced impairment of CD4+ T-cell function occurs.

In this paper, a patient with respect to the strength or weakness of its immune system is considered as a patient with the weak, moderate, or strong immune system and is indicated by 𝑊, 𝑀, or 𝑆, respectively. The initial condition of 𝑣(0)=𝑣₀ varies in different patients. For that reason, a primary response is provoked when the immune system encounters HIV for the first time and, in this stage, a number of viruses depending on the strength or weakness of the immune system are eliminated, and the proposed model describes the changes in the immune cells level and the viral load after this stage which is called the secondary immune response. Therefore, a stronger immune response implies a lower 𝑣₀ and vice versa. The values of the model parameters and 𝑣₀ corresponding to patients 𝑊, 𝑀, and 𝑆 are shown in Figure 4 as triangular fuzzy numbers. These parameters were chosen to be consistent with biological plausibility. Moreover, we assume that, at time 𝑡=0, the level of CD4+ T-cells is normal and there is no CTL-mediated immune response in all patients, that is, 𝑥(0)=𝑥₀=(100,0) and 𝑧(0)=𝑧₀=(0,0). We must note that 𝑥(𝑡) denotes CD4+ T-cells level in percentage at time 𝑡. The derivative in the second form (ii) leads to solutions with decreasing length of their support which leads us to the conclusion that the uncertainty decreases with the time lapse which is not consistent with real situation. Moreover, the existence of these solutions implies that the initial conditions should be fuzzy. Therefore, we consider only the solutions with the derivative in the first forms (i) which are more consistent with real situation. Consequently, as mentioned in Section 2, the fuzzy model (7) is transformed to the following system of ordinary differential equations (ODEs): 𝑥̇‌_𝛼(𝑡)=𝜆_𝛼−𝜎_𝛼𝑥_𝛼(𝑡)−𝑐_𝛼𝑣_𝛼(𝑡),𝑥̇_𝛼(𝑡)=𝜆_𝛼−𝜎_𝛼𝑥_𝛼(𝑡)−𝑐_𝛼𝑣_𝛼(𝑡),𝑣̇_𝛼(𝑡)=𝑘_𝛼𝑣_𝛼(𝑡)−𝑎_𝛼𝑧_𝛼(𝑡),𝑣̇_𝛼(𝑡)=𝑘_𝛼𝑣_𝛼(𝑡)−𝑎_𝛼𝑧_𝛼(𝑡),𝑧̇_𝛼(𝑡)=ℎ_𝛼𝑥_𝛼(𝑡)−𝜏_𝛼𝑣_𝛼(𝑡),𝑧̇_𝛼(𝑡)=ℎ_𝛼𝑥_𝛼(𝑡)−𝜏_𝛼𝑣_𝛼(𝑡),𝑥_𝛼(0)=𝑥_0𝛼,𝑥_𝛼(0)=𝑥_0𝛼,𝑣_𝛼(0)=𝑣_0𝛼,𝑣_𝛼(0)=𝑣_0𝛼,𝑧_𝛼(0)=𝑧_0𝛼,𝑧_𝛼(0)=𝑧_0𝛼. (8) For each 𝛼∈[0,1], the ODEs (8) are linear; hence, the exact solutions in discrete times are obtained using the ode45 code in MATLAB. However, it would be appropriate to propose explicit solutions as a function of 𝛼 and 𝑡. The next section is devoted to this topic.

Figure 4 

The values of model parameters corresponding to patients 𝐖 (- -), 𝐌 (⋯), and 𝐒 (—).

3.1. The Approximate Explicit Solutions Based on a Fitting Method

The proposed method is based on the fact that a linear combination of suitable functions of α and 𝑡 can generate the best fit to the exact values obtained by the ode45 in the least squares sense. The following discussion shows that these functions can be exponential. The ODEs (8) can be written in a matrix form as: 𝑋̇_𝛼(𝑡)=𝐴_𝛼𝑋_𝛼(𝑡)+𝐵_𝛼,𝑋_𝛼(0)=𝑋_0𝛼. (9) By the variation of constants formula for ordinary differential equations, the solution of the initial value problem (9) is 𝑋_𝛼(𝑡)=𝑒^𝐴𝛼𝑡𝑋_0𝛼+∫^𝑡₀𝑒^{𝐴𝛼(𝑡−𝜏)}𝐵_𝛼𝑑𝜏. Since the six-dimensional matrix 𝐴_𝛼 depends on α, the calculation of 𝑒^𝐴𝛼𝑡 becomes difficult. But this matrix can be written as 𝐴_𝛼=𝐴₁+(1−𝛼)𝐴₂, where 𝐴₁ and 𝐴₂ are α-independent matrices. We have 𝑒^𝐴𝛼𝑡=𝑒^𝐴1𝑡𝑒^{(1−𝛼)𝐴2𝑡}+𝑂(𝑡), where 𝑂(𝑡) is a function that lim_𝑡→0𝑂(𝑡)/𝑡=0. Assuming that ℓ_𝑗s and 𝜍_𝑗𝑠,𝑗=1,…,6, are eigenvalues of 𝐴₁ and 𝐴₂, respectively, there are invertible matrices 𝑃 and 𝑄 such that 𝐴₁=𝑃𝐷₁𝑃⁻¹ and 𝐴₂=𝑄𝐷₂𝑄⁻¹, where 𝐷₁=diag(ℓ₁,…,ℓ₆) and 𝐷₂=diag(𝜍₁,…,𝜍₆). Therefore, for small 𝑡,𝑒^𝐴𝛼𝑡𝑒≈𝑃^ℓ1𝑡0⋯00𝑒^ℓ2𝑡⋯⋮⋮⋮⋱00⋯0𝑒^ℓ6𝑡𝑃⁻¹𝑒×𝑄^{(1−𝛼)𝜁1𝑡}0⋯00𝑒^{(1−𝛼)𝜁2𝑡}⋯⋮⋮⋮⋱00⋯0𝑒^{(1−𝛼)𝜁6𝑡}𝑄⁻¹. (10) As a result, the fitting functions are chosen as 𝑒^{(ℓ𝑖+(1−𝛼)𝜍𝑗)𝑡}, 𝑖,𝑗=1,…,6. Therefore, by choosing the numbers ℓ_𝑗, 𝑗=1,…,𝑛, and 𝜍_𝑖, 𝑖=1,…,𝑚, where 𝑛,𝑚∈{1,2,…,6}, an approximate solution can be found in the following form:𝑆_𝛼(𝑡)≈𝐾𝐸_𝛼(𝑡), (11) where 𝑆_𝛼(𝑡)=[𝑥_𝛼(𝑡),𝑥_𝛼(𝑡),𝑣_𝛼(𝑡),𝑣_𝛼(𝑡)]^𝑇, 𝐸_𝛼(𝑡)=[𝑒^{(ℓ1+(1−𝛼)𝜍1)𝑡},…,𝑒^{(ℓ1+(1−𝛼)𝜍𝑚)𝑡},…,𝑒^{(ℓ𝑛+(1−𝛼)𝜍1)𝑡},…,𝑒^{(ℓ𝑛+(1−𝛼)𝜍𝑚)𝑡}]^𝑇 and 𝐾=[𝑘_𝑖𝑗]_4×𝑛𝑚 denotes the coefficients matrix that can be found using the lsqnonlin code of the optimization toolbox in MATLAB. Obviously, 𝑧_𝛼(𝑡) and 𝑧_𝛼(𝑡) are obtained from the 5th and the 6th equations in (8) using the approximate values 𝑥_𝛼(𝑡), 𝑥_𝛼(𝑡), 𝑣_𝛼(𝑡), and 𝑣_𝛼(𝑡).

3.2. Dynamic Behavior of the Immune Cells Level and the Viral Load in Patient 𝑊

Figure 5 shows the level of immune cells and the HIV viral load of patient 𝑊 during the time interval [0, 1800]. The darker color shows the curve with the higher possibility. Figure 5(a) shows that the gradual declines in CD4+ T-cells level correspond to the low possibilities, while the rapid declines and the progression to full blown AIDS after a gradual decay, have the high possibilities of occurrence. Figure 5(b) shows a rapid increase in the viral load during the later stages of the disease. Moreover, we observe that, at each time, a higher viral load corresponds to a higher possibility. Figure 5(c) shows a clear correlation between CTLs level in the blood and HIV progression. As the viral load increases upon initial infection, CTLs increase in order to decrease the virus. But ultimately the level of these cells begins to decrease, which is due to virus-induced impairment of CD4+ T-cell function, with the high possibilities after about the 1450th day. Besides, a lower CTLs level has a higher possibility of occurrence and vice versa. From Figures 5(a) and 5(b), there is an inverse correlation between the HIV viral load and the level of CD4+ T-cells. Following the proposed method in Section 3.1, choose the eigenvalues of the corresponding matrices 𝐴₁ and 𝐴₂ as 𝜍₁,𝜍₂=−0.2600×10⁻³,0.6537×10⁻¹¹,ℓ₁,ℓ₂,ℓ₃,ℓ₄,ℓ₅=−0.2801×10⁻²,0.2801×10⁻²,−0.2145×10⁻⁶,−0.1128×10⁻⁶,0.6760×10⁻². (12) Then, we have an approximate explicit solution in the form of (11) where the corresponding coefficients matrix 𝐾 is𝐾=.0.668410.1635−0.22450.2229−0.710045.2900−0.702845.2924−0.00090.0005−3.140613.97790.1989−0.19740.676343.89730.683243.89940.0004−0.0008−0.01230.0109−0.00020.00020.0044−0.00380.0044−0.0038−0.00000.00030.8273−0.82220.0219−0.0217−0.16240.1608−0.16230.16090.0004−0.0001 (13)

(a)

(b)

(c)

Figure 5 

CD4+ T-cells level (a), the viral load (b), and CTLs level (c) versus time in patient W.

A comparison between the exact and approximate solutions which is shown in Figure 6 confirmed the effectiveness of this approach.

(a)

(b)

Figure 6 

Exact (—) and approximate (×) CD4+ T-cells level (a) and the viral load (b) in patient 𝑊.

3.3. Dynamic Behavior of the Immune Cells Level and the Viral Load in Patient 𝑆

Figure 7 shows the changes in the immune cells level and the viral load in patient 𝑆. With respect to Figure 7(a), CD4+ T-cells level decreases gradually during the 1800 days from infection with the high possibilities. Moreover, an increase in CD4+ T-cells as well as the progression to AIDS arises with a low possibility. The HIV viral load is low and a lower viral load has a higher possibility of occurring, as shown in Figure 7(b). A high CD4+ T-cell count and a low HIV viral load lead to establishment of a lasting CTL response which is shown in Figure 7(c). A high HIV viral load and a low CD4+ T-cells level impair the immune response where this arises with the low possibilities as shown in Figure 7 as the light curves. The selected numbers ℓ_𝑖𝑠 and 𝜍_𝑗𝑠 and the corresponding coefficients matrix 𝐾 representing the approximate explicit solutions in the form of (11) are as𝜍₁,𝜍₂=0.2600×10⁻³,0.6537×10⁻¹¹,ℓ₁,ℓ₂,ℓ₃,ℓ₄,ℓ₅=−0.2800×10⁻²,0.2800×10⁻²,−0.1159×10⁻⁶,0.6240×10⁻²,−0.2378×10⁻⁶,𝐾=.−3.320313.91750.1287−0.12650.766243.92930.0003−0.00060.772343.94120.83039.7266−0.14310.1412−0.806045.5176−0.00080.0005−0.799645.53000.3895−0.39020.0110−0.0109−0.07240.07250.0004−0.0002−0.07240.07260.0076−0.00990.0002−0.0002−0.00510.0063−0.00000.0002−0.00510.0063 (14)

(a)

(b)

(c)

Figure 7 

CD4+ T-cells level (a), the viral load (b), and CTLs level (c) versus time in patient S.

3.4. Dynamic Behavior of the Immune Cells and the Viral Load in Patient 𝑀

Figure 8 shows the changes in the immune cells level and the viruses in patient 𝑀. From Figure 8(a), the uncertainty of CD4+ T-cells level increases and new possibilities, varying fromm an increment to normal level to rapid progression to full blown AIDS, arise after a gradual decay. But the most possible scenario is between where the level of these cells decreases at a moderate rate. A moderate viral load occurs with a high possibility, as shown in Figure 8(b). A moderate CD4+ T-cells level as well as a moderate viral load implies a moderate CTLs level, as shown in Figure 8(c). Besides, from this figure, a high (low) viral load and a low (high) CD4+ T-cells number decrease (increase) CTLs level, and this happens with a low possibility.

(a)

(b)

(c)

Figure 8 

CD4+ T-cells level (a), the viral load (b), and CTLs level (c) versus time in patient M.

Here, a representation for solutions is given in the form of (11) whereℓ₁,ℓ₂,ℓ₃,ℓ₄,ℓ₅=−0.2800×10⁻²,0.2800×10⁻²,0.2257×10⁻⁶,−0.1143×10⁻⁶,0.6500×10⁻²,𝜍₁,𝜍₂,𝜍₃=0,−0.2600×10⁻³,0.2600×10⁻³,𝐾=.45.7977−20.3001−14.92030.8002−0.6250−0.177226.51149.60398.573426.51889.61698.5785−0.0004−0.00060.000625.4227−12.6502−2.18780.4622−0.0933−0.370531.24918.39455.039631.25498.40975.0415−0.00000.0005−0.0008−0.2693−0.30820.5693−0.0159−0.00160.0174−0.00680.0869−0.0772−0.00690.0869−0.0773−0.00030.00010.0005−0.05170.2515−0.2007−0.00350.0074−0.0040−0.0033−0.03630.0406−0.0033−0.03630.0406−0.00020.00040.0000 (15)

3.5. A Comparison between the Immune Cells Level and the Viral Load in Patients 𝑊, 𝑀, and 𝑆

From Figures 5(b), 7(b), and 8(b), there is an inverse correlation between the viral load and the immune system strength. CD4+ T-cell and CTL levels in patient 𝑆 are more than the level of these cells in patient W, as shown in Figures 5 and 7. A high viral load in patient 𝑊 leads to a virus-induced impairment of CD4+ T-cell function. Therefore, CTLs level in this patient is less than the level of these cells in patient S, as shown in Figures 5(c) and 7(c). A comparison between CTLs level in patients 𝑀 and 𝑆 shows that it is possible that CTLs level in patient 𝑀 be slightly more than the level of these cells in patient S, which can be due to this fact that patient 𝑀 has more antigens that are required to stimulate CTLs (see Figures 7(b) and 8(b)), and, thus, a higher level of CD4+ T-cells in patient 𝑀 is possible as shown in Figures 7(a) and 8(a). Figure 9 shows CD4+ T-cell and CTL levels and the viral load of patients W, M, and 𝑆 on the 1800th day. Here, we perform a comparison between the immune cells level and the HIV viral load of patients based on their weighted center of gravity (WCOG).

(a)

(b)

(c)

Figure 9 

The fuzzy numbers indicating CD4+ T-cells level (a), the viral load (b), and CTLs level (c) on the 1800th day in patients 𝑊 (- -), 𝑀 (⋯) and, 𝑆 (—).

For this end, the WCOG of 𝑥(𝑡), 𝑣(𝑡), and 𝑧(𝑡) is denoted by 𝑥(𝑡), 𝑣(𝑡), and 𝑧(𝑡), respectively; hence,𝑥(𝑡)=¹₀𝛼𝑥_𝛼(𝑡)+𝑥_𝛼(𝑡)𝑣(𝑡)=𝑑𝛼,₁₀𝛼𝑣_𝛼(𝑡)+𝑣_𝛼(𝑡)𝑧(𝑡)=𝑑𝛼,₁₀𝛼𝑧_𝛼(𝑡)+𝑧_𝛼(𝑡)𝑑𝛼. (16) With respect to (11), an approximate explicit formula for the WCOG of 𝑥(𝑡) and 𝑣(𝑡) is given as 𝑆(𝑡)≈𝑁𝐹(𝑡), where 𝑆(𝑡)=[𝑥(𝑡),𝑣(𝑡)]^𝑇and 𝑁 is a matrix with two rows that its first and second rows are obtained by summing the first two rows and the last two rows of the corresponding coefficients matrix K, respectively. Moreover,𝐹(𝑡)=¹₀𝛼𝐸_𝛼(𝑡)=𝑒𝑑𝛼_ℓ1𝑡𝜍₁𝑡𝑒^𝜍1𝑡𝜍₁𝑡−1𝜍₁𝑡𝑒−1,…,_ℓ1𝑡𝜍_𝑚𝑡𝑒^𝜍𝑚𝑡𝜍_𝑚𝑡−1𝜍_𝑚𝑡𝑒−1,…,_ℓ𝑛𝑡𝜍₁𝑡𝑒^𝜍1𝑡𝜍₁𝑡−1𝜍₁𝑡𝑒−1,…,_ℓ𝑛𝑡𝜍_𝑚𝑡𝑒^𝜍𝑚𝑡𝜍_𝑚𝑡−1𝜍_𝑚𝑡−1^𝑇. (17) The WCOG of 𝑥(𝑡), 𝑣(𝑡), and 𝑧(𝑡) corresponding to patients 𝑊, 𝑀, and 𝑆 is shown in Figure 10. From Figures 10(a) and 10(c), CD4+ T-cells and CTLs level are proportional to the strength of patient’s immune system such that a stronger immune system leads to a higher level of these cells. Moreover, there is an inverse correlation between the viral load and the strength of the immune system as shown in Figure 10(b).

(a)

(b)

(c)

Figure 10 

The WCOG of CD4+ T-cells level (a), the viral load (b) and CTLs level (c) versus time in patients 𝑊 (- -), 𝑀 (⋯), and 𝑆 (—).

4. Fuzzy Optimal Control Problem

In this section, we formulate a fuzzy optimal control problem that identifies the parameter 𝑘 in (8), with a function of the fuzzy control variable 𝑢. In particular, we will replace the parameter 𝑘 with the function 𝑘Θ_𝐻𝑢. This choice then identifies the control variable 𝑢(𝑡) with the rate of inhibition of virus reproduction, which is modeled as a simple function of drug dosage. Therefore, we have the fuzzy set control differential equations (FSCDEs) aṡ̇𝑥=𝜆Θ𝜎⊗𝑥Θ𝑐̃⊗𝑣,𝑣=𝑘Θ_𝐻̇𝑢⊗𝑣Θ𝑎⊗𝑧,𝑧=ℎ⊗𝑥Θ𝜏⊗𝑣. (18) This paper aims to propose a drug regimen that minimizes both the viral load and the drug costs. Here, we assume that the cost of the treatment is proportional to 𝑢²(𝑡) at time 𝑡. Therefore, the fuzzy functional 𝐽(̃𝑣,𝑢)=∫^𝑡𝑓_𝑡0(𝑤⊗𝑣(𝑡)⊕𝑢²(𝑡))𝑑𝑡 should be minimized, where the positive fuzzy number 𝑤 is used to set the relative importance between minimizing the viral load and the systemic cost to the body. Let 𝑇=[𝑡₀,𝑡_𝑓], and assume 𝐶(𝑇) be the set of all continuous fuzzy number valued functions on 𝑇. Assuming that the minimum and the maximum of allowable drug dosage are denoted by fuzzy numbers 𝑢^min and 𝑢^max, respectively, then we are seeking a 𝑢^∗∈𝑈 such that 𝐽(̃𝑣^∗,𝑢^∗)≤𝐽(̃𝑣,𝑢), for all 𝑢∈𝑈, where 𝑈={𝑢∈𝐶(𝑇)∶𝑢^min≤𝑢(𝑡)≤𝑢^max, for all 𝑡∈𝑇}. By using the derivative in the first form (i), the FSCDEs (18) is converted to the following control system:𝑥̇‌_𝛼(𝑡)=𝜆_𝛼−𝜎_𝛼𝑥_𝛼(𝑡)−𝑐_𝛼𝑣_𝛼(𝑡),𝑥̇_𝛼(𝑡)=𝜆_𝛼−𝜎_𝛼𝑥_𝛼(𝑡)−𝑐_𝛼𝑣_𝛼(𝑡),𝑣̇_𝛼(𝑡)=𝑘_𝛼−𝑢_𝛼(𝑡)𝑣_𝛼(𝑡)−𝑎_𝛼𝑧_𝛼(𝑡),𝑣̇_𝛼(𝑡)=𝑘_𝛼−𝑢_𝛼(𝑡)𝑣_𝛼(𝑡)−𝑎_𝛼𝑧_𝛼(𝑡),𝑧̇_𝛼(𝑡)=ℎ_𝛼𝑥_𝛼(𝑡)−𝜏_𝛼𝑣_𝛼(𝑡),𝑧̇_𝛼(𝑡)=ℎ_𝛼𝑥_𝛼(𝑡)−𝜏_𝛼𝑣_𝛼(𝑡),𝑥_𝛼(0)=𝑥_0𝛼,𝑥_𝛼(0)=𝑥_0𝛼,𝑣_𝛼(0)=𝑣_0𝛼,𝑣_𝛼(0)=𝑣_0𝛼,𝑧_𝛼(0)=𝑧_0𝛼,𝑧_𝛼(0)=𝑧_0𝛼. (19) With respect to Definition 5, a fuzzy function 𝑢^∗∈𝑈 is viewed as an optimal solution, if, for each 𝛼∈[0,1], the pair (𝑢^∗_𝛼,𝑢^∗_𝛼)∈𝑈 minimizes the functionals [∫^𝑡𝑓_𝑡0(𝑤⊗𝑣(𝑡)⊕𝑢²(𝑡))𝑑𝑡]_𝛼=∫^𝑡𝑓_𝑡0(𝑤_𝛼𝑣_𝛼(𝑡)+𝑢²_𝛼(𝑡))𝑑𝑡 and [∫^𝑡𝑓_𝑡0(𝑤⊗𝑣(𝑡)⊕𝑢²(𝑡))𝑑𝑡]_𝛼=∫^𝑡𝑓_𝑡0(𝑤_𝛼𝑣_𝛼(𝑡)+𝑢²_𝛼(𝑡))𝑑𝑡 simultaneously, where 𝑣_𝛼(𝑡) and 𝑣_𝛼(𝑡) are the solution of ODEs (19) corresponding to control pair (𝑢_𝛼,𝑢_𝛼), and 𝑈 is the set of all measurable control pairs (𝑢_𝛼,𝑢_𝛼) that 𝑢_𝛼^min≤𝑢_𝛼(𝑡)≤𝑢_𝛼^max,𝑢_𝛼^min≤𝑢_𝛼(𝑡)≤𝑢_𝛼^max, for all 𝑡∈𝑇. Therefore, we restrict our attention to optimizing the functional 𝐽(𝑣_𝛼,𝑣_𝛼,𝑢_𝛼,𝑢_𝛼)=∫^𝑡𝑓_𝑡0(𝑤_𝛼𝑣_𝛼(𝑡)+𝑤_𝛼𝑣_𝛼(𝑡)+𝑢²_𝛼(𝑡)+𝑢²_𝛼(𝑡))𝑑𝑡 over the set 𝑈. We now proceed to compute candidates for an optimal solution by applying the Pontryagin’s Maximum Principle [24] and begin by defining the Lagrangian to be 𝐿_𝛼𝑥_𝛼,𝑥_𝛼,𝑣_𝛼,𝑣_𝛼,𝑧_𝛼,𝑧_𝛼,𝑢_𝛼,𝑢_𝛼,𝜆_1𝛼,…,𝜆_6𝛼=𝑤_𝛼𝑣_𝛼(𝑡)+𝑤_𝛼𝑣_𝛼(𝑡)+𝑢²_𝛼+𝑢²_𝛼+𝜆_1𝛼𝜆_𝛼−𝜎_𝛼𝑥_𝛼−𝑐_𝛼𝑣_𝛼+𝜆_2𝛼𝜆_𝛼−𝜎_𝛼𝑥_𝛼−𝑐_𝛼𝑣_𝛼+𝜆_3𝛼𝑘_𝛼−𝑢_𝛼𝑣_𝛼−𝑎_𝛼𝑧_𝛼+𝜆_4𝛼𝑘_𝛼−𝑢_𝛼𝑣_𝛼−𝑎_𝛼𝑧_𝛼+𝜆_5𝛼ℎ_𝛼𝑥_𝛼−𝜏_𝛼𝑣_𝛼+𝜆_6𝛼ℎ_𝛼𝑥_𝛼−𝜏_𝛼𝑣_𝛼−𝜔₁₁𝑢_𝛼−𝑢_𝛼^𝑚𝑖𝑛−𝜔₁₂𝑢_𝛼^𝑚𝑎𝑥−𝑢_𝛼−𝜔₂₁𝑢_𝛼−𝑢_𝛼^𝑚𝑖𝑛−𝜔₂₂𝑢_𝛼^𝑚𝑎𝑥−𝑢_𝛼, (20) where 𝜔_𝑖𝑗≥0 are the penalty multipliers satisfying 𝜔₁₁(𝑢^∗_𝛼−𝑢_𝛼^min)=𝜔₁₂(𝑢_𝛼^max−𝑢^∗_𝛼)=0 and 𝜔₂₁(𝑢^∗_𝛼−𝑢_𝛼^min)=𝜔₂₂(𝑢_𝛼^max−𝑢^∗_𝛼)=0. Thus, the Maximum Principle gives the existence of adjoint variables 𝜆_𝑗𝛼, 𝑗=1,…,6, satisfying 𝜆̇_1𝛼=−𝜕𝐿_𝛼𝜕𝑥_𝛼=𝜎_𝛼𝜆_2𝛼−ℎ_𝛼𝜆_5𝛼,𝜆̇_2𝛼=−𝜕𝐿_𝛼𝜕𝑥_𝛼=𝜎_𝛼𝜆_1𝛼−ℎ_𝛼𝜆_6𝛼,𝜆̇_3𝛼=−𝜕𝐿_𝛼𝜕𝑣_𝛼=−𝑤_𝛼+𝑐_𝛼𝜆_2𝛼−𝑘_𝛼−𝑢_𝛼𝜆_3𝛼+𝜏_𝛼𝜆_6𝛼,𝜆̇_4𝛼=−𝜕𝐿_𝛼𝜕𝑣_𝛼=−𝑤_𝛼+𝑐_𝛼𝜆_1𝛼−𝑘_𝛼−𝑢_𝛼𝜆_4𝛼+𝜏_𝛼𝜆_5𝛼,𝜆̇_5𝛼=−𝜕𝐿_𝛼𝜕𝑧_𝛼=𝑎_𝛼𝜆_4𝛼,𝜆̇_6𝛼=−𝜕𝐿_𝛼𝜕𝑧_𝛼=𝑎_𝛼𝜆_3𝛼, (21) where 𝜆_𝑗𝛼(𝑡_𝑓)=0, 𝑗=1,…,6, are the transversality conditions. The Lagrangian is minimized with respect to 𝑢_𝛼 and 𝑢_𝛼 at the optimal pair (𝑢^∗_𝛼,𝑢^∗_𝛼). So the partial derivatives of the Lagrangian with respect to 𝑢_𝛼 and 𝑢_𝛼 are zero. Since, 𝜕𝐿/𝜕𝑢_𝛼=2𝑢_𝛼−𝜆_4𝛼𝑣_𝛼−𝜔₁₁+𝜔₁₂=0, we have 𝑢_𝛼=0.5(𝑣_𝛼𝜆_4𝛼+𝜔₁₁−𝜔₁₂). To determine an explicit expression for the optimal control without 𝜔₁₁ and 𝜔₁₂, we consider the following three cases. (i) On the set {𝑡∣𝑢_𝛼^min<𝑢^∗_𝛼(𝑡)<𝑢_𝛼^max}, we set 𝜔₁₁(𝑡)=𝜔₁₂(𝑡)=0; hence, 𝑢_𝛼=0.5𝑣_𝛼𝜆_4𝛼. (ii) On the set {𝑡∣𝑢^∗_𝛼(𝑡)=𝑢_𝛼^max}, we set 𝜔₁₁(𝑡)=0; hence, 𝑢_𝛼=𝑢_𝛼^max=0.5(𝑣_𝛼𝜆_4𝛼−𝜔₁₂) which implies that 0.5𝑣_𝛼𝜆_4𝛼≥𝑢_𝛼^max. (iii) On the set {𝑡∣𝑢^∗_𝛼(𝑡)=𝑢_𝛼^min}, we set 𝜔₁₂(𝑡)=0; hence, 𝑢_𝛼=𝑢_𝛼^min=0.5(𝑣_𝛼𝜆_4𝛼+𝜔₁₁) which implies that 0.5𝑣_𝛼𝜆_4𝛼≤𝑢_𝛼^min. Combining all the three cases in compact form gives𝑢_𝛼(𝑡)=max𝑢_𝛼^min,min0.5𝑣_𝛼(𝑡)𝜆_4𝛼(𝑡),𝑢_𝛼^max. (22) Using similar arguments, we also obtain the following expression for the second optimal control function:𝑢_𝛼(𝑡)=max𝑢_𝛼^min,min0.5𝑣_𝛼(𝑡)𝜆_3𝛼(𝑡),𝑢_𝛼^max. (23) We point out that the optimality system consists of the state system (19) with the initial conditions, the adjoint or costate system (21) with the terminal conditions, together with the expressions (22) and (23) for the control functions. We show the optimal controls and corresponding states and costates satisfying (19) and (21) by a star superscript ∗. Obviously, if we assume that the set [𝜆^∗_4𝛼(𝑡),𝜆^∗_3𝛼(𝑡)]is a valid α-level set of a positive fuzzy number valued function, say Λ^∗₂(𝑡), then, from (22) and (23),the optimal fuzzy control function 𝑢^∗∈𝐶(𝑇) can be written as𝑢^∗(𝑡)=max𝑢^min,min0.5⊗𝑣^∗(𝑡)⊗Λ^∗₂(𝑡),𝑢^max. (24) Moreover, if we assume that the sets [𝜆^∗_2𝛼(𝑡),𝜆^∗_1𝛼(𝑡)]and [𝜆^∗_6𝛼(𝑡),𝜆^∗_5𝛼(𝑡)]are valid α-level sets of fuzzy number valued functions, say Λ^∗₁(𝑡) and Λ^∗₃(𝑡), respectively, and the sets [𝜆̇^∗_1𝛼(𝑡),𝜆̇^∗_2𝛼(𝑡)], [𝜆̇^∗_3𝛼(𝑡),𝜆̇^∗_4𝛼(𝑡)], and [𝜆̇^∗_5𝛼(𝑡),𝜆̇^∗_6𝛼(𝑡)]are valid α-level sets, then it is easy to see that (ΘΛ^∗₂)⊗_𝑍𝑎 and Λ^∗₃⊗_𝑍ℎ∈𝐹(𝑅). Moreover, we can verify that the fuzzy functions Λ^∗₁(𝑡), Λ^∗₂(𝑡), Λ^∗₃(𝑡), and 𝑢^∗ satisfy the following system of FDEs, with the terminal conditions Λ_𝑗(𝑡_𝑓)=(0,0), 𝑗=1,2,3, and using the derivative in the second form (ii) (see Definition 9): Λ̇₁=ΘΘΛ₁⊗𝜎⊕Λ₃⊗_𝑍ℎ,Λ̇₂=Θ𝑤⊕ΘΛ₁⊕Λ⊗𝑐₂̃⊗𝑘Θ_𝐻𝑢⊕ΘΛ₃⊗𝜏,Λ̇₃=ΘΘΛ₂⊗_𝑍𝑎. (25) Therefore, if we insure that the assumptions mentioned above are satisfied, finding the optimal fuzzy control function can be equivalent to solving a two-point fuzzy boundary value problem (FBVP) which consists of fuzzy system (18) with the initial conditions and the fuzzy system (25) with the final conditions, together with the expression 𝑢(𝑡)=max(𝑢^min,min(0.5⊗𝑣(𝑡)⊗Λ₂(𝑡),𝑢^max)) for the control function. Here, we verify that Λ^∗₁(𝑡), Λ^∗₂(𝑡), Λ^∗₃(𝑡), and 𝑢^∗ satisfy the second equation in (25) and the rest can be verified similarly. Since the functions 𝜆^∗_3𝛼 and 𝜆^∗_4𝛼 are positive from the two last equations in (21), we conclude that 𝜆^∗_5𝛼 and 𝜆^∗_6𝛼 are negative functions. Moreover, from the two first equations in (21), we have𝜆^∗_1𝛼(𝑡)=^𝑡𝑓_𝑡ℎ_𝛼cosh𝜎_𝛼𝜎_𝛼(𝜏−𝑡)𝜆^∗_5𝛼(𝜏)+𝜎_𝛼𝜎_𝛼ℎ_𝛼sinh𝜎_𝛼𝜎_𝛼(𝜏−𝑡)𝜆^∗_6𝛼(𝜏𝜆)𝑑𝜏,_∗2𝛼(𝑡)=^𝑡𝑓_𝑡𝜎_𝛼𝜎_𝛼ℎ_𝛼sinh𝜎_𝛼𝜎_𝛼(𝜏−𝑡)𝜆^∗_5𝛼(𝜏)+ℎ_𝛼cosh𝜎_𝛼𝜎_𝛼(𝜏−𝑡)𝜆^∗_6𝛼(𝜏)𝑑𝜏, (26) which implies that the functions 𝜆^∗_1𝛼 and 𝜆^∗_2𝛼 are negative; hence,Θ𝑤⊕ΘΛ^∗₁⊗𝑐̃⊕Λ^∗₂⊗𝑘Θ_𝐻𝑢^∗⊕ΘΛ^∗₃⊗𝜏_𝛼=−𝑤⊕ΘΛ^∗₁⊗𝑐̃⊕Λ^∗₂⊗𝑘Θ_𝐻𝑢^∗⊕ΘΛ^∗₃⊗𝜏_𝛼=−𝑤_𝛼−𝜆^∗_2𝛼𝑐_𝛼+𝜆^∗_3𝛼𝑘_𝛼−𝑢^∗_𝛼−𝜆^∗_6𝛼𝜏_𝛼=𝜆̇^∗_3𝛼. (27) We summarize our results in the following theorem.

Theorem 13. Assume that Λ^∗₂(𝑡) is a positive fuzzy function on T. A candidate for fuzzy optimal control is given by (24), if the fuzzy functions 𝑢^∗(𝑡), 𝑥^∗(𝑡), 𝑧^∗(𝑡), and 𝑣^∗(𝑡) satisfy (18) using the derivative in the first form (i) and 𝑢^∗(𝑡), Λ^∗₁(𝑡), Λ^∗₃(𝑡), and Λ^∗₂(𝑡) satisfy (25) using the derivative in the second form (ii).

We have solved the optimality system corresponding to patient 𝑊 with the fuzzy weight 𝑤=(0.24×10⁻⁵,0.24×10⁻⁶) and the fuzzy bounds 𝑢^min=(0,0,0.5×10⁻³) and 𝑢^max=(0.15×10⁻²,0.25×10⁻³,0), during the time interval [0, 1800] by using the gradient method [24]. Numerical results show that the optimal adjoint variables satisfy the assumptions mentioned above. Figure 11 shows the optimal fuzzy adjoint variables, while the optimal fuzzy control function 𝑢^∗ and its WCOG as a real-valued output 𝑢^∗ indicating an optimal drug regimen are depicted in Figure 12. Moreover, the optimal fuzzy states 𝑥^∗(𝑡), 𝑧^∗(𝑡), and 𝑣^∗(𝑡) indicating the immune cells level and the HIV viral load in presence of treatment are shown in Figure 13. The values of 𝑢^∗ on the 425th and the 1750th days are, respectively, shown in Figures 14 and 15, for the sake of clarity.

(a)

(b)

(c)

Figure 11 

Optimal fuzzy adjoint variables Λ∗1(a), Λ∗2(b), and Λ∗3(c) versus time.

(a)

(b)

Figure 12 

Optimal fuzzy control (a) and its WCOG (b) versus time.

(a)

(b)

(c)

Figure 13 

CD4+ T-cells level (a), the viral load (b), and CTLs level (c) in presence of treatment versus time in patient W.

Figure 14 

𝑢max(- -), 0.5⊗𝑣∗(425)⊗Λ∗2(425) (–·–), 𝑢∗(425)=min(0.5⊗𝑣∗(425)⊗Λ∗2(425),𝑢max) (—).

Figure 15 

𝑢min (- -), 0.5⊗𝑣∗(1750)⊗Λ∗2(1750) (–·–), 𝑢∗(1750)=max(0.5⊗𝑣∗(1750)⊗Λ∗2(1750),𝑢min) (—).

The proposed treatment by Figure 12, which its intensity is decreasing during the time interval [0, 1800], reduces the proliferation rate of viruses considerably. Therefore, the HIV viral load is very low as shown in Figure 13(b). A very low viral load slows the destruction of CD4+ T-cells which are due to contact with virus particles; hence, a CD4+ T-cells level below 70% is not possible in a treated patient as shown in Figure 13(a). From Figure 13(c), a high level of CD4+ T-cells and a low viral load lead to establishment of a lasting CTL-mediated immune response.

5. Conclusion

In this paper, we proposed a fuzzy mathematical model of HIV dynamic. Simulation results show that the proposed three-dimensional FDEs can describe the ambiguous immune cells level and the HIV viral load which are due to existing patients with various strength of their immune system. Moreover, we utilized the proposed fuzzy model and studied a fuzzy optimal control problem minimizing both the viral load and drug costs. Using the Pontryagin’s Maximum Principle leads us to the conclusion that the fuzzy optimal control function may not exist in general, but an optimality system containing fuzzy state and fuzzy adjoin equations is derived under certain assumptions. Motivated by these results, we tend to exploit necessary and sufficient conditions to the existence of fuzzy solutions for linear fuzzy optimal control problems. We expect to address these problems in further works.

Acknowledgment

The authors would like to thank the anonymous referees for their valuable comments that improved the quality of this paper.