Rachid Lounes

Hector De Arazoza

Ying-Hen Hsieh
email: hsieh@mail.cmu.edu.tw

IN CUBA,1986-2008: MODELING CHANGES IN DETECTION WITH TEMPORALLY VARYING PARAMETERS June 6, 2012

Keywords: HIV/AIDS epidemic, Cuba, contact tracing, detection, compartmental model

Background HIV/AIDS prevalence in Cuba has remained consistently the lowest in the Caribbean region. The Cuban HIV/AIDS program features a detection system which has experienced many changes over the years, including the gradual discontinuation of sanatorial system in the 1990's and increased detections through the family doctors after 2000.

Methodology/Major

Findings A nonlinear compartmental model with several types of detections is proposed. Analysis reveals that local asymptotic stability for disease-free equilibrium (DFE) can be achieved if: (i) random screening is sufficiently effective; (ii) infection by detected HIV-positive individuals is minimal. If the average number of infections by a known infective exceeds unity, the endemic equilibrium is always unstable and the total number of infectives could increase without bound, provided that the initial infective population sizes are sufficiently large. On the other hand, if the number of infections by a known infective is less than one, then either the DFE or the endemic equilibrium is globally asymptotically stable, leading to a more manageable epidemic, even if the disease is not eradicated. We make use of the Cuban HIV/AIDS data to fit the model during two separate time periods in 1986-2008 to reflect the implementation of different types of detections. The reproduction numbers for each time period are then computed from the two sets of estimated parameter values.

Conclusions

The results indicate that the HIV/AIDS epidemic in Cuba is currently near an endemic equilibrium. Random screening is still the most important mean of surveillance and detection, while infections by the known infectives must be monitored closely and kept in check, thus highlighting the importance of education and behavior change. The combined use of theoretical analysis and data-oriented application in this work provides new insight into the dynamics of a sexually transmitted disease and its epidemiology.

Introduction

The Cuban HIV/AIDS program was first initiated in 1983, although the first HIV case in Cuba, a heterosexual male returning from abroad, did not occur until December 1985, while the first AIDS case was diagnosed in Cuba in April 1986 and passed away later that month. HIV/AIDS prevalence in Cuba has consistently remained the lowest in the Caribbean region, which has the second highest HIV prevalence rate in the world after Sub-Saharan Africa [START_REF] Inciardi | HIV/AIDS in the Caribbean Basin[END_REF]. The UNAIDS Epidemiological Fact Sheet on HIV and AIDS for Cuba (2008 update) reports an HIV prevalence of less than 0.1 % for adults. With 11.2 million inhabitants and 6.1 million of age 15-49. Life expectancy in Cuba is 75 years and 80 years for men and women, respectively [START_REF] Who | World Health Organisation: Cuba[END_REF], the highest in Latin America.

The Cuban HIV/AIDS program features a detection system [START_REF] De Arazoza | The HIV/AIDS epidemic in Cuba: Description and tentative explanation of its low HIV prevalence[END_REF] that allows for the detection of HIV cases from several sources, some of which were initiated at the beginning of the program while others were introduced later. Some have been discontinued in time for various reasons. The earliest and most prominent detections came from random screening of blood donors, persons with other sexually transmitted infections (STI), persons who are admitted to hospital with suspicion of HIV infection or subject to specific procedures like dialysis, and persons who volunteered to be tested. In addition, there were detections through persons having received a recommendation for HIV testing from his/her general practitioner (family doctor) and through sexual partner tracing. Other minor sources of detection include testing of pregnant women and prison inmates that were implemented for only some periods of time during the years 1986-2008, when more than 34 million tests were performed citeCub4. Since 2002, the total number of tests performed has stabilized to 1.6-1.7 million test every year [START_REF] Minsap | CUBA: Plan Estratégico Nacional ITS/VIH/SIDA 2007-2011[END_REF]. Since 1990, each time a person is tested for HIV, she/he is informed that such a test is going to be performed. From 1986, in keeping with the "partner notification program", a person tested to be HIV-positive is invited to give names and contact details of his/her sexual partners during the past two years. These partners are then traced and a recommendation for voluntary HIV testing is made.

It is important to note that the detection system has experienced changes over the years. The impact of this change in detection, both qualitatively and quantitatively, on the dynamics of the HIV epidemic, and more importantly on the disease spread, is largely unknown.

The primary objective of the Cuban Program to control HIV/AIDS epidemic is active search to detect persons infected with HIV long before they show any signs of AIDS. The focus of our modeling is not to model how new infections by HIV are generated, but how the HIV-infected persons are detected and how it changes over time. In a primarily numerical investigation [START_REF] De Arazoza | Modeling Detection of HIV in Cuba[END_REF], a model was proposed to fit the Cuban data mainly to estimate the basic reproduction numbers for Cuban HIV. In this work we will give a complete analysis of this model and further make use of the analytical results to fully ascertain the temporal changes in the effectiveness of the different types of detection.

Methods

The Model

We will briefly present the model introduced in [START_REF] De Arazoza | Modeling Detection of HIV in Cuba[END_REF] Model Variables:

1. X(t): number of HIV infected persons that do not know they are infected at time t, 2. Y (t): number of HIV infected persons that know they are infected at time t,

Z(t): number of persons with AIDS at time t.

The variable Y (t) is divided into 3 subcategories according to three method of detection: Y 1 (t) denotes random search, Y 2 (t) is contact tracing and Y 3 (t) is detection through the family doctor system. These three subclasses are important for our investigation of the different types of detection with each one being assigned a different detection parameter, which will be used when fitting the model to the Cuban data. The model is given as follows:

dX dt = σX -k 2 XY -k 3 X 2 + λ ′ Y, dY 1 dt = k 1 X -γ ′ Y 1 , dY 2 dt = k 2 X Y -γ ′ Y 2 , ( 1
)
dY 3 dt = k 3 X 2 -γ ′ Y 3 , dZ dt = βX + β ′ Y -µ ′ Z.
where

Y (t) = Y 1 (t) + Y 2 (t) + Y 3 (t) and γ ′ = β ′ + µ and σ = λ -k 1 -β -µ.
To simplify the model equations, we combined the unknown infectives Y (t) in our system. The model dynamics is then described by the following system:

dX dt = σX -k 3 X 2 -k 2 XY + λ ′ Y, dY dt = (k 1 + k 3 X)X -γ ′ Y + k 2 X Y, ( 2
)
dZ dt = βX + β ′ Y -µ ′ Z.
We consider the system only in the region D = {X ≥ 0, Y ≥ 0, Z ≥ 0}. Clearly, D is positively invariant under the flow induced by [START_REF] Who | World Health Organisation: Cuba[END_REF]. The model flow diagram is given in Figure 1.

Unknown HIV-infected X(t) Known HIV-infected Y(t) AIDS Z(t) ȜԞ +Ȝ'Y k 1Ԟ+k2XY +k3X 2 ß' ß µ µ µ' Figure 1: Model flow diagram.
We give the following remarks regarding model (2):

1. In (2) there are three ways for individuals to go from "unknown HIV infected" (X) to "known HIV infected" (Y ). One is through the nonlinear term k 2 XY for contact tracing, where the individual is found through his contacts with persons that are known to live with HIV. The term k 3 X 2 models the detection through family doctors. The third way they can be detected is through the term k 1 X which models all the other "random" ways of searching for seropositives. It is important to note that 1/k 1 can be viewed as the mean time from infection to detection for the persons found through a random screening. Alternatively, we could express all manners in which an unknown infective can be detected to be in the general form of "F (X, Y ) X", where F (X, Y ) is a recruitment function from the class X into the class Y . For our present study we will take the function F (X, Y ) as a polynomial of degree 1:

k 1 + k 2 Y + k 3 X.
2. We assume that the known HIV-infected persons are infectious, but at a much lower rate than those that do not know they are infected due to education or change of behavior. This has been shown, using Cuban HIV data of 1986-2000, that the infection rate of the known HIV-infected persons is about 5% of that of persons who did not yet know that they are infected [START_REF] De Arazoza | A non linear model for a sexually transmitted disease with contact tracing[END_REF]. In this case the infection rate of known HIV-positives λ ′ is taken as a fraction of the infection rate of the undetected HIV-positive persons λ.

3.

We assume that once a person develops AIDS it is no longer infectious.

Random detection is carried out at with a constant effort (in general) it includes the testing of all blood donations at hospitals and blood banks, testing of all pregnant women, testing of all persons that will follow medical procedures etc. . These tests are carried out automatically (in a way) by the Health System. Contact Tracing detections is carried out by the epidemiology department at each municipality throughout the Cuban territory. It requires an effort from the system of interviewing the persons already detected to obtain from them (if they are willing to do it) the information on their sexual contacts in the past two years, then the epidemiology dept. has to find the persons and talk to them and convince them to have a test. This is a time consuming and needs manpower to be effective. Nevertheless it is the more effective source of detection. While other sources of detection get less than 1% (even less than 0.1%) positive cases from the total number of tests performed, Contact Tracing has around 3% of positive detections with a lot less number of tests performed. In terms of efficiency Contact Tracing is the most efficient source of detection, but it is probably the most expensive. The General Practitioners source of detection has become a very important part of the system detecting around a third of new detections in a year. It is a system that relies in the structure of the Cuban Health System and it also takes part in National education campaigns to get the individuals to get themselves tested for HIV. In this way the Family Doctors influence the detection system in more ways than the data reflects.

Preliminary Analysis

We will first consider the following two cases; Case 1. If λ ′ -γ ′ = 0, the system has a disease-free equilibrium P 0 = (0, 0, 0) if σ + k 1 = λ -γ ̸ = 0 and a set of endemic equilibria

P * = (X * , Y * , Z * ) of the form (λ ′ -k 2 X * )Y * + (σ + k 3 X * )X * = 0 and Z * = βX * +β ′ Y * µ ′ if σ + k 1 = λ -γ = 0.
Case 2. If λ ′ -γ ′ ̸ = 0, the system has two equilibria, P 0 = (0, 0, 0) is the disease-free equilibrium and P * = (X * , Y * , Z * ) is the endemic equilibrium, where

X * = A B , Y * = X * × C 1 C 2 , Z * = βX * + β ′ Y * µ ′ . ( 3
)
with

A = σ γ ′ + λ ′ k 1 , C 1 = σ + k 1 = λ -γ, C 2 = γ ′ -λ ′ and B = k 2 C 1 + k 3 C 2 .
Set

R 1 = λ ′ γ ′ , R 2 = λ γ and R 0 = λ/(k 1 γ) + λ ′ k 1 /[γ ′ (k 1 + γ)] = k 1 k 1 + γ R 1 + γ k 1 + γ R 2 . ( 4
)
A, C 1 and C 2 can be rewritten as

A = γ ′ (k 1 + γ)(R 0 -1), C 1 = γ(R 2 -1), , and
C 2 = γ ′ (1 -R 1 ).
The parameters R 0 , R 1 and R 2 are the threshold parameters or reproduction numbers of the model, and play a significant role in the analysis of the behaviour of trajectories for [START_REF] Granich | Cuba's national AIDS program. The first decade[END_REF].

Since Case 1 requires specific values of the parameters and hence is of little practical importance, we shall suppose in what follows that λ ′ -γ ′ ̸ = 0.

The endemic equilibrium is feasible (i.e. P * ∈ D) if and only if

A × B = (R 0 -1) ( γk 2 (R 2 -1) + γ ′ k 3 (1 -R 1 )
) > 0 and ( 5)

C 1 × C 2 = γγ ′ (R 2 -1)(1 -R 1 ) > 0 (6)
More precisely, the endemic equilibrium is feasible if and only if

R 0 -1, R 2 -1 and 1 -R 1 (7)
have the same sign.

The Jacobian matrix of the linear approximation of the system in a neighborhood of an equilibrium point P = (X, Y, Z) is given by

J(X, Y, Z) =   σ -2k 3 X -k 2 Y λ ′ -k 2 X 0 k 1 + 2k 3 X + k 2 Y -γ ′ + k 2 X 0 β β ′ -µ ′   .
There is one eigenvalue -µ ′ that is always strictly negative, which is associated with the variable Z, to study the stability of the equilibria of system (2), we only need to consider the system

dX dt = σX -k 3 X 2 -k 2 XY + λ ′ Y, dY dt = (k 1 + k 3 X)X -γ ′ Y + k 2 X Y. ( 8
)
in the region

D ′ = {(X, Y ) | X ≥ 0, Y ≥ 0} ⊂ D. D ′ ,
which is positively invariant under the flow induced by [START_REF] Granich | Cuba's national AIDS program. The first decade[END_REF]. We will denote J 1 the Jacobian matrix for the system (8) that is formed by the first two rows and columns of the matrix J(X, Y, Z). Moreover, we will also denote Q 0 = (0, 0) and Q * = (X * , Y * ) to be the respective disease-free and endemic equilibria in D ′ .

We will give the proof for the stability of the equilibria.

Simulating the evolution of the epidemic.

We will use the model to simulate the dynamics of the epidemic under several possible changes in the detection system. For this we will take the values we found for the parameters in [START_REF] De Arazoza | Modeling Detection of HIV in Cuba[END_REF] for the period 1999 -2008, together with the values for the known parameters:

We will start in 1999, taking at each simulation a random value in the confidence interval of each parameter. For initial point for the variables X we will take the interval [943, 950], obtained from the end points for the period 1986 -1998, and for Y and Z the known values for persons detected and living with HIV (1612) and living with AIDS (326). Also as a control we will consider the number of persons that have died due to AIDS, this in the form of a new equation to add to the system (2), which we will use in the following form:

dX dt = σX -k 3 X 2 -k 2 XY + λ ′ Y, dY dt = (k 1 + k 3 X)X -γ ′ Y + k 2 X Y, ( 9
)
dZ dt = βX + β ′ Y -µ ′ Z, dM dt = µ ′ Z.
As initial value for M we will use 705.

Results

Local stability of

Q 0 and Q * .
The Jacobian matrix at the point Q 0 is given by:

J 1 (Q 0 ) = ( σ λ ′ k 1 -γ ′ ) Q 0 ,
and therefore P 0 , is locally asymptotically stable (LAS for short) if and only if the trace of J 1 is strictly negative and its determinant strictly positive, i.e.

Q 0 LAS ⇐⇒ σ -γ ′ < 0 and -(σγ ′ + k 1 λ ′ ) > 0 ⇐⇒ σ -γ ′ < 0 and A < 0 ⇐⇒ R 0 < 1. When R 0 > 1, Q 0 is unstable as a saddle equilibrium point. Remember that Q * is X * = γ ′ (k 1 + γ) B (R 0 -1), Y * = X * C 1 C 2 = X * γ ′ λ -γ 1 -R 1 = X * γγ ′ R 2 -1 1 -R 1 .
Q * , and therefore P * , is LAS if and only if it exits (i.e. ( 5) and ( 6) are satisfied or equivalently iff ( 7) is satisfied), the trace of J 1 (Q * ) is strictly negative and its determinant strictly positive, i.e. Q * exists and

τ 1 = Tr(J 1 (Q * )) = σ -γ ′ + k 2 (X * -Y * ) -2k 3 X * < 0 and ∆ 1 = det(J 1 (Q * )) > 0. ( 10
)
After some straight forward calculations we obtain that ∆

1 = A and condition ∆ 1 > 0 is equivalent to R 0 > 1. Thus when Q 0 is LAS, if Q * exists, it is unstable.Subsequently, when R 0 > 1, we have also R 2 > 1 and R 1 < 1 to ensure X * > 0, and Y * > 0. Let us suppose R 0 > 1, R 1 < 1 and R 2 > 1.
Then x = X * > 0 and y = Y * > 0 and

xyτ 1 = (σx -k 3 x 2 -k 2 xy)y -γ ′ xy -k 3 x 2 y + k 2 x 2 y = -λ ′ y 2 -k 3 x 2 y + [y(k 2 x -γ ′ )]x, since σx -k 3 x 2 -k 2 xy = -λ ′ y by (8.1) = -λ ′ y 2 -k 3 x 2 y -(k 1 x + k 3 x 2 )x, since y(k 2 x -γ ′ ) = -k 1 x -k 3 x 2 by (8.2) (11)
which is negative and hence τ 1 < 0. Finally, Q * exists and is LAS iff

R 0 > 1 and R 1 < 1 and R 2 > 1.
Remark: Q 0 being LAS implies that either Q * does not exist in our domain, which is the case if either B > 0 (X * does not exists) or

B < 0 and (R 2 -1)(1 -R 1 ) < 0 (Y * does not exists); or Q * exists, where R 2 < 1 and R 1 > 1, but Q * is unstable.

Global stability of

Q 0 and Q * .
Lemma 1 There are no periodic orbits in D ′ .

Proof:

Let g : (X, Y ) -→ 1 XY , then ∂ ∂X {g(X, Y )X ′ } + ∂ ∂Y {g(X, Y )Y ′ } = - ( λ ′ X 2 + k 3 Y + k1+k 3 X Y 2
) keeps the same sign in D ′ , and using the Dulac criteria we conclude that there are no periodic orbits in the set.

∀t ≥ 0 and ∀k ≥ 0, we set P t = (X(t), Y (t)) and

P k = (X(t k ), Y (t k )) with t 0 = 0. Theorem 1 If R 0 < 1 and R 1 < 1 then Q 0 is globally asymptotically stable (GAS for short) in D ′ . Moreover, when R 0 > 1, Q 0 is unstable. Proof: Q 0 is the unique equilibrium point in D ′ . Let V be the function defined on D ′ by ∀(X, Y ) ∈ D ′ , V (X, Y ) = γ ′ X + λ ′ Y. V is a Lyapunov function for the point Q 0 on D ′ and ∀(X, Y ) ∈ D ′ , V ′ (X, Y ) = γ ′ X [(k 1 + γ)(R 0 -1) + (k 3 X + k 2 Y )(R 1 -1)] . (12) R 1 < 1 then V ′ (X, Y ) < 0 for all (X, Y ) ∈ D ′ . Moreover, V ′ (X, Y ) = 0 iff X = 0, or R 0 = 1 and R 1 = 1. Along the axis {X = 0}, X ′ ≥ 0, we conclude that the largest invariant subset in {(X, Y ) ∈ D ′ | V ′ (X, Y ) = 0} is the singleton {Q 0 }, and from LaSalle's invariant principle we conclude that Q 0 is GAS in D ′ . We have shown previously that if R 0 > 1, Q 0 is a saddle equilibrium point. Therefore, Q 0 is unstable when R 0 > 1. Theorem 2 If R 0 < 1, R 1 > 1 and R 2 < 1
, then the basin of attraction of Q 0 is a triangle formed by the axes and a line that goes through the point Q * and has slope

- λ 2 + γ ′ -k 2 X * λ ′ -k 2 X * .
Proof: Q * is a saddle point with a stable and an unstable manifold both of dimension 1. Let λ 1 be the negative eigenvalue, λ 2 be the positive eigenvalue and E s the eigenspace associated to λ 1 and W s the manifold that is tangent to Q * + E s at each point. E s is the straight line with slope

- λ 2 + γ ′ -k 2 X * λ ′ -k 2 X * .
We know that

1. λ ′ -k 2 X * = λ ′ B-k 2 A B = γ ′ (R 1 -1)(k 2 σ-λ ′ k 3 ) B > 0 because R 1 -1 > 0 and both σ and B are negative. 2. γ ′ -k 2 X * = γ ′ B-k 2 A B = γ ′ (1 -R 1 )(k 1 k 2 + γ ′ k 3 > 0 because 1 -R 1 < 0 and B < 0. Then λ 2 + γ ′ -k 2 X * > 0.
Hence the slope is negative. This means that there is a triangular region formed by the axes and E s that forms the basin of Q 0 . A trajectory that starts in this region can not leave the region because the vector field at the axes points inwards and it cannot cross E s . Therefore this region is invariant and all trajectories starting inside the region must have Q 0 as it is a ω-limit set.

Theorem 3

Suppose that R 0 > 1, R 1 > 1 and R 2 > 1. Then Q 0 is unstable, Q * does not exist and trajectories are unbounded.

Proof:

We have seen before that when R 0 > 1, Q 0 is unstable, and [START_REF] De Arazoza | A non linear model for a sexually transmitted disease with contact tracing[END_REF]

implies that Q * is in D ′ . Let N = X + Y . Then N ′ = X ′ + Y ′ = γ(R 2 -1)X + γ ′ (R 1 -1)Y . Denote a = min(γ(R 2 -1), γ ′ (R 1 -1)). Then a > 0 and N ′ ≥ a(X +Y ) = aN so we have N (t) ≥ N (0)e at , therefore N is not bounded. Let X ′ = f (X, Y ), Y ′ = g(X, Y ), φ : x -→ x k 3 x -σ λ ′ -k 2 x if λ ′ -k 2 x ̸ = 0 and ψ : x -→ x k 3 x + k 1 γ ′ -k 2 x if γ ′ -k 2 x ̸ = 0.
Let us note by C φ and C ψ the curves that represent the functions φ and ψ. Then

X ′ = f (X, Y ) = { (λ ′ -k 2 X)(Y -φ(X)) if λ ′ -k 2 X ̸ = 0 X(σ -k 3 X) otherwise, (13)
and

Y ′ = g(X, Y ) = { (γ ′ -k 2 X)(ψ(X) -Y ) if γ ′ -k 2 X ̸ = 0 X(k 1 + k 3 X) otherwise. ( 14
)
Let x such that λ ′ -k2x ̸ = 0 and γ ′ -k 2 x ̸ = 0, then we have

ψ(x) -φ(x) = x -Bx + A (λ ′ -k 2 x)(γ ′ -k 2 x) (15) Theorem 4 Suppose that R 0 > 1, R 1 > 1 and R 2 < 1.
Then Q 0 is unstable, Q * does not exist and trajectories are unbounded.

Proof:

In this case, we have σ < 0,

λ ′ > γ ′ and λ ≤ γ, A = γ ′ (k 1 + γ)(R 0 -1) > 0 and B = γk 2 (R 2 -1) + γ ′ k 3 (1 -R 1 ) < 0.
Then we have -Bx + A > 0, ∀x ≥ 0. Also, we have seen before that when R 0 > 1, Q 0 is unstable, and [START_REF] De Arazoza | A non linear model for a sexually transmitted disease with contact tracing[END_REF]

leads that Q * is not in D ′ .
Lemma 2

• If 0 < x < γ ′ k 2 , then 0 < φ(x) < ψ(x). • If γ ′ k 2 < x < λ ′ k 2 , then ψ(x) < 0 < φ(x). • If x > λ ′ k 2 , then φ(x) < ψ(x) < 0.
The proof follows directly from [START_REF] Granich | Universal voluntary HIV testing with immediate antiretroviral therapy as a strategy for elimination of HIV transmission: a mathematical model[END_REF]. Remark: When φ and ψ are both nonnegative,

C φ ∩ C ψ = {Q 0 }.
For (x, y) ∈ D ′ , we set

D 1 ={(x, y), x > 0, 0 < y < φ(x)}, D 2 ={(x, y), x > 0, φ(x) < y < ψ(x)}, D 3 ={(x, y), 0 ≤ x < γ ′ k 2 , y > ψ(x)}
and subsequently yield the following lemma.

Lemma 3

1. If P t ∈ D 1 then X ′ < 0 and Y ′ > 0. Moreover, by continuity, if P t ∈ C φ , X ′ = 0 and Y ′ > 0.

If

P t ∈ D 2 , then X ′ > 0 and Y ′ > 0. Moreover, by continuity, if P t ∈ C ψ , X ′ > 0 and Y ′ = 0, if P t ∈ C φ , X ′ = 0 and Y ′ > 0. 3. If P t ∈ D 3 , then X ′ > 0 and Y ′ < 0. Moreover, by continuity, if P t ∈ C ψ , X ′ > 0 and Y ′ = 0.
For proof of Lemma 3:, see Lemma 2 and (15).

First, we remark that by Lemma 3, D 2 is positively invariant under the flow induced by ( 8), but the trajectories leave all bounded subregions since we have for Y = M, Y ′ > 0, ∀M > 0. By using Lemma 3, wa have the following cases:

1. If P 0 ∈ D 1 , then X decreases and Y increases and the trajectory must cross C φ at a time t 1 > 0 and enter region D 2 . But in this region, X also increases and the trajectory follows the curve C φ . Subsequently,

X(t) -→ t→+∞ λ ′ k 2 .
2. If P 0 ∈ D 2 , we have the same conclusion as in the previous case.

3. If P 0 ∈ D 3 , then X increases and Y decreases, hence the trajectory must cross C ψ at a time t 1 > 0 and enter region D 2 . Thus we have the same conclusion as previous cases.

For global stability of Q * , we suppose that R 0 > 1, R 1 < 1, and R 2 > 1. Then Q * exists and is LAS and Q 0 is unstable.

Lemma 4 1. X * ∈]0, γ ′ k 2 [. 2. X * ∈                ]0, λ ′ k 2 [ if σ k 3 < 0 < λ ′ k 2 ] σ k 3 , λ ′ k 2 [ if 0 < σ k 3 < λ ′ k 2 ] λ ′ k 2 , σ k 3 [ if λ ′ k 2 < σ k 3 < γ ′ k 2 ] λ ′ k 2 , γ ′ k 2 [ if γ ′ k 2 < σ k 3 Proof: First observe that λ > γ, λ ′ < γ ′ , and B = k 2 γ(R 2 -1) + k 3 γ ′ (1 -R 1 ) > 0.
Hence we have

X * - γ ′ k 2 = γ ′ (R 1 -1)(k 1 k 2 + γ ′ k 3 ) k 2 B , X * - λ ′ k 2 = γ ′ (1 -R 1 )(k 2 σ -λ ′ k 3 ) k 2 B , X * - σ k 3 = (λ -γ)(λ ′ k 3 -σk 2 ) k 3 B .
Recall that X ′ and Y ′ can be rewritten as in ( 13) and ( 14), respectively. Also ( 15) is equivalent to

ψ(x) -φ(x) = xB X * -x (λ ′ -k 2 x)(γ ′ -k 2 x) . ( 16
)
The sign of this difference is the sign of the product (X * -x)(λ ′ -k 2 x)(γ ′ -k 2 x). When φ and ψ are both nonnegative, C φ ∩ C ψ = {Q 0 , Q * }. We now give the main theorems for the global stability of Q * .

Theorem 5 Suppose that R 0 > 1, R 1 < 1, and R 2 > 1. If σ ≤ 0, then Q * is GAS in D ′ . Proof: For (x, y) ∈ D ′ set D 1,1 = {x > X * , 0 < y < ψ(x)} D 1,2 = {0 < x < X * , 0 < y < φ(x)} D 1 = D 1,1 ∪ D 1,2 D 2 = {0 < x < X * , φ(x) < y < ψ(x)} D 3,1 = {X * < x < λ ′ k 2 , ψ(x) < y < φ(x)} D 3,2 = { λ ′ k 2 < x < γ ′ k 2 , y > ψ(x)} D 3 = D 3,1 ∪ D 2,2 D 4,1 = {0 < x < X * , y > ψ(x)} D 4,2 = {X * < x < λ ′ k 2 , y > φ(x)} D 4 == D 4,1 ∪ D 4,2
We now give a lemma which will be used in the proof of Theorem 5.

Lemma 5 Suppose σ ≤ 0.

(1)

• If 0 < x < λ ′ k 2 , then φ(x) ≥ 0 and ψ(x) ≥ 0. • If λ ′ k 2 < x < γ ′ k 2 , then φ(x) ≤ 0 and ψ(x) ≥ 0. • If x > γ ′
k 2 then φ(x) ≤ 0, and ψ(x) ≤ 0.

(2

) In ]0, γ ′ k 2 [, we have ψ(x) -φ(x)        > 0 if 0 < x < X * or λ ′ k 2 < x < γ ′ k 2 = 0 if x = X * < 0 if X * < x < λ ′ k 2 .
(3)

• If P t ∈ D 1 , then X ′ < 0 and Y ′ > 0. Moreover, by continuity, if P t ∈ C ψ (D 1 ), X ′ < 0 and Y ′ = 0, if P t ∈ C φ (D 2 ), X ′ = 0 and Y ′ > 0. • If P t ∈ D 2 , then X ′ > 0 and Y ′ > 0. Moreover, by continuity, if P t ∈ C ψ , X ′ > 0 and Y ′ = 0, if P t ∈ C φ , X ′ = 0 and Y ′ > 0.
• If P t ∈ D 3 , then X ′ < 0 and Y ′ < 0. Moreover, by continuity, if P t ∈ C ψ , X ′ < 0 and Y ′ = 0 and if P t ∈ C φ , X ′ = 0 and Y ′ < 0.

• If P t ∈ D 4 , then X ′ > 0 and Y ′ < 0. Moreover, by continuity, if P t ∈ C ψ (D 4,1 ), X ′ > 0 and Y ′ = 0, if P t ∈ C φ (D 4,2 ), X ′ = 0 and Y ′ < 0.
The proof of Lemma 5 follows from Lemma 4 and ( 16). We first remark that by Lemma 5, D 2 is a bounded region, positively invariant under the flow induced by [START_REF] Granich | Cuba's national AIDS program. The first decade[END_REF]. D 3 is also positively invariant under the same flow and trajectory entering D 3 delimits a bounded subregion of D 3 , which becomes positively invariant. Indeed, on the line

Y = M, Y ′ < 0, ∀M > Y * . Moreover D 2 ∩ D 3 = {Q * }.
By using Lemma 5, wa have the following situations: Next, for σ > 0, we consider three separate cases. The proof of the following theorem for the first case where σ k 3 < λ ′ k 2 is similar to the proof of Theorem 5.

Theorem 6 Suppose that R 0 > 1, R 1 < 1, and R 2 > 1. If σ > 0 and σ k 3 < λ ′ k 2 , then Q * is GAS in D ′ .
For the second case of γ ′ k 2 < σ k 3 , we have the following theorem:

Theorem 7 Suppose that R 0 > 1, R 1 < 1 and R 2 > 1. If σ > 0 and γ ′ k 2 < σ k 3 , then Q * is GAS in D ′ . Proof: For (x, y) ∈ D ′ set D 1,1 = {X * < x < γ ′ k 2 , φ(x) < y < ψ(x)} D 1,2 = { γ ′ k 2 < x < σ k 3 , y > φ(x)} D 1,3 = {x ≥ σ k 3 , y > 0} D 1 = D 1,1 ∪ D 1,2 ∪ D 1,3 D 2,1 = {0 < x < X * , y < ψ(x)} D 2,2 = {X * < x < γ ′ k 2 , y < φ(x)} D 2 = D 2,1 ∪ D 2,2 D 3,1 = { λ ′ k 2 < x < X * , y > φ(x)} D 3,2 = {X * < x < γ ′ k 2 , y > ψ(x)} D 3 = D 3,1 ∪ D 3,2 D 4,1 = {x < λ ′ k 2 , y > ψ(x)} D 4,2 = { λ ′ k 2 < x < X * , ψ(x) < y < φ(x)} D 4 = D 4,1 ∪ D 4,2 Lemma 6 Suppose σ > 0 and γ ′ k 3 < σ k 3 . (i) • If 0 < x < λ ′ k 2 , then φ(x) ≤ 0 and ψ(x) ≥ 0. • If λ ′ k 2 < x < γ ′ k 2 , then φ(x) ≥ 0 and ψ(x) ≥ 0. • If γ ′ k 2 < x < σ k 3 , then φ(x) ≥ 0 and ψ(x) ≤ 0. • If x > σ k 3 , then φ(x) ≤ 0 and ψ(x) ≤ 0. (ii) For x > 0, we have ψ(x) -φ(x)        > 0 if 0 < x < λ ′ k 2 or X * < x < γ ′ k 2 = 0 if x = X * < 0 if λ ′ k 2 < x < X * or x > γ ′ k 2 .
(iii)

• If P t ∈ D 1 , then X ′ < 0 and Y ′ > 0. Moreover, by continuity, if P t ∈ C ψ , X ′ < 0 and Y ′ = 0, if P t ∈ C φ , X ′ = 0 and Y ′ > 0. • If P t ∈ D 2 , then X ′ > 0 and Y ′ > 0. Moreover, by continuity, if P t ∈ C ψ , X ′ > 0 and Y ′ = 0, if P t ∈ C φ , X ′ = 0 and Y ′ > 0. • If P t ∈ D 3 , then X ′ < 0 and Y ′ < 0. Moreover, by continuity, if P t ∈ C ψ , X ′ < 0 and Y ′ = 0, if P t ∈ C φ , X ′ = 0 and Y ′ < 0. • If P t ∈ D 4 , then X ′ > 0 and Y ′ < 0. Moreover, by continuity, if P t ∈ C ψ , X ′ > 0 and Y ′ = 0 and if P t ∈ C φ , X ′ = 0 and Y ′ < 0.
Proof: See Lemma 4 and (16). By using Lemma 6, wa have the following situations:

• If P 0 ∈ D 4 then X increases and Y decreases and at a time t 1 > 0, the trajectory must cross C ψ and enter D 2 .

• If P 0 ∈ D 3 then X decreases and Y decreases and at a time t 1 > 0, the trajectory must cross C φ and enter D 4 .

• If P 0 ∈ D 2 then X increases and Y increases and at a time t 1 > 0, the trajectory must cross C φ and enter D 1 . The trajectory stays inside a rectangle of vertices {(X 1 , Y 2 ), (X 1 , Y 3 ), (X 2 , Y 2 ), (X 2 , Y 3 )} and as there are no closed orbits the ω-limit set is the point Q * .

• If P 0 ∈ D 1 ,
To complete the the global stability analysis of Q * , we only need the following theorem for the remaining case where λ ′ k 2 < σ k 3 ≤ γ ′ k 2 , the proof of which is similar to that of Theorem 7:

Theorem 8 Suppose that R 0 > 1, R 1 < 1 and R 2 > 1. If σ > 0 and λ ′ k 2 < σ k 3 ≤ γ ′ k 2 , then Q * is GAS in D ′ .
Using Theorems 5-8 together, we conclude that as long as

R 0 > 1, R 1 < 1 and R 2 > 1, Q * is always GAS in D ′ .
The results in this section can be summarized in the following table, that was presented without proof in [START_REF] De Arazoza | Modeling Detection of HIV in Cuba[END_REF].

Table 2: Asymptotic states for the model. "GAS" denotes equilibrium is globally asymptotically stable, "LAS" denotes it is locally asymptotically stable and "N E" denotes equilibrium does not exist in the domain .

R 0 R 1 R 2 Q 0 Q * (X, Y ) → < 1 < 1 < 1 GAS NE Q 0 < 1 < 1 > 1 GAS NE Q 0 < 1 > 1 < 1 LAS unstable Q 0 or unbounded > 1 < 1 > 1 unstable GAS Q * > 1 > 1 < 1 unstable NE unbounded > 1 > 1 > 1 unstable NE unbounded

Application to the Cuban HIV/AIDS data.

We make use of model ( 2) to fit the data of reported HIV and AIDS cases in Cuba. As mentioned earlier in the Introduction, we divide the period 1986-2008 into two separate time periods, namely 1986-1999 and 1999-2008, to take into account of the introduction of the family doctors into the detection system in the 1990's. The initial known HIV-positive and AIDS population sizes in 1986, Y (0) and Z(0), for fitting the first time period are known from the dataset. Subsequently, the corresponding initial values in 1999 for fitting the second time period are taken from the model fitting results of the first time period. Moreover, estimate for λ, β, β ′ , µ, µ ′ and, X(0) has been obtained in previous work ([7, 10, 11, 12]) and are given in Tables 3 and4 (µ and β are omitted from Table 4 because they are the same as in Table 3). For the period 1999-2008 we have The family doctors program was started in the 1990's, first as a pilot project where the family doctors typically did not prescribe HIV testing. It is only after 1999 when detection through family doctors started to take on an important role in the number of yearly detections, at more than 30% of the new detections in a year. Another significant difference between the two periods chosen is λ ′ . In the first period (1986-1999), the sanatorial system was believed to have played an important role in preventing HIV transmission from persons that had already been detected to be living with HIV. Hence we suppose that there is practically no transmission from the persons living with HIV (Y ) and assume λ ′ = 0 for the first period but not for the second period. We fit the model to the data to obtain values for k 1 , k 2 , and k 3 and λ ′ during the second time period only, by minimizing an error function, As traditional optimization methods failed to work properly, we utilize a genetic algorithm approach to find an initial point for starting the optimization method by using a gradient method. In order to compute the standard errors for the parameters, 200 fitting runs were made using different values of the known parameters taken randomly from their confidence interval. Using PET, a software written on MATLAB [START_REF] Barrios | Estimación de Parámetros en Modelos Epidemiológicos de VIH-SIDA en Cuba[END_REF], we obtain the least-square estimates for k 1 and k 2 for each of the two periods of the Cuban HIV epidemic from 1986 to 2008, as well as that of k 3 and λ ′ for the second period, by fitting the Cuban HIV data of the persons known to live with HIV, Y (t), to the model as described previously. The model fit for Y (t) is given in Figure 2. Here the resulting numbers for 1999 obtained from the first stage of estimation using data from 1986-1999 was used as the initial values for the second stage of estimation using the 1999-2008 data. The estimated mean values of the parameters for each of the two periods with 95% confidence intervals, obtained from the 100 best fits, are also given in Tables 3 and4, respectively, together with the intervals used for the parameters estimated from the data base (λ, µ, µ ′ , β, β ′ ).

We further compute the theoretical values of the number of the unknown persons living with HIV, X(t), from the model fit. The theoretical values of X(t) from 1999 to 2008 are given in Figure 3. By comparing the estimation results for the two periods (see Tables 3 and4), we conclude that detection by random screening (k 1 ) improved significantly after 1999, perhaps reflecting the steeper increase in reported cases after 2000 (see Figure 2), while detection via contact tracing (k 2 ) was at a similar level throughout the entire course of the epidemic. Detection by family doctors (k 3 ) was slightly higher than that of contact tracing after 1999 but of similar magnitude. Both the analytical result (of the dynamics) and the data-fitting parameter estimates indicate that random screening was the most effective route of detection, while contact tracing and family played mainly secondary and complementary roles, as had been previously proposed in [START_REF] Hsieh | Correspondence to Universal voluntary HIV testing and immediate antiretroviral therapy[END_REF]. The estimates of the model parameters also allow us to calculate the three reproduction numbers, R 0 , R 1 , and R 2 , for each of the two time periods along with the 95% confidence intervals (see Table 5).

During both time periods, the Cuba HIV epidemic is in the case of R 0 > 1, R 1 < 1 and R 2 > 1. Hence Q 0 is unstable and Q * is globally asymptotically stable, with all trajecto- ries approaching Q * asymptotically. Hence we can conclude that, given no drastic changes in the future, the HIV epidemic in Cuba will tend (in the long term) towards an endemic steady state which we can estimate from our parameter estimates using the expression for P * = (X * , Y * , Z * ) given in Equation (2). That is, assuming no drastic changes in the prevention, transmission, detection, or treatment of HIV in Cuba in the future, the number of persons in Cuba living with HIV that do not know they are infected would approach and linger around 2700. From Figure 3, the theoretical number of unknown persons in Cuba currently living with HIV is estimated to be around 2400 in 2008 through our model fit. Therefore, we speculate with optimism that, at the endemic steady state, the number of persons living with HIV that represent the main core for the spread of the disease in Cuba will not increase drastically in the future.

Simulations

We considered several possible scenarios where the detection system is modified. These were: We compare the values for the different variables of system (9) at the end of the year 2010, taking as base value those given by the system using the values for the parameters fitted to the real values given in Table [START_REF] Minsap | CUBA: Plan Estratégico Nacional ITS/VIH/SIDA 2007-2011[END_REF]. Then we repeat the process changing the values of the parameters according to the scenarios we propose and compare (in percentage) the changes given by the system. We have obtained the following results. The percentage is an increase or decrease according to the sign This case seems to be a good one, the epidemic does not increase very much, the increase of the undetected HIV+ (X) is around 25%, and we can expect that the epidemic will increase, but will not explode in the near future.

As we can see the the size of the epidemic (X + Y + Z) increases by almost 34% in 6 years. It is important to note that even if the detected HIV+ increases (Y ) the unknown HIV+ (X) increases by almost 68%, these persons that do not know that they carry the HIV are the ones that infect the most, making it a very dangerous situation for the future. Also it is important to note that the number of persons that are detected and living with HIV (X + Y + Z) increases more than 25%, this will make a significant impact on the resources that have to be allocated for treatment.The number of deaths due to AIDS also increase, but not dramatically. In this scenario (we assume that only the random detection is carried on), the epidemic is really growing with an increase of almost 50% (X + Y + Z) and an explosion of the unknown infections of 145%. The future in this case is rather bleak.

Next we simulate an overall deterioration of the detection system. this we will do by taking the parameter values n(k 1 , k 2 and k 3 ) at a fraction of their value, we will use a constant p that will be a multiplier to each k i . We can see that any general reduction in the detection effort will lead to an explosion of the epidemic (X + Y + Z), an increase of the undetected cases (X), a decrease in detection, but also an increase in the number of AIDS cases and of the deaths due to AIDS.

Another indicator is what is called coverage and efficiency of the detection system [START_REF] De Arazoza | The HIV/AIDS epidemic in Cuba: Description and tentative explanation of its low HIV prevalence[END_REF]. The efficiency of the system can be measured by the percentage of persons living with HIV that has been detected by the system. From the different scenarios that we have simulated we compute the percentage of the epidemic that is known and we give the results in Table 10.

The efficiency obtained in [START_REF] De Arazoza | The HIV/AIDS epidemic in Cuba: Description and tentative explanation of its low HIV prevalence[END_REF] for the year 2005 was 79.6%, if we use the theoretical values of the model to compute the efficiency of the detection system in 2010 we obtain 76.6%, if we use the number of detected persons living with HIV in 2010 (real data) and the undetected persons living with HIV from the value given by the model we get an efficiency of 85.8%.

Conclusions

For the model proposed we have identified three threshold parameters of epidemiological importance, namely R 0 , R 1 , and R 2 , for the dynamic behavior of our model. A summary of their respective roles in the asymptotic state of the system is given in Table 1. Previous modeling studies with secondary or dual reproduction numbers include, among others, [START_REF] Hsieh | HIV Model with Secondary Contact Tracing: Impact of the Partner Notification Program in Cuba[END_REF][START_REF] Hsieh | Community treatment of HIV-1: Initial and asymptotic dynamics[END_REF][START_REF] Busenberg | A model for HIV in Asia[END_REF][START_REF] Hsieh | Behavior Change and Treatment of Core Group and Bridge Population: Its Effect on the Spread of HIV/AIDS[END_REF][START_REF] Hsieh | The effect of density-dependent treatment/behavior change on the transmission dynamics of HIV[END_REF][START_REF] Hsieh | Basic reproduction number HIV model incorporating commercial sex and behavior change[END_REF].

The biological interpretation of these parameters is given as follows. R 1 = λ ′ /(β ′ + µ) is the number of infections caused by an infective after he/she has been tested positive for HIV and before of AIDS-related illness (ARI), and R 2 = λ γ is the number of infections caused by an infective that is not detected during the asymptomatic period, i.e., if an infective either develops ARI without prior detection or pass away before onset of AIDS. Finally,

R 0 = λ/(k 1 + γ) + λ ′ k 1 /[γ ′ (k 1 + γ)]
gives the number of infections by an infective while living with HIV who had been detected through random screening at his/her asymptomatic stage. R 0 can be considered as the well-known basic reproduction number of the system, provided that we consider the detection of HIV-positive individuals in Cuba as part of routine disease surveillance instead of an intervention measure, especially in light of the more realistic model assumption that those detected to be HIV-positive can still infected others, albeit at a lower level, before the onset of ARI [START_REF] De Arazoza | A non linear model for a sexually transmitted disease with contact tracing[END_REF]. Mathematically, R 0 determines whether the DFE (Q 0 ) is locally asymptotically stable, which is typical of its theoretical significance in asymptotic analysis. R 1 and R 2 are secondary reproduction numbers, which help us to determine the asymptotic behavior of the system pertaining to the endemic equilibrium. In the present model, whether R 1 is larger than unity determines the existence of an unstable endemic equilibrium Q * , provided R 0 < 1. Moreover, when R 0 > 1, R 1 < 1, and R 2 > 1, we have the existence of a globally asymptotically stable equilibrium Q * .

As the threshold parameters of the proposed model, these reproduction numbers and the conditions of asymptotic behaviors we have obtained are epidemiologically meaningful. Our results indicate that random screening is most important as a surveillance tool for HIV, since the number of infections due to an infective detected through random screening, as quantified by R 0 , determines whether the DFE is (locally) asymptotically stable. In other words, if the averaged total number of infections by an infective detected through random screening exceeds one, then there will always be an epidemic. On the other hand, as intervention measures, local asymptotic stability for DFE can be achieved (i.e., R 0 to be brought down to less than one) if: (i) random screening is sufficiently effective (i.e., k 1 is sufficiently large), and (ii) infections by detected HIV-positive individuals are minimal (i.e., λ ′ is sufficiently small). We note that the assumption that the parameters do not change in the long term is not realistic, when we take into account the increasing knowledge regarding the virus and its treatment, vaccines, as well as other social/politial and public health aspects that could affect the dy-namics of the epidemic. For example, it is commonly accepted that existing therapy reduces the probability of transmission for HIV, in some studies it has been assumed to be as low as 1% of the untreated infectives (see e.g., [START_REF] Granich | Universal voluntary HIV testing with immediate antiretroviral therapy as a strategy for elimination of HIV transmission: a mathematical model[END_REF]). In terms of our model, this means that coefficient λ ′ is further reduced which will make R 1 smaller and subsequently altering the value of the asymptotic point. However, therapy delays the onset of AIDS, this implies that β ′ will become smaller, and subsequently causing R 1 to become larger. So long-term predictions (as in the case of asymptotic behavior) make sense only as an indicator of how the epidemic is going to behave with no significant changes in the relevant factors. In this sense, from Table 2 we know that if, for some reasons in the future we have R 1 > 1, trajectories could become unbounded and this would mean that the epidemic is out of control, if no drastic measures are taken. Therefore, from the public health policy perspective it is very important to monitor any changes in the value of R 1 .

Our results further highlight the importance of education of the known infectives for the purpose of achieving behavior change, in light of the second threshold parameter R 1 = λ ′ /γ ′ . If there is no significant behavior change by the known infectious to practice less risky sexual habits and their infection rate is still high (i.e., λ ′ is sufficiently large) resulting in the average number of infections by a known infective (R 1 ) exceeds unity, then the endemic equilibrium is always unstable and there is always a possibility for the total number of infectious (i.e., X +Y ) to increase without bound. This is true even if R 0 < 1, provided that the initial population sizes are outside the domain of attraction of the DFE. This scenario of adverse impact of public health measures, which had been shown previously to be theoretically possible in [START_REF] Hsieh | Community treatment of HIV-1: Initial and asymptotic dynamics[END_REF][START_REF] Anderson | Potential of community-wide chemotherapy or immunotherapy to control the spread of HIV-1[END_REF][START_REF] Hsu | Modeling intervention measures and public response during SARS outbreak[END_REF], is only possible provided k 1 is sufficiently small compared to γ, given that R 0 is a convex combination of R 1 and R 2 (see Equation ( 3)). In other words, an ill-designed detection system without the necessary behavior changes by the detected infectives might adversely lead to the epidemic increasing without bounded if (i) random screening if not sufficiently comprehensive (i.e., k 1 is too small); (ii) lack of education program to change behavior (i.e., λ ′ is too high); and (iii) HIV prevalence was too high when the system was first implemented (initial infective populations outside of the domain of attraction of the DFE). This result further highlights the importance of universal testing in high-prevalence regions, where it has been shown that low infectiousness that accompanied treatment of detected HIV-positives must be drastically reduced to have a lasting impact on HIV prevalence [START_REF] Granich | Universal voluntary HIV testing with immediate antiretroviral therapy as a strategy for elimination of HIV transmission: a mathematical model[END_REF][START_REF] Hsieh | Correspondence to Universal voluntary HIV testing and immediate antiretroviral therapy[END_REF]. In contrast, if through adequate education to achieve sufficiently impacting change behavior of the known infectious (so that λ ′ is sufficiently low) and subsequently the average number of infections by a known infective is lowered down to less than one, then either the DFE Q 0 or the endemic equilibrium Q * will globally asymptotically stable, leading to a more manageable epidemic for the public health authority, even if it is not possible to fully eradicate the disease.

Another aspect is the evaluation ,through the model, of drastic changes in the detection system as shown in tables (6 -9). Any reduction of the detection system (9) produces an increase in the number of undetected persons living with HIV that will change the future of the epidemic as the persons living with HIV that do not know they are infected are more likely to produce new infections, and this will make the total size of the epidemic to increase without control. We can also see that the efficiency of the detection system decreases [START_REF] Lounes | A two-type model for the Cuban national programme on HIV/AIDS[END_REF]. If the detection system is modified, and some of the detection sources are no longer considered, for example eliminate Contact Tracing (k 2 = 0) or the intensive detection through the Family Doctors (k 3 = 0) the efficiency decreases but it is still of the order of 70%. If only random detection is kept (k 2 = k 3 = 0), efficiency could drop to a level of 60%.

Our model shows that the detection system is an intrinsic part of the low prevalence of the Cuban HIV/AIDS epidemic, any change in the detection system could create chaos in the management of the epidemic due to the short term increase of the number of undetected cases and the reduction of the efficiency of the detection system.

Rapatski et al. [START_REF] Rapatski | Mathematical epidemiology of HIV/AIDS in Cuba during the period 1986-2000[END_REF] had suggested that for models of Cuban HIV-AIDS epidemic, it is more realistic for the model parameters to be varying with time. However there are several obstacles with the approach of nonconstant coefficients, not the least of which is due to the fact that there is no reliable estimation for the number of undetected persons living with HIV. That is why in this paper we chose to use parameters that are constant in a certain time interval, or alternatively as step functions (with 2 steps for the present work) that also produce a change in the model itself, by changing the recruitment function from X to Y , F (X, Y ), from a linear polynomial in Y for the period 86-99, to a linear polynomial in X and Y for the period 99-08. The results show that we can obtain better knowledge of the HIV/AIDS epidemic in Cuba, both in terms of the comparative roles played by different methods of detection at different time periods, and in terms of the current status and possible future outlook of the epidemic.

Model Parameters: 1 .

 1 λ: rate of recruitment of new HIV infected persons infected by X, 2. λ ′ : rate of recruitment of new HIV infected persons infected by Y , 3. k 1 : rate at which the unknown HIV infected persons are detected by the system, independently of other HIV-positive persons (through "random" screening), 4. k 2 : rate at which unknown HIV-infected persons are detected by the system through contact tracing, 5. k 3 : rate at which unknown HIV infected persons are detected by the system, through the doctor, 6. β: rate at which the undetected HIV-positive persons develop AIDS, reciprocal of the mean incubation 7. β ′ : rate at which the detected HIV-positive persons develop AIDS, the reciprocal of the mean time it takes to go from Y to Z, 8. µ: mortality rate of the sexually active population, 9. µ ′ : mortality rate of the population with AIDS.

1 .

 1 If P 0 ∈ D 1 , then X decreases and Y increases and at a time t 1 > 0, the trajectory must cross C ψ and enter region D 3 or cross C φ and enter the bounded positively invariant regionD 2 containing only Q * (noting that Q 0 is unstable). If the trajectory enters region D 3 , D 3 ∩ {x > 0, Y = Y (t 1)} determines a positively invariant subregion of D 3 and its contains only Q * . By Lemma 1 there are no periodic orbits, then the Poincaré-Bendixson Theorem entails that P t tends towards Q * .2. If P 0 ∈ D 2 or P 0 ∈ D 3 we have the same conclusion as in the previous case. 3. If P 0 ∈ D 4 then X increases and Y decreases and at a time t 1 > 0, the trajectory must cross C φ and enter region D 3 or cross C ψ and enter region D 2 . The first case allows us to conclude that P t tends towards Q * .

Figure 2 :

 2 Figure 2: Phase plane portrait for the years 1986-2008. The dots (X Model (t), Y Data (t)) are the real data, the solid line (X Model (t), Y Model (t)) denotes the model-generated curve.

Figure 3 :

 3 Figure 3: Model-generated curve for unknown HIV-infected in Cuba during 1986-2008.

1 .

 1 Model flow diagram. 2. Phase plane portrait for the years 1986-2008. The dots (X Model (t), Y Data (t)) are the real data, the solid line (X Model (t), Y Model (t)) denotes the model-generated curve. 3. Model-generated curve for unknown HIV-infected in Cuba during 1986-2008.

Table 1 :

 1 Estimated mean values with the 95% confidence intervals for the known parameters and the fitted parameters in the period 1999-2008. UCI and LCI denote the respective upper and lower bounds for the 95% confidence intervals.

					.
		Known Parameters		Fitted Parameters
		λ	µ ′	β ′	λ ′	k 1	k 2	k 3
	Mean .4137 .116370 .107261 .08274 .2195 2.728×10 -5 6.801×10 -5
	LCI	.4134 .116366 .107258 .08265 .2192 2.725×10 -5 6.783×10 -5
	UCI	.4140 .116373 .107264 .08283 .2198 2.732×10 -5 6.818×10 -5

  then X decreases and Y increases and at a time t 1 > 0, the trajectory must cross C ψ and enter D 3 with X * < X 1 < X 0 and Y 0 < Y 1 . In D 3 , X and Y both decrease and at a certain time t 2 > t 1 , the trajectory must cross C φ and enter D 4 with X 2 < X * and Y 2 < Y 1 . In D 4 , X increases and Y decreases and at a certain time t 3 > t 2 , the trajectory must cross C ψ and enter D 2 with X 2 < X 3 < X * and Y 3 < Y 2 . In D 2 , X increases and Y increases and at a certain time t 4 > t 3 , the trajectory must cross C φ and enter D 1 with X * < X 4 < X 1 and Y 3 < Y 4 < Y 2 . This establishes a circulation flow.

Table 3 :

 3 Estimated mean values with the 95% confidence intervals for the known parameters and the fitted parameters in the period 1986-1999 (k 3 = 0). UCI and LCI denote the respective upper and lower bounds for the 95% confidence intervals.

				.
			Known Parameters		Fitted Parameters
		λ	µ	µ ′	β	β ′	k 1	k 2
	Mean .4139	.001	.3957 .0924 .1168 .1338 2.413×10 -5
	LCI	.4107 .00097 .3947 .0894 .1145 .1335 2.395×10 -5
	UCI	.4172 .00103 .3968 .0953 .1192 .1342 2.431×10 -5

Table 4 :

 4 Estimated mean values with the 95% confidence intervals for the known parameters and the fitted parameters in the period 1999-2008. UCI and LCI denote the respective upper and lower bounds for the 95% confidence intervals.

					.
		Known Parameters		Fitted Parameters
		λ	µ ′	β ′	λ ′	k 1	k 2	k 3
	Mean .4137 .116370 .05682 .08274 .2195 2.728×10 -5 6.801×10 -5
	LCI	.4134 .116366 .05571 .08265 .2192 2.725×10 -5 6.783×10 -5
	UCI	.4140 .116373 .05797 .08283 .2198 2.732×10 -5 6.818×10 -5

Table 5 :

 5 Estimated mean values with the 95% confidence intervals for R 0 , R 1 , and R 2 . UCI and LCI denote the respective upper and lower bounds for the 95% confidence intervals.

			.
		1986-1999	1999-2008
		LCI	UCI	LCI	UCI
	R 0 1.824 1.837 1.858 1.859
	R 1	0	0	0.762 0.765
	R 2 4.446 4.665 4.428 4.434

Table 6 :

 6 Evolution from 2004 to 2010, k 2 = 0.

	X	Y	Z	M	X + Y + Z Y + Z Y + Z + M
	24.49 -10.22 0.24 0.02	1.45	-5.61	-3.93

Table 7 :

 7 Evolution from 2004 to 2010, k 3 = 0.

	X	Y	Z	M	X + Y + Z Y + Z Y + Z + M
	67.97 29.36 16.13 3.50	33.94	23.54	17.57

Table 8 :

 8 Evolution from 2004 to 2010 k 2 = 0 and k 3 = 0.

	X	Y	Z	M	X + Y + Z Y + Z Y + Z + M
	144.95 17.64 20.76 4.22	48.55	19.02	14.61

Table 9 :

 9 Evolution from 2004 to 2010 pk 1 , pk 2 and pk 3 .

	p	X	Y	Z	M	X + Y + Z Y + Z Y + Z + M
	0.75 39.72	-4.64	3.81 0.82	8.61	-0.92	-0.40
	0.50 107.72 -13.75 9.62 2.00	22.61	-3.46	-1.83
	0.25 240.34 -31.23 19.74 3.88	49.64	-8.78	-5.01

Table 10 :

 10 Efficiency in 2010(Model).

	case	Efficiency (%)
	k 3 = 0	70.61
	k 2 = 0	71.23
	k 2 = k 3 = 0	61.35
	p * k i ; p = 0.75	69.84
	p * k i ; p = 0.50	60.28
	p * k i ; p = 0.25	46.68

Acknowledgments.

This work was carried out during visits to the University of Paris Descartes by YHH and H de A. The authors received support from the French "Agence National pour la Recherche" project "Viroscopy". H de A also received support from the Spanish AECID, from their project PCI D/023835/09. YHH is also supported by the National Science Council of Taiwan. For all the support we have received from all the different sources we are grateful.

Conflict of Interest

We declare that we have no conflict of interest.