# An Efficient Therapy Strategy under a Novel HIV Model

Authors: Chunming Zhang; Xiaofan Yang; Wanping Liu; Lu-Xing Yang
Journal: Discrete Dynamics in Nature and Society
 (2011)

---

## Abstract

By incorporating the chemotherapy into a previous model describing the interaction of the immune system with the human immunodeficiency virus (HIV), this paper proposes a novel HIV virus spread model with control variables. Our goal is to maximize the number of healthy cells and, meanwhile, to minimize the cost of chemotherapy. In this context, the existence of an optimal control is proved. Experimental results show that, under this model, the spread of HIV virus can be controlled effectively.


---

## Body

## 1. Introduction

Numerous studies have been devoted to the description and understanding of the spread of infectious diseases (especially, the acquired immunodeficiency syndrome (AIDS)) [1–18]. Mathematical modeling of the human immunodeficiency virus (HIV) viral dynamics has offered many insights into the pathogenesis and treatment of HIV [1, 2, 4–10, 12–16, 18]. Consequently, many mathematical models have been developed to depict the relationships among HIV, etiological agent for AIDS and CD4+T lymphoblasts, which are the targets for the virus [13]. Some of these models investigate how to avoid an excessive use of drugs because it might be toxic to human body and, hence, cause damages [1, 4–6, 8–11, 14, 15, 17, 18].Recently, Sedaghat et al. [13] proposed a model, which describes the law governing the transition of two populations of target cells, the T cells (the abbreviation of the CD4+T lymphoblasts) and the M cells (say, macrophages, T cells in a lower state of activation, or another cell type), in the effect of free virus (see Figure 1). The T cells produce most of the plasma virus and are responsible for the first-phase decay, while the M cells are responsible for the second-phase decay. T cells are classified into three categories: TU cells (uninfected T cells), T1 cells (early-stage infected T cells), and T2 cells (late-stage infected T cells). Let TU, T1 and T2 denote the numbers of TU cells, T1 cells, and T2 cells, respectively. Likewise, M cells are classified into three categories: Mu cells (uninfected M cells), M1 cells (early-stage infected M cells), and M2 cells (late-stage infected M cells). Let Mu, M1, and M2 denote the numbers of MU cells, M1 cells and M2 cells, respectively. Besides, let V denote the number of free viruses. Sedaghat et al. [13] made the following reasonable assumptions.(A1)
TU cells are produced with constant rate θT. MU cells are produced with constant rate θM.(A2)
TU cells become T1 cells with constant rate βT. MU cells become M1 cells with constant rate βM.(A3)
T1 cells become T2 cells with constant rate kT. M1 cells become M2 cells with constant rate kM.(A4)
These cells die with constant ratesδTU, δT1, δT2, δMU, δM1, and δM2 respectively.(A5)
Free viruses (V) are cleared at a rate c, produced by T2 cells with a burst size of NT, and produced by M2 cells with a burst size of NM, respectively.Figure 1
The HIV model.Under these assumptions, Sedaghat et al. [13] deduced the following system of ordinary differential equations:(1.1)dTUdt=θT-δTUTU-βTTUV,dT1dt=βTTUV-(δT1+kT)T1,dT2dt=kTT1-δT2T2,dMUdt=θM-δMUMU-βMMUV,dM1dt=βMMUV-(δM1+kM)M1,dM2dt=kMM1-δM2M2,dVdt=NTT2+NMM2-cV.
For a highly simplified version of this system, Sedaghat et al. [13] derived its analytic solution.It is well known [5, 6, 8–11, 13, 15, 17] that there are mainly two categories of anti-HIV drugs: the reverse transcriptase inhibitors (RTIs), which prevent new HIV infection by disrupting the conversion of viral RNA into DNA inside of T cells, and the protease inhibitors (PIs), which reduce the number of virus particles produced by actively-infected T cells.In consideration of this, this paper introduces a novel HIV model by incorporating the drug dosage into the above-mentioned model. Our goal is to maximize the number of healthy cells and, meanwhile, to minimize the cost of chemotherapy. In this context, the existence of an optimal control strategy is proved. Experimental results show that, under this model, the spread of HIV virus can be controlled effectively.


## 2. Presentation of a New Model

For our purpose, let us introduce the following notations (see Figure2): u1(t):
the dosage of RTI at timet, which is assumed to take values in the interval [0,1];u2(t):
the dosage of PI at timet, which is assumed to take values in [0,1];γ1:
the capability of preventingTU cells from becoming T1 cells with per unit dosage of RTI;γ2:
the capability of preventingMU cells from becoming M1 cells with per unit dosage of RTI;α1:
the capability of preventingT2 cells from producing viruses with per unit dosage of PI;α2:
the capability of preventingM2 cells from producing viruses with per unit dosage of PI.Figure 2
The HIV model with therapy strategy.Next, let us consider the following assumptions.(A6)
Due to the effect of RTIs,TU cells become T1 cells with rate βT(1-u1(t))γ1, and MU cells become M1 cells with rate βM(1-u1(t))γ2, where γ1 and γ2 are constants.(A7)
Due to the effect of PIs, Free viruses (V) are produced by T2 and M2 cells with a burst size of α1(1-u2(t))NT and α2(1-u2(t))NM, respectively, where α1 and α2 are constants.Under assumptions (A1)–(A7), we can derive the following system of ordinary differential equations:(2.1)dTUdt=θT-δTUTU-βTVTU(1-u1)γ1,dT1dt=βTVTU(1-u1)γ1-(δT1+kT)T1,dT2dt=kTT1-δT2T2,dMUdt=θM-δMUMU-βMVMU(1-u1)γ2,dM1dt=βMVMU(1-u1)γ2-(δM1+kM)M1,dM2dt=kMM1-δM2M2,dVdt=α1NTT2(1-u2)+α2NMM2(1-u2)-cV.Our target is to maximize the objective functional by increasing the number of healthy T and M cells and minimizing the cost based on the percentage effect of the chemotherapy given. For that purpose, we introduce the followingobjective functional(2.2)J(u1(t),u2(t))=∫t0t1{B1TU+B2MU-[A1u12+A2u22]}dt,
where B1,B2 represent the benefit of per TU cell and per MU cell, respectively, and A1,A2 represent the cost of per unit RTI and per unit PI, respectively. Our goal is to obtain an optimal control pair (u1*,u2*) such that(2.3)J(u1*,u2*)=max⁡{J(u1,u2):(u1,u2)∈U},
where 𝒰 is the admissible control set defined by(2.4)U=U1×U2,U1=U2={u(t):umeasurable,0≤u(t)≤1,t∈[t0,t1]}.


## 3. Existence of an Optimal Control Pair

For our purpose, let us introduce the following four assumptions.(A8)
The set of control and corresponding state variables is nonempty.(A9)
The admissible control set𝒰 is closed and convex.(A10)
All the right hand sides of equations of system (2.1) are continuous, bounded above by a sum of bounded control and state, and can be written as a linear function of u with coefficients depending on time and state.(A11)
There exist positive constantsc1,c2 and β>1 such that the integrand (denoted by L(y,u,t)) of the objective functional (2.2) is concave and satisfies the condition L(y,u,t)≤c1-c2(u12+u22)β/2.In what follows, it is always assumed that assumptions (A1)–(A7) hold.Theorem 3.1.
Consider system (2.1) with initial conditions, and the objective functional (2.2). There exists u*=(u1*,u2*) such that
(3.1)J(u1*,u2*)=max⁡u∈UJ(u1,u2).Proof.
It suffices to verify the assumptions (A8)–(A11) with respect to the seven ODEs of system (2.1).
Since the coefficients involved in the system are bounded, and each state variable of the system is bounded on the finite time interval, it follows by a result (see AppendixA) from [19] we can obtain the existence to the solution of the system (2.1).
The control set𝒰=U1×U2 is obviously closed and convex, because both U1 and U2 are closed and convex sets.
By definition, each right hand side of the ODEs of system (2.1) is continuous and can be written as a linear function of u with coefficients depending on time and states. The fact that all state variables TU, T1, T2, MU, M1, M2, V, and 𝒰 are bounded on [t0,t1], implies the rest of assumption (A10).
It is easy to see thatL(y,u,t) is concave in 𝒰. By setting c1=max⁡{B1TU+B2MU}, c2=inf⁡(A1,A2) and β=2, we can derive
(3.2)L(y,u,t)=B1TU+B2MU-[A1u12+A2u22]≤c1-c2(u12+u22).
The proof is complete.


## 4. Optimally Controlling Chemotherapy

In this section, we discuss the theorem related to the characterization of the optimal control. This result depends on thePontryagin’s Maximum Principle, which gives necessary conditions for the optimal control. First, we rewrite the system (2.1) in the following vector notation: (4.1)dy(t)dt=A(y,u,t);∀t>t0,∀u∈U,y(t0)=y0,
where y(t) and A(y,u,t) are given by (4.2)y(t)=(TU(t),T1(t),T2(t),MU(t),M1(t),M2(t),V(t))T,A(y,u,t)=(g1(y,u,t),g2(y,u,t),…,g6(y,u,t),g7(y,u,t))T.TheHamiltonian associated with our problem is (4.3)H(y,u,p,t)=L(y,u,t)+λ(t)TA(y,u,t),
where the adjoint vector λ(t) is defined by the adjoint equation (4.4)dλ(t)dt=-Ayλ(t)-Ly,λ(t1)=0.Here(4.5)Ay=(∂g1∂TU∂g2∂TU∂g3∂TU∂g4∂TU∂g5∂TU∂g6∂TU∂g7∂TU∂g1∂T1∂g2∂T1∂g3∂T1∂g4∂T1∂g5∂T1∂g6∂T1∂g7∂T1∂g1∂T2∂g2∂T2∂g3∂T2∂g4∂T2∂g5∂T2∂g6∂T2∂g7∂T2∂g1∂MU∂g2∂MU∂g3∂MU∂g4∂MU∂g5∂MU∂g6∂MU∂g7∂MU∂g1∂M1∂g2∂M1∂g3∂M1∂g4∂M1∂g5∂M1∂g6∂M1∂g7∂M1∂g1∂M2∂g2∂M2∂g3∂M2∂g4∂M2∂g5∂M2∂g6∂M2∂g7∂M2∂g1∂V∂g2∂V∂g3∂V∂g4∂V∂g5∂V∂g6∂V∂g7∂V)=(E,F),where(4.6)E=(-δTU-βTV(1-u1)γ1βTV(1-u1)γ100-(δT1+kT)kT00-δT2000000000000),F=(00000000000α1NT(1-u2)-δMU-βMV(1-u1)γ2βMV(1-u1)γ2000-(δM1+kM)kM000-δM2α2NM(1-u2)000-c).
In addition, the Ly in system (4.3) is (4.7)Ly=(∂L∂TU,∂L∂T1,∂L∂T2,∂L∂MU,∂L∂M1,∂L∂M2,∂L∂V)T,=(B1,0,0,B2,0,0,0)T.Next, adding the penalty term will give us the optimality condition(4.8)ξ(y,u,λ,t)=H(y,u,λ,t)+Γ(u(t))ω(t),
where Γ is an operator from ℝ2 to ℝ4 defined by (4.9)Γ(u(t))=(1-u1(t),u1(t),1-u2(t),u2(t)),ω(t)=(ω11(t)ω12(t)ω21(t)ω22(t)),
where all ωij,i,j=1,2 are nonnegative penalty multipliers satisfying the following conditions: (4.10)(1-u1*(t))ω11(t)=u1*(t)ω12(t)=(1-u2*(t))ω21(t)=u2*(t)ω22(t)=0.According to thePontryagin’s Maximum Principle, if the control u*(t) and the corresponding state y*(t) constitute an optimal pair, there exists an adjoint vector λ(t) defined system (4.4) such that the function ξ(y,u,λ,t) defined by (4.8) reaches its maximum on the set 𝒰 at the point u*. This gives the following result.Theorem 4.1.
Given an optimal control pairu*(t)=(u1*(t),u2*(t)) and a solution y*(t)=(TU*(t),T1*(t),T2*(t),MU*(t),M1*(t),M2*(t),V*(t)) of the corresponding system, then there exist seven adjoint variables λ1(t),λ2(t),…,λ7(t) satisfying
(4.11)dλ1dt=[δTU+βTV(1-u1)γ1]λ1-βTV(1-u1)γ1λ2-B1,dλ2dt=(δT1+kT)λ2-kTλ3,dλ3dt=δT2λ3-α1NT(1-u2)λ7,dλ4dt=[δMU+βMV(1-u1)γ2]λ4-βMV(1-u1)γ2λ5-B2,dλ5dt=(δM1+kM)λ5-kMλ6,dλ6dt=δM2λ6-α2NM(1-u2)λ7,dλ7dt=cλ7+βTTU(1-u1)γ1λ1-βTTU(1-u1)γ1λ2+βMMU(1-u1)γ2λ4-βMMU(1-u1)γ2λ5,
with the final conditions
(4.12)λ1(t1)=λ2(t1)=⋯=λ7(t1)=0.
Furthermore, u1*(t)=min⁡{max⁡{0,R1(t)},1},u2*(t)=min⁡{max⁡{0,R2(t)},1}, where
(4.13)R1(t)=V*2A1(βTTU*γ1(λ1-λ2)+βMMU*γ2(λ4-λ5)),R2(t)=-λ72A2(α1NTT2+α2NMM2).Proof.
According to the previous section, an optimal couple(y*(t),u*(t)) exists for maximizing the objective functional (2.2) subject to the system (2.1). Therefore, by Pontryagin’s Maximum Principle, there exists a vector λ(t)=(λ1(t),…,λ7(t))T satisfying
(4.14)λ(t)dt=-∂H∂y=-Ly-Ayλ(t).
That yields
(4.15)λ1(t)dt=-(∂g1(y*,u*,t)∂TU,…,∂g7(y*,u*,t)∂TU)λ(t)-∂L(y*,u*,t)∂TU,λ2(t)dt=-(∂g1(y*,u*,t)∂T1,…,∂g7(y*,u*,t)∂T1)λ(t)-∂L(y*,u*,t)∂T1,λ3(t)dt=-(∂g1(y*,u*,t)∂T2,…,∂g7(y*,u*,t)∂T2)λ(t)-∂L(y*,u*,t)∂T2,λ4(t)dt=-(∂g1(y*,u*,t)∂MU,…,∂g7(y*,u*,t)∂MU)λ(t)-∂L(y*,u*,t)∂MU,λ5(t)dt=-(∂g1(y*,u*,t)∂M1,…,∂g7(y*,u*,t)∂M1)λ(t)-∂L(y*,u*,t)∂M2,λ6(t)dt=-(∂g1(y*,u*,t)∂M2,…,∂g7(y*,u*,t)∂M2)λ(t)-∂L(y*,u*,t)∂M2,λ7(t)dt=-(∂g1(y*,u*,t)∂V,…,∂g7(y*,u*,t)∂V)λ(t)-∂L(y*,u*,t)∂V.
Through simple calculations, we derive system (4.11). The Pontryagin’s Maximum Principle gives the following necessary conditions to obtain the optimal pair (y*,u*):
(4.16)∂ξ(y*,u*,λ,t)∂u1=0,∂ξ(y*,u*,λ,t)∂u2=0,
where ξ(y*,u*,λ,t)=H(y*,u*,λ,t)+Γ(u*(t))ω(t). From (4.10) and (4.16), we have
(4.17)∂ξ(y*,u*,λ,t)∂u1=0⟹∂L(y*,u*,t)∂u1+∂(λ(t)TA(y*,u*,t))∂u1+∂(Γ(u*(t))ω(t))∂u1=0,
which implies
(4.18)u1*(t)=V*2A1(βTTU*γ1(λ1-λ2)+βMMU*γ2(λ4-λ5)).
On the other hand,
(4.19)∂ξ(y,u*,λ,t)∂u2=0⟹∂L(y*,u*,t)∂u2+∂(λ(t)TA(y,u*,t))∂u2+∂(Γ(u*(t))ω(t))∂u2=0,
which indicates
(4.20)u2*(t)=-λ72A2(α1NTT2+α2NMM2).
Now from the constraint condition, the following three cases arise.Case 1.
t∈{t:0<u1*(t)<1} and ω11(t)=ω12(t)=0. Then u1*(t)=R1(t).Case 2.
t∈{t:u1*(t)=0} and ω11(t)=0. Then 0=u1*(t)=R1(t)+ω12(t), which implies u1*(t)≥R1(t) because ω12(t)≥0.Case 3.
t∈{t:u1*(t)=1} and ω12(t)=0. Then u1*(t)=R1(t)-ω11(t), which leads to 1=u1*(t)≤R1(t), owing to ω11≥0.Hence, we have u1*(t)=min⁡{max⁡{0,R1(t)},1}. Similarly, we can get that u2*(t)=min⁡{max⁡{0,R2(t)},1}.
The proof is complete.Now, the optimality system is given by incorporating the optimal control pair in the state system coupled with the adjoint system. Thus, we have(4.21)dy*(t)dt=A(y*,u*,t);∀t>t0,dλ(t)dt=-Ay*λ(t)-Ly*,y*(t0)=y0*,λ(t1)=0.
We substitute the expressions u*=(u1*,u2*) in the above system. The uniqueness of the solution of the optimality system can be derived by a standard method (refer to [6] for more details on the proof).


## 5. Numerical Algorithm and Results

The resolution of the optimal system is created improving the Gauss-Seidel-like implicit finite-difference method developed by [7] and denoted by GSS1 method. It consists on discretizing the interval [t0,t1] at the points tk=kl+t0(k=0,1,…,n), where l is the time step.In the following, we define the state and adjoint variablesTU(t), T1(t), T2(t), MU(t), M1(t), M2(t), V(t), λ1(t)~λ7(t) and the controls u1(t) and u2(t) in terms of nodal points TUk, T1k, T2k, MUk, M1k, M2k, Vk, λ1k~λ7k, u1k, u2k as the state and adjoint variables and the controls at initial time t0, while TUn, T1n, T2n, MUn, M1n, M2n, Vn, λ1n~λ7n, u1n, u2n as the state and adjoint variables and the controls at final time t1. As it is well known that the approximation of the time derivative by its first-order forward-difference is given, for the first state variable TU, by (5.1)dTU(t)dt=lim⁡l→0TU(t+l)-TU(t)l.We use the scheme developed by Gumel et al. [7] in the following way: (5.2)TUk+1-TUkl=θT-δTUTUk+1-βTγ1Vk(1-u1k)TUk+1.
Analogously, we have (5.3)T1k+1-T1kl=βTγ1Vk(1-u1k)TUk+1-(δT1+kT)T1k+1,T2k+1-T2kl=kTT1k+1-δT2T2k+1,MUk+1-MUkl=θM-δMUMUk+1-βMγ2Vk(1-u1k)MUk+1,M1k+1-M1kl=βMγ2Vk(1-u1k)MUk+1-(δM1+kM)M1k+1,M2k+1-M2kl=kMM1k+1-δM2M2k+1,Vk+1-Vkl=α1NTT2k+1(1-u2k)+α2NMM2k+1(1-u2k)-cVk+1.By applying an analogous technology, we approximate the time derivative of the adjoint variables by their first-order backward-difference and we use the appropriated scheme as follows:(5.4)λ1n-k-λ1n-k-1l=[δTU+βTVk+1(1-u1k)γ1]λ1n-k-1-βTVk+1(1-u1k)γ1λ2n-k-B1,λ2n-k-λ2n-k-1l=(δT1+kT)λ2n-k-1-kTλ3n-k,λ3n-k-λ3n-k-1l=δT2λ3n-k-1-α1NT(1-u2k)λ7n-k,λ4n-k-λ4n-k-1l=[δMU+βMVk+1(1-u1k)γ2]λ4n-k-1-βMVk+1(1-u1k)γ2λ5n-k-B2,λ5n-k-λ5n-k-1l=(δM1+kM)λ5n-k-1-kMλ6n-k,λ6n-k-λ6n-k-1l=δM2λ6n-k-1-α2NM(1-u2k)λ7n-k,λ7n-k-λ7n-k-1l=cλ7n-k-1+βTTUk+1(1-u1k)γ1(λ1n-k-1-λ2n-k-1)+βMMUk+1(1-u1k)γ2(λ4n-k-1-λ5n-k-1).Hence, we can establish an algorithm to solve the optimality system and then to compute the optimal control pair by employing the GSS1 method (5.2)–(5.4) that we denote by IGSS1 method here (see Appendix B).


### 5.1. Numerical Results

By making some parameter value choices, computer simulation experiments are done to verify the effectiveness of our new model by comparing the disease progression before and after introducing the two optimal control variablesu1*(t), u2*(t). For the following parameters and initial values:θT=10, θM=10, δTU=0.02, δT1=0.5, δT2=1, δMU=0.0495, δM1=0.0495, δM2=0.0495, βT=0.00008, βM=0.00008, kT=0.1, kM=0.1, NT=100, NM=100, TU0=1000, T10=0, T20=0, MU0=1000, M10=0, M20=0, V0=1, c=0.03.The experimental results obtained are listed in Table1 (in which “before treatment’’ and “after treatment’’ are denoted by BT and AT, resp.).Table 1Time (days)TU BTTU ATMU BTMU ATV BTV AT01000100010001000112980.2369980.2012926.3417926.149415.3344419.5681364955.4799953.1210856.9665852.1262283.2647334.658886839.4257889.1900754.0975754.10064710.942315.448098270.1388833.1982372.5725671.291138782.17297.340081011.92899784.045341.26735601.137199155.03280.27154200.511601628.53181.021524388.0657255003.3208.54835300.379269597.73720.756860313.0676335762.8155.17956400.334284574.20370.666976268.0361376288.7115.46816500.322352555.23090.643089240.5969387774.585.919158For more clearness, it is better to present these comparative results by the following graphs. When the viruses attack the human body, uninfected T and M cells decrease (see Figures3 and 4).Figure 3
Uninfected T cells.Figure 4
Uninfected M cells.The viruses do not stop to proliferate and so its abundance dramatically increases (see Figure5). However, after introducing the optimal controls, the situation changes. A few days later, the effect of chemotherapy starts to appear; which explains the growth of uninfected T and M cells and the diminishing of viruses (see Figure 6).Figure 5
Virus population before optimal controls.Figure 6
Virus population after optimal controls.Finally, the optimal controlsu1*(t), u2*(t) for drug administration are presented through Figures 7 and 8.Figure 7
Optimal control variableu1*(t).Figure 8
Optimal control variableu2*(t).


## 5.1. Numerical Results

By making some parameter value choices, computer simulation experiments are done to verify the effectiveness of our new model by comparing the disease progression before and after introducing the two optimal control variablesu1*(t), u2*(t). For the following parameters and initial values:θT=10, θM=10, δTU=0.02, δT1=0.5, δT2=1, δMU=0.0495, δM1=0.0495, δM2=0.0495, βT=0.00008, βM=0.00008, kT=0.1, kM=0.1, NT=100, NM=100, TU0=1000, T10=0, T20=0, MU0=1000, M10=0, M20=0, V0=1, c=0.03.The experimental results obtained are listed in Table1 (in which “before treatment’’ and “after treatment’’ are denoted by BT and AT, resp.).Table 1Time (days)TU BTTU ATMU BTMU ATV BTV AT01000100010001000112980.2369980.2012926.3417926.149415.3344419.5681364955.4799953.1210856.9665852.1262283.2647334.658886839.4257889.1900754.0975754.10064710.942315.448098270.1388833.1982372.5725671.291138782.17297.340081011.92899784.045341.26735601.137199155.03280.27154200.511601628.53181.021524388.0657255003.3208.54835300.379269597.73720.756860313.0676335762.8155.17956400.334284574.20370.666976268.0361376288.7115.46816500.322352555.23090.643089240.5969387774.585.919158For more clearness, it is better to present these comparative results by the following graphs. When the viruses attack the human body, uninfected T and M cells decrease (see Figures3 and 4).Figure 3
Uninfected T cells.Figure 4
Uninfected M cells.The viruses do not stop to proliferate and so its abundance dramatically increases (see Figure5). However, after introducing the optimal controls, the situation changes. A few days later, the effect of chemotherapy starts to appear; which explains the growth of uninfected T and M cells and the diminishing of viruses (see Figure 6).Figure 5
Virus population before optimal controls.Figure 6
Virus population after optimal controls.Finally, the optimal controlsu1*(t), u2*(t) for drug administration are presented through Figures 7 and 8.Figure 7
Optimal control variableu1*(t).Figure 8
Optimal control variableu2*(t).


## 6. Conclusions

By incorporating the chemotherapy into a previous model describing the interaction of the immune system with the human immunodeficiency virus (HIV), this paper has proposed a novel HIV virus spread model with control variables. Our goal is to maximize the number of healthy cells and, meanwhile, to minimize the cost of chemotherapy. In this context, the existence of an optimal control has been proved. Experimental results show that, under this model, the spread of HIV virus can be controlled effectively.Our next work is to study other kinds of models, especially those with impulsive drug effect.



---

*Source: 828509-2011-10-24.xml*