1. Introduction

The treatment of Human Immunodeficiency Virus-1 (HIV) infection presently faces extraordinary opportunities and challenges in achieving durable efficacy in previously untreated subjects. In fact, even though it does not allow a complete elimination of the virus from several tissue and blood cell reservoirs, antiretroviral treatment has been shown to decrease viral replication as detected by HIV ribonucleic acid (RNA) in plasma [27,35]. The virological control obtained on treatment, documented by undetectable levels of HIV RNA in plasma, has been shown to improve the immunological setting of patients, increasing the number of CD4+ cells, the subset of lymphocytes mainly affected by the virus [23].

Since high rates of HIV RNA control have been shown on treatment both in recent cohort studies and in clinical trials [36], this achievement has been associated with a marked decrease in several opportunistic infections and clinical comorbidities [40,12]. Besides, a successful virological control has also beneficial effects in reducing the spread of HIV infection through sexual contacts, as documented by studies evaluating antiretroviral treatment as a tool for the prevention of HIV transmission to healthy partners [14].

The decrease in CD4 cell counts during HIV infection in the absence of therapy has been used in the past to obtain a marker of immune depletion in order to advise in favor or against the introduction of antiretroviral therapy. Since AIDS-related opportunistic infections generally appear when CD4 are lower than 200 cells$/\mu$L, previous guidelines suggested the initiation of antiretroviral treatment when CD4 decreased below 350 cells$/\mu$L [41]. On the contrary, recent clinical evidences suggest that the early introduction of antiretroviral therapy at diagnosis, regardless of CD4 cell counts, may have benefits on clinical outcomes [37]. Indeed, it is expected that this change in antiretroviral prescription, known as 'test and treat' strategy, will also reduce the spread of HIV infection to the uninfected population [42].

According to estimates by the Italian National Institute of Health, about 120,000 individuals are living with HIV in Italy [7]. A recent investigation reports that 94,146 individuals were in care at public clinical centres at the end of 2012 and 82,501 were receiving antiretroviral treatments [6]. The proportion of patients in therapeutic failure is precisely unknown, even if clinical trial data show that the prevalence of HIV undetectability on treatment can reach 80-90$\%$. The future trends of HIV spread in the Italian population will be determined by the burden of subjects with detectable HIV RNA, thus depending on several factors, such as the appropriate use of treatment options and the effort to anticipate the initiation of antiretroviral therapy. An appropriate prediction of the future trends in HIV epidemics will likely help to enlighten the possible impact of different public health and clinical strategies on the spread of the infection.

Since the seminal paper by May and Anderson [26], a large number of studies have been devoted to mathematical models of HIV epidemic. We mention some recent investigations [17,16,18] that establish global dynamics properties of classes of models of interest for the HIV epidemic. Other mathematical models have been developed either to assess the impact of larger treatment availability on the spread of HIV [21,11], showing a beneficial impact on the reduction of the infection incidence rate, and to assess the effects of virus mutation and of drug resistance onset [2,34,24].

Many papers aim at describing the epidemic evolution in single countries or in single geografic regions (see e.g [33,28,44]). Concerning Italy, Arcà et al. built a detailed multi-stage ODE model for the HIV transmission in Latium (the Italian region around Rome), with eight groups differentiated by sex and routes of infection [1]. With focus on the epidemic among drug users in Latium, Iannelli et al. [19] proposed an age-structured model validated against data on the number of new AIDS cases. Camoni et al. [7], by considering different high-risk subpopulations, and gathering a number of different statistical information, gave an estimate of the history of the epidemic in Italy, taking also into account the practice of therapy. Nevertheless, no predictions are available about the possible impact of an early antiretroviral treatment and 'test and treat' strategies on HIV epidemic in specific European Countries, such as Italy.

In order to forecast the possible impact of new strategies in antiretroviral management, a differential equation model is developed in the present paper to describe the evolution of HIV epidemic in Italy for the years $2003 \div 2025$. Unlike the above mentioned studies, we try to describe the HIV epidemic in Italy by considering a single population of susceptibles, thus making no reference to high-risk groups. The main motivation of our choice is that of keeping the model simple, also in view of future developments about the effects of virus resistance to antiretroviral drugs. The model presented here takes into account the extent of immigration and its effect, on both the general demography and the dynamics of the infected subpopulations.

In Section 2, the ODE model is formulated, and in Section 3 the parameter values are assessed, in part by means of best fitting of available data over the decade $2003 \div 2012$. Predictions until year 2025 are reported in Section 4, showing the benefits of early initiation of the antiretroviral therapy. Section 5 reports some results of the analysis of the model in its time-invariant version, in particular the global stability of the equilibrium points is established either in the absence and in the presence of HIV infection among immigrants. Some concluding remarks (Section 6) complete the work.

---

2. Model formulation

The present model describes the dynamics of HIV infection in Italy by assuming that the heterosexual/homosexual activity and the exchange of needles among drug users are the only significant modalities of HIV transmission. Therefore, the mother-child transmission and the transmission by blood transfusion are disregarded. This choice is supported by recent reports of the Italian National Institute of Health [4,5], indicating that the percentages of new diagnoses in 2012 and 2013 related to transmission routes not involving inter-individual contacts are less than $2\%$.

The HIV epidemic in Italy is still mainly concentrated in high risk subpopulations (injected drug users, male homosexuals, female sex workers, multi-partners heterosexual male). However, there are evidences of an important spreading of the infection outside these groups. Indeed, 32$\%$ of infected people at the end of 2012 was estimated to belong to the remaining general population [7].

On the basis of this evidence and for the sake of simplicity, we assume in our model that the population susceptible of infection, $S(t)$, consists of all the uninfected Italians having age in the range 20 $\div$ 70 years. We choose these age limits for susceptibles since the $97.6\%$ of the new diagnoses recorded in 2012 and 2013 falls in that age range [4,5]. Concerning the Italian demography, we will explicitly take into account the presence of immigration, since immigration in the last years is comparable with the birth rate of the resident population [31]. As a consequence, the role of immigration in the dynamics of HIV infected population cannot be disregarded a priori.

To represent the intra-host disease progression, the untreated infected population is distributed over four compartments in cascade according to the CD4 counts. In particular, according to [1,34], we distinguish the infected population $I(t)$ into the following subpopulations:

• $I_{1}(t)$: number of individuals in the primary infection stage;

• $I_{2}(t)$: infected individuals with CD4 count $> 350$ cells$/\mu$L;

• $I_{3}(t)$: infected individuals with CD4 count in the range 200$\!\div\!$350 cells$/\mu$L;

• $I_{4}(t)$: infected individuals with CD4 count $< 200$ cells$/\mu$L.

Compartment $I_{4}$ represents individuals in the AIDS stage. A further compartment, $I_{5}(t)$, is added to account for infected people under HAART (highly active antiretroviral therapy), introduced in Italy since 1996. Compartment $I_{5}$ includes both patients achieving complete and partial viral suppression [34].

Under the above assumptions, we can describe the dynamics of susceptible and infected subpopulations by the following ODE system (see the block diagram of Figure 1):

$\dot{S}(t) \!=\! \Lambda(t) + (1 - \alpha) \Phi(t)\!-\mu(t)S(t)\!-\!\!\left(\sum\limits_{n=1}^5\!\beta_{n} I_{n}(t)\!\right)\!S(t) \\ \displaystyle \dot{I}_{1}(t) \!=\! a_{1} \alpha \Phi(t)\!+\!\!\left(\sum\limits_{n=1}^5\!\beta_{n} I_{n}(t)\!\right)\!S(t)\!-\large(\theta_{1}\!+\!\mu(t)\large) I_{1}\!(t), \\ \displaystyle \dot{I}_{2}(t) \!=\! a_{2} \alpha \Phi(t)\!+\!\theta_{1} I_{1}\!(t)\!-\large(\theta_{2}\!+\!\mu(t)\!+\!\delta_{2}\large) I_{2}(t), \\ \displaystyle \dot{I}_{3}(t) \!=\! a_{3} \alpha \Phi(t)\!+\!\theta_{2} I_{2}\!(t)\!-\large(\theta_{3}\!+\!\mu(t)\!+\!\delta_{3}\large) I_{3}(t)\!+\!\xi(t) I_{5}(t), \\ \displaystyle \dot{I}_{4}(t) \!=\! \theta_{3} I_{3}\!(t)\!-\large(\theta_{4}\!+\!\mu(t)\!+\!\delta_{4}\large) I_{4}(t), \\ \displaystyle \dot{I}_{5}(t) \!=\! \delta_{2} I_{2}(t)\!+\!\delta_{3} I_{3}(t)\!+\!\delta_{4} I_{4}(t)\!-\!(\!\xi(t)\!+\!\mu(t)) I_{5}(t),$

(1)

where $\Lambda(t)$ is the net input rate in the susceptible population due to the balance between the inflow and the outflow through the selected age boundaries; $\Phi(t)$ is the rate of adult immigration; $\alpha$ is the HIV prevalence among immigrants; $a_{1}$, $a_{2}$, $a_{3}$ are the fractions subdividing the flow $\alpha \Phi(t)$ into the compartments $I_{1}$, $I_{2}$, $I_{3}$ respectively ($a_{1} + a_{2} + a_{3} = 1$); $\mu(t)$ is the per capita loss rate owing to causes not related to the infection (namely, natural death plus emigration); $\beta_{n}$, $n=1, \ldots, 5$, is the infectivity of individuals in compartment $I_{n}$; $\theta_{n}$, $n = 1, 2, 3$ is the per capita disease progression rate from compartment $I_{n}$ to compartment $I_{n+1}$; $\theta_{4}$ is the per capita death rate from AIDS.

Note that, since immigration of people in AIDS stage appears rather unlikely, no input from immigration is assumed into compartment $I_{4}$. Precise data about the presence of immigrants who are under HAART at their arrival in Italy are lacking and official statistical data sources do not report this information. However, observations deriving from clinical experience at the hospital "Luigi Sacco" in Milan indicate a negligible percentage of treated immigrants at the arrival in Italy, mainly concentrated among transexual males from South America. In addition, a survey on a limited number of patients at a clinical centre in Bologna (Italy) evidences that only 14.6$\%$ of immigrants were aware to be infected at their arrival and reports no individuals under treatment [25]. Based on these indications, we suppose no input from immigration into compartment $I_{5}$. Concerning the fractions $a_1$, $a_2$, and $a_3$, we assume

$\displaystyle a_{n} = \frac{1/\theta_{n}}{1/\theta_{1} + 1/\theta_{2} + 1/\theta_{3}}, n = 1, 2, 3,$

(2)

that is $a_{n}$ is taken proportional to the average transit time of the corresponding stage.

The possibility of starting the treatment from different stages of the disease is taken into account by the per capita treatment rates, $\delta_{n}$, of infected individuals in compartment $I_{n}$, $n = 2, 3, 4$. The possibility that treatment is initiated from the primary infection stage is instead disregarded, since at present it is difficult to have a diagnosis of HIV positivity at this stage. $\xi(t)$ is the per capita exit rate from compartment $I_{5}$, and it represents voluntary interruptions or failures of treatment. Although it may be questioned whether the disease in patients exiting from therapy actually evolves as in naive infected individuals, we adopted this simpler view. Therefore we assumed that people exiting from $I_{5}$ return, for simplicity, into one of the stages of the active disease progression. We identify this stage as the stage $I_{3}$, in view of some clinical observations on the time interval typically elapsed from interruption of successful therapy and the falling of CD4 counts below the level of $200$ cells$/\mu$L (2 years in average, [29]).

We assumed in model 1 a bilinear form of the incidence rate, instead of the more usual standard incidence where the infective population is normalized to the total population size, which is $N(t) = S(t) + \sum_{n=1}^{5} I_{n}(t)$ in our case. We note that $N(t)$ is well approximated by $N_{[20, 70]}(t)$, that is the Italian population in the age range $20 \div 70$ years, since $N$ counts in addition only the infected older than 70 years. Thus, denoting by $\bar{N}_{[20, 70]}$ the average value of $N_{[20, 70]}(t)$ over the years 2003 $\div$ 2012, the coefficients $\beta_{n}$, $n = 1, \ldots, 5$, are of the order of the inverse of $\bar{N}_{[20, 70]}$, i.e. $\beta_{n} = O(1/\bar{N}_{[20, 70]})$. The choice of a bilinear infection rate brings about some advantages for the model analysis, as we will see later. Moreover, in our case, the bilinear infection rate is a very good approximation of the standard incidence, and then the dynamics under the two incidence forms are expected to be very similar. In fact, the Italian population in the age range $20\div 70$ years was almost constant from January 1st 2003 to December 31st 2012, changing of only $2\%$ from 38.378 to 39.139 millions of individuals [31].

---

3. Parameter tuning and data fitting

The behaviour of the model has been tested by comparison with epidemiological data over the years $2003 \div 2012$. Precisely, we considered the following data:

• number of infected individuals over the years $2003 \div 2012$, plus the number of treated individuals at the end of 2012, estimated by the UNAIDS tool EPP/Spectrum [7];

• number of treated individuals at the end of 2012 from data of public clinical centres [6];

• number of new AIDS cases over the years $2003 \div 2012$ [5];

• number of deaths from AIDS over the years $2003 \div 2011$ [5].

Data reported by the National Institute of Statistics [31] were used to estimate the time-course of the demographic quantities $\Lambda$, $\Phi$, and $\mu$. The rate of immigration in Italy (per day) is shown in Figure 2 as a piecewise constant function. Although, the age profile of immigrants as a function of time is actually unknown, a survey conducted in 2001 suggests that a percentage of about $70\%$ of immigrants was in the age range $20 \div 70$ years [3]. Thus, we adopted such a percentage to calculate $\Phi(t)$ from the immigration rate of Figure 2. Since no reliable data about the prevalence of HIV infection among immigrants were available, we kept for $\alpha$ the upper bound of the prevalence values estimated in the Italian population at the end of 2012 [7]. The per capita loss rate $\mu(t)$, was taken as a piecewise constant function where the generic constant $\mu_k$ represents the mean loss rate over the year $k$. The value $\mu_k$ was estimated according to

$\mu_k = \frac{M_k}{\bar{N}_{[20, 70]_{k}} \cdot 365},$

where $M_k$ is the number of deaths (not caused by AIDS) plus the number of emigrants during the year $k$, in the age range of interest, and $\bar{N}_{[20, 70]_{k}}$ is the average number of inhabitants of the selected ages during year $k$. The function $\mu(t)$ is plotted in Figure 3. It can be noted that $\mu(t)$ changes only slightly with time, thus we assumed $\mu(t)$ constant and equal to its average value $\bar{\mu} = 1.1 \cdot 10^{-5}$ day$^{-1}$.

Since $N_{[20, 70]}(t)$ is the number of Italians, both non-infected and infected, in the $20 \div 70$ age range, we can write the following approximate balance equation

$\frac{d N_{[20, 70]}(t)}{dt} = \Lambda(t) + \Phi(t) - \mu(t) N_{[20, 70]}(t). \label{eqN}$

(3)

that does not include the loss due to the ageing of infected people beyond age $70$-which can be neglected with respect to the same loss of non-infected individuals taken into account by $\Lambda$-and the deaths from AIDS, totally negligible compared to the non-AIDS deaths.

Eq. 3 allows the estimation of $\Lambda(t)$. For simplicity, we estimated the average value, $\bar{\Lambda}$, of $\Lambda(t)$ over the decade, by a least square fitting of the solution of Equation 3, with $\mu(t) = \bar{\mu}$ and the calculated $\Phi(t)$, to the data on $N_{[20, 70]}(t)$. The estimated $\bar{\Lambda}$ is equal to $-156.36$ individuals per day and the obtained best fitting of $N_{[20, 70]}(t)$ data is reported in Figure 4. The negative value of $\bar{\Lambda}$ reflects the current scarce natality in Italy and the progressive increase of the mean population age.

Concerning the disease progression, the parameters $\theta_{n}$, $n = 1, \ldots, 4$, have been chosen according to [45,34] and are reported in Table 1. Note that $\sum_{n=1}^{4} \theta_{n}^{-1} \simeq 11$ years, which is the typical duration of the disease for untreated individuals [45]. It is well known that the infectivity depends on the viral load of the infectant [32,22] and then on the stage of the disease, with the largest infectivity occurring at the primary stage. We fixed the ratios $\beta_{n}/\beta_{2}$, $n = 1, 3, 4$ following [34], according to the reported viral loads of the different stages. As for the infectivity of treated subjects, different estimates can be found in the literature (see e.g. [15,34,10]). Following again [34], we assumed that patients with complete viral suppression have negligible infectivity. The percentage of treated patients in Italy with viral suppression has been estimated in [6] using CD4 count $> 350$ cells$/\mu$L as a proxy of therapy success, and it amounted to 62.4$\%$ in 2012. From this percentage and the viral load range in treated patients with incomplete viral suppression reported in [34], we could estimate for $\beta_{5}/\beta_{2}$ a value of about $0.20$. This ratio could obviously be lower provided that individuals under therapy kept a behaviour safer than that of untreated people, but a positive impact of ART on risky behaviours is not demonstrated so far [8,20]. A ratio $\beta_{5}/\beta_{2} = 0.20$ was then assumed (see Table 1). However, it must be stressed that this value could be an overestimated and then a precautionary value, particularly when used for the predictions over the period 2013 $\div$ 2025. Indeed, the percentage of viral suppressed patients reported in [6] is probably a pessimistic estimate, in view of the therapeutic achievements of the current clinical practice.

We remark here that the parameters $\theta_{n}$ reflect the intra-host natural history of the disease and so they can reasonably be thought independent of the infection routes. Moreover, whereas the quantitative relationship between probability of transmission per contact and viral load of the infectant can certainly depend on the transmission modality, it is likely that the ratio between the probabilities of transmission is related to the ratio of the viral load of the corresponding infectants in a way that does not depend on the transmission route. Indeed, this independence has been assumed, at least for the sexual transmission in [43,34].

The exit rate from the treatment $\xi(t)$ has been assumed variable in the interval $2003 \div 2013$, with $1/\xi(t)$ ranging linearly from $15$ to $30$ years. The assumption of a decreasing $\xi(t)$ is in agreement with the observed reduction of resistance in Italy in the reference decade [13]. The parameters characterizing $1/\xi(t)$, that is the mean residence time in the compartment $I_{5}$, were calibrated using clinical information and a preliminary trial and error fitting of the number of new AIDS cases and of deaths from AIDS.

The remaining parameters $\beta_{2}$, $\delta_{2}$, $\delta_{3}$, $\delta_{4}$ were estimated from all the available data by means of a weighted least squares procedure, using weights equal to the inverse of the squared data. In order to reduce the effect of the transient response on the fitting curves, the ODE system was integrated starting from year 1996 with an initial condition for $I_{5}$ equal to zero, according to the introduction of HAART in that year. The number of deaths from AIDS in each year was fitted by the integral of $\theta_{4} I_{4}(t)$ over the year of interest, whereas the number of new AIDS cases by the integral of $\delta_{4} I_{4}(t)$. The quantity $\delta_{4} I_{4}(t)$ can be intended as the notification rate of new AIDS cases since the usual practice consists in treating a patient as soon as AIDS is diagnosed. Figures 5-6 show the fitting curves and Table 1 reports all the values of the parameters.

The best fitting of the available data suggests that, in the decade $2003 \div 2012$, the per capita treatment rate of AIDS patients was the largest, followed by the treatment rate of patients with CD4 count between 350 and 200 cells/$\mu$l. Early treatment of patients with higher CD4 counts resulted negligible ($\delta_{2} \approx 10^{-19}$ day$^{-1}$).

---

4. Predictions for the years 2013-2025

We start by predicting the HIV epidemic evolution up to 2025, assuming that $\Lambda(t)$, $\mu(t)$, $\Phi(t)$ and $\xi(t)$ are constant and equal to the corresponding values taken at the end of $2012$. All the remaining parameters are set to the values reported in Table 1. We denote as the reference prediction the system evolution obtained under the mentioned parameter setting.

Next, to investigate the effect of changing the treatment eligibility criterion, we hypothesize two scenarios: a) $\delta_{3}$ is increased to the estimated value of $\delta_{4}$ reported in Table 1; b) both $\delta_{2}$ and $\delta_{3}$ are increased to the value of $\delta_{4}$. For both scenarios, we can observe that anticipating the initiation of treatment results in a lower number either of new AIDS cases and AIDS deaths (see Figure 7). Concerning the evolution of infected people, only scenario b) shows a significant effect (see Figure 8). In particular, it leads to a rapid increase in the number of treated individuals, whereas the number of total infected grows with reduced slope. In scenario b), the model predicts for the year 2025 a 46$\%$ reduction either in the new cases of AIDS and in the number of AIDS deaths, compared to the values of the reference prediction, and a 5$\%$ reduction in the total number of infected individuals. Figure 9 reports the time course of the incidence rate from which we notice in case of scenario b) a 33$\%$ decrease of newly infected subjects at year 2025. We can note that the decreasing trend observed in the period 2003 $\div$ 2012 for either the number of new AIDS cases and the number of AIDS deaths, as well as for the incidence rate, is not maintained over the whole prediction period. This behaviour is likely due to the fact that $\xi(t)$ is assumed no longer decreasing after 2012.

Finally, we investigated the effect of changing the treatment exit rate and the infectivity of treated subjects. In particular, we first assumed $1/\xi(t)$ linearly increasing in 2013 $\div$ 2025 with the same slope of the previous decade. This case, in which $\xi(t)$ is decreasing, represents a continuative improvement of the treatment efficiency. Next, we assumed $1/\xi(t)$ decreasing in 2013 $\div$ 2025, still linearly but with opposite slope of the previous decade. This simulation, with increasing $\xi(t)$, is intended to mimic a relapse of resistance in the prediction period. The results are reported in Figures 10, 11. No apparent effect can be noticed on the evolution of the infected individuals, whereas, as expected, the number of treated individuals decreases when $\xi(t)$ is supposed to increase. By contrast, either for decreasing and for increasing $\xi(t)$, the number of new cases of AIDS and the number of AIDS deaths are remarkably affected (see Figure 11). In particular, for both these quantities, the decreasing trend established over $2003 \div 2013$ is maintained if $\xi(t)$ is decreasing whereas it is dramatically reverted if $\xi(t)$ is supposed to increase.

The effect of changing the infectivity of treated subjects is shown in Table 2. The values of some significant quantities at January 1st 2025 are reported for $\pm 50\%$ variations of the baseline value of the ratio $\beta_{5}/\beta_{2}$, and keeping all the other parameters fixed to the values of Table 1. Note that with $\beta_{5}/\beta_{2} = 0.1$, the HIV incidence rate shows a 34$\%$ reduction compared to its reference value. The reduction of the infectivity of treated subjects appears to work in synergy with the early treatment: when $\beta_{5}/\beta_{2} = 0.1$ and $\delta_{2} = \delta_{3} = \delta_{4}$, the incidence rate at 2025 is reduced indeed by 63$\%$ (see also Figure 9, dashed line).

---

5. Analysis of asymptotic properties

In this section, we establish the stability properties of the time-invariant version of the model described by Eqs. 1, where the functions $\Lambda(t)$, $\mu(t)$, $\Phi(t)$, $\xi(t)$ are assumed constant. Namely, we will study the ODE system:

$\dot{S}(t) \!=\! \Lambda + (1 - \alpha) \Phi\!-\mu S(t)\!-\!\!\left(\sum\limits_{n=1}^5\!\beta_{n} I_{n}(t)\!\right)S(t) \\ \displaystyle \dot{I}_{1}(t) \!=\! a_{1} \alpha \Phi\!+\!\!\left(\sum\limits_{n=1}^5\!\beta_{n} I_{n}(t)\!\right)S(t)\!-\large(\theta_{1}\!+\!\mu\large) I_{1}\!(t), \\ \displaystyle \dot{I}_{2}(t) \!=\! a_{2} \alpha \Phi\!+\!\theta_{1} I_{1}\!(t)\!-\large(\theta_{2}\!+\!\mu\!+\!\delta_{2}\large) I_{2}(t), \\ \displaystyle \dot{I}_{3}(t) \!=\! a_{3} \alpha \Phi\!+\!\theta_{2} I_{2}\!(t)\!-\large(\theta_{3}\!+\!\mu\!+\!\delta_{3}\large) I_{3}(t)\!+\!\xi I_{5}(t), \\ \displaystyle \dot{I}_{4}(t) \!=\! \theta_{3} I_{3}\!(t)\!-\large(\theta_{4}\!+\!\mu\!+\!\delta_{4}\large) I_{4}(t), \\ \displaystyle \dot{I}_{5}(t) \!=\! \delta_{2} I_{2}(t)\!+\!\delta_{3} I_{3}(t)\!+\!\delta_{4} I_{4}(t)\!-\!(\!\xi\!+\!\mu) I_{5}(t),$

(4)

with $0 \leq \alpha \leq 1$ and with the assumption

$\Lambda + (1 - \alpha) \Phi > 0. \label{condition}$

(5)

The above assumption prevents the susceptible population from the extinction. Indeed, if $S=0$, 5 implies $\dot S >0$. Note that 5 is satisfied with $\Lambda = \bar{\Lambda}$, $\Phi = \Phi(t=2013)$ (see Table 1).

We observe that the set

$\Omega = \left\{(S, I_{1}, \ldots, I_{5}) \in I\!\!R^{6}_{+} \, | \, S + \sum\limits_{n=1}^{5} I_{n} \leq \displaystyle \frac{\Lambda + \Phi}{\mu} \right\}, \label{Omega}$

(6)

where $I\!\!R^{6}_{+}$ denotes the non-negative orthant of $I\!\!R^{6}$, is invariant for the system 4. This property can be easily established by examining the directions of the vector field defined by 4 on the boundary of $I\!\!R^{6}_{+}$ and taking into account that

$\dot{N}(t) =\dot S(t) + \sum\limits_{n=1}^5 \dot I_n(t) = \Lambda + \Phi - \mu N(t) - \theta_{4} I_{4}(t) \nonumber \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \leq \Lambda + \Phi - \mu N(t), \label{InequalityN}$

(7)

for non-negative $I_{4}(t)$. The above inequality, in fact, implies

$N(t) \leq N(t_0) {\rm e}^{- \mu (t-t_0)} + \displaystyle \frac{\Lambda + \Phi}{\mu} \left(1 - {\rm e}^{- \mu (t-t_0)}\right), \label{UpperBoundN}$

(8)

which in turn implies $N(t)$ bounded by $(\Lambda + \Phi)/{\mu}$ for any initial state in $\Omega$.

Also, we remark that if an equilibrium point of system 4 exists in $I\!\!R_+^6$, it necessarily belongs to $\Omega$. Let $E$ be such an equilibrium, and $N^{*}$, $I_{4}^{*}$ be the corresponding values of the total population and of the number of individuals in the AIDS stage respectively. Since $E$ is an equilibrium, we have

$\Lambda + \Phi - \mu N^* - \theta_{4} I_{4}^{*}=0,$

(9)

whereas if $E$ were such that $N^{*} > (\Lambda + \Phi)/\mu$ we would get the contradiction:

$\Lambda + \Phi - \mu N^* - \theta_{4} I_{4}^{*} < 0.$

---

5.1. Existence and stability of equilibria for $\mathbf{\alpha = 0}$

Denoting by $E = (S^{*}, I_{1}^{*},$ $\ldots, I_{5}^{*})$ a generic equilibrium point of system 4 with $\alpha = 0$, i.e. without infected among immigrants, we have that $E$ must satisfy the following algebraic system

$\begin{array}{l} \displaystyle \Lambda + \Phi = \mu S^{*} + \left(\sum\limits_{n = 1}^{5}\beta_{n} I_{n}^{*}\right)S^{*}, \\ \displaystyle Q_{1}I_{1}^{*} = \left(\sum\limits_{n = 1}^{5}\beta_{n} I_{n}^{*}\right)S^{*}, \\ Q_{2} I_{2}^{*} = \theta_{1} I_{1}^{*}, \\ Q_{3} I_{3}^{*} = \theta_{2} I_{2}^{*} + \xi I_{5}^{*}, \\ Q_{4} I_{4}^{*} = \theta_{3} I_{3}^{*}, \\ Q_{5} I_{5}^{*} = \delta_{2} I_{2}^{*} + \delta_{3} I_{3}^{*} + \delta_{4} I_{4}^{*}, \end{array} \label{EqulibriumSystem}$

(10)

where

$\begin{array}{l} Q_{1} = \theta_{1} + \mu, \\ Q_{2} = \theta_{2} + \mu + \delta_{2}, \\ Q_{3} = \theta_{3} + \mu + \delta_{3}, \\ Q_{4} = \theta_{4} + \mu + \delta_{4}, \\ Q_{5} = \xi + \mu. \end{array} \label{Q_n_quantities}$

(11)

From the last four equations in 10 we obtain:

$I_{n}^{*} \!=\! m_{n} I_{1}^{*}, n = 2, \ldots, 5, \label{X_Y_Z_T_EquilibriumRealtions}$

(12)

where

$\begin{array}{l} m_{2} \!=\! \displaystyle \frac{\theta_{1}}{Q_{2}}, \\ m_{3} \!=\! \displaystyle \frac{\theta_{2} + \delta_{2} (\xi/Q_{5})}{Q_{3} - \delta_{3} (\xi/Q_{5}) - \theta_{3} (\delta_{4}/Q_{4}) (\xi/Q_{5})} m_{2}, \\ m_{4} \!=\! \displaystyle \frac{\theta_{3}}{Q_{4}} m_{3}, \\ m_{5} \!=\! \displaystyle \frac{\delta_{2}}{Q_{5}} m_{2} + \frac{\delta_{3}}{Q_{5}} m_{3} + \frac{\delta_{4}}{Q_{5}} m_{4}. \end{array} \label{m_n_quantities}$

(13)

Introducing the notation

$\eta = \beta_{1} + \beta_{2} m_{2} + \beta_{3} m_{3} + \beta_{4} m_{4} + \beta_{5} m_{5}, \label{eta}$

and taking into account 12, the first two equations of the system 10 can be rewritten as

$\begin{array}{l} \Lambda + \Phi - \mu S^{*} - \eta I_{1}^{*} S^{*} = 0, \\ I_{1}^{*} (\eta S^{*} - Q_{1}) = 0. \end{array} \label{PS_system_NOphi}$

(14)

System 14 admits two different solutions

$\begin{array}{l} S^{*} \!=\! \displaystyle \frac{\Lambda + \Phi}{\mu}, \, \, I_{1}^{*} \!=\! 0 \end{array} \ \ \ \ {\rm or} \ \ \ \ \ S^{*} \!=\! \displaystyle \frac{Q_{1}}{\eta}, \, \, I_{1}^{*} \!=\! \displaystyle \left(\!\frac{\Lambda + \Phi}{Q_{1}}\!-\!\frac{\mu}{\eta}\!\right), \label{2equilibria}$

(15)

meaning that system 4 has in $\Omega$ the disease-free equilibrium $E_{df}$

$S^{*} \!=\! \displaystyle \frac{\Lambda + \Phi}{\mu}, I_{n}^{*} \!=\! 0, \, \, n = 1, \ldots, 5, \label{E_DiseaseFree}$

(16)

and, if and only if $(\Lambda + \Phi) \eta > \mu Q_{1}$, the positive endemic equilibrium $E_{en}$

$S^{*} \!=\! \displaystyle \frac{Q_{1}}{\eta}, I_{n}^{*} \!=\! \displaystyle m_{n}\!\left(\!\frac{\Lambda + \Phi}{Q_{1}}\!-\!\frac{\mu}{\eta}\!\right), \, \, n = 1, \ldots, 5, \label{E_Endemic}$

(17)

where $m_{1} = 1$.

The stability analysis of the disease-free equilibrium has been performed by means of the reproduction number $R$, computed exploiting the method proposed in [38], after verifying that our dynamical model satisfies all the structural hypothesis required. Then, let us define the vectors

$f = \left[\!\! \begin{array}{c} \displaystyle \!\left(\sum\limits_{n = 1}^{5} \beta_{n} I_{n}(t)\!\right)S(t) \\ 0 \\ 0 \\ 0 \\ 0 \end{array} \!\!\right], \label{f}$

(18)

$v = \left[ \begin{array}{c} Q_{1} I_{1}(t) \\ - \theta_{1} I_{1}(t) + Q_{2} I_{2}(t) \\ - \theta_{2} I_{2}(t) - \xi I_{5}(t) + Q_{3} I_{3}(t) \\ - \theta_{3} I_{3}(t) + Q_{4} I_{4}(t) \\ - \delta_{2} I_{2}(t) - \delta_{3} I_{3}(t) - \delta_{4} I_{4}(t) + Q_{5} I_{5}(t) \end{array} \right], \label{v}$

(19)

where the generic entry $f_{i}$ represents the rate of appearance of new infections in compartment $i$ while the entry $v_{i}$ represents the sum of the transfer rates of individuals exiting from (positive sign) and entering (negative sign) compartment $i$. Let us also compute the Jacobian matrices of vectors $f$, $v$ at the equilibrium $E_{df}$

$F = \nabla{f}|_{E_{df}} = \displaystyle \frac{\Lambda + \Phi}{\mu} \! \left[ \begin{array}{ccccc} \beta_{1}&\beta_{2}&\beta_{3}&\beta_{4}&\beta_{5} \\ 0&0&0&0&0 \\ 0&0&0&0&0 \\ 0&0&0&0&0 \\ 0&0&0&0&0 \\ \end{array} \right], \label{F}$

(20)

$V = \nabla{v}|_{E_{df}} = \left[ \begin{array}{ccccc} Q_{1}&0&0&0&0 \\ - \theta_{1}&Q_{2}&0&0&0 \\ 0&- \theta_{2}&Q_{3}&0&- \xi \\ 0&0&- \theta_{3}&Q_{4}&0 \\ 0&- \delta_{2}&- \delta_{3}&- \delta_{4}&Q_{5} \end{array} \right]. \label{V}$

(21)

According to [9], the reproduction number is defined as the spectral radius of the next generation matrix $F V^{-1}$. So, taking into account the expression of the inverse matrix of $V$,

$V^{-1} = \frac{1}{Q_{1}} \left[ \begin{array}{ccccc} 1&*&*&*&* \\ m_{2}&*&*&*&* \\ m_{3}&*&*&*&* \\ m_{4}&*&*&*&* \\ m_{5}&*&*&*&* \end{array} \right], \label{V-1}$

(22)

where only the first column of $V^{-1}$ is needed for the computation of $R$, and recalling the expression of $F$ given by 20, we obtain

$R = \rho \left(F V^{-1}\right) = \displaystyle \frac{\Lambda + \Phi}{\mu} \frac{\eta}{Q_{1}}, \label{R}$

(23)

where $\rho(\cdot)$ is the spectral radius operator. Note that condition 5 guarantees $R > 0$.

It is useful to rewrite the expression of the endemic equilibrium in terms of $R$, given by 23, as

$S^{*} \!=\! \displaystyle \frac{\Lambda + \Phi}{\mu} \frac{1}{R}, I_{n}^{*} \!=\! \displaystyle m_{n}\!\displaystyle \frac{Q_{1}}{\mu (R - 1)}, \, \, n = 1, \ldots, 5, \label{E_Endemic_2}$

(24)

showing that the endemic equilibrium exists if and only if $R>1$.

Thanks to the bilinear form of the incidence rate, it can be verified that model 4 with $\alpha = 0$ belongs to the class of models studied by Guo et al. in [18]. Then, their stability analysis of the disease-free equilibrium and of the possible endemic equilibrium applies to our model. We have the following results.

Theorem 5.1. When $\alpha = 0$ and $\Lambda + \Phi > 0$, $R$ is positive and the disease-free equilibrium $E_{df}$ exists for any value of $R$. For $R \leq 1$, $E_{df}$ is globally asymptotically stable in $\Omega$, while for $R > 1$ is unstable.

proof. See the proof of Theorem 4.1 in [18]. It can be easily verified that our model satisfies all the basic assumptions ($H_{1}$)-($H_{8}$), as well as the additional hypotheses ($A_{1}$), ($A_{2}$) assumed by the Authors to prove the global stability of the disease-free equilibrium for $R \le 1$ and its instability for $R>1$ (see Appendix).

Theorem 5.2. When $\alpha = 0$ and $\Lambda + \Phi > 0$, the endemic equilibrium $E_{en}$ exists if and only if $R > 1$, and it is globally asymptotically stable in $\Omega \setminus \{E_{df}\}$.

Proof. The proof of the existence comes directly from Eq. 24. Indeed, vector 24 is positive if and only if $R > 1$. For the global stability of the endemic equilibrium, see the proof of Theorem 5.1 in [18]. The additional hypothesis assumed in [18] hold in our case (see Appendix).

In view of the analysis of model 4 with $\alpha > 0$, however, we report here how the proof of the general result in [18] can be specialized in our case. The global stability of $E_{en}$ can be indeed proved considering the following Lyapunov function

$L(S, I_{1}, \dots, I_{5}) = c_{1} \left(\!S - S^{*} - S^{*} \ln\left(\!\frac{S}{S^{*}}\!\right)\!\!\right) \!+\! \sum\limits_{n = 1}^{5} c_{n} \left(\!I_{n} - I_{n}^{*} - I_{n}^{*} \ln\left(\!\frac{I_{n}}{I_{n}^{*}}\!\right)\!\!\right), \label{LyapunovFunction}$

(25)

which, for positive coefficients $c_{n}$, has the following properties:

$\begin{array}{l} L(S, I_{1}, \dots, I_{5}) = 0 \, \Leftrightarrow \, (S, I_{1}, \dots, , I_{5}) = (S^{*}, I_{1}^{*}, \dots, I_{5}^{*}), \\ L(S, I_{1}, \ldots, I_{5}) > 0, \forall (S, I_{1}, \dots, I_{5}) \neq (S^{*}, \dots, I_{5}^{*}). \end{array} \label{LyapunovFunction_Properties}$

(26)

From Eq. 25 it follows

$\dot{L}(S, I_{1}, \ldots, I_{5}) = c_{1} \left(1 - \frac{S^{*}}{S}\right) \dot{S} + \sum\limits_{n = 1}^{5} c_{n} \left(1 - \frac{I_{n}^{*}}{I_{n}}\right) \dot{I}_{n}. \label{DerivLyapunovFunction}$

(27)

By exploiting either the equations of the dynamic system 4 and the equilibrium equations in 10, and introducing also the compact notation $s = S/S^{*}$, $i_{n} = I_{n}/I_{n}^{*}$, $n = 1, \ldots, 5$, we obtain

$\begin{split} \displaystyle &\dot{L}(S, I_{1}, \ldots, I_{5}) \\ \displaystyle =&\, c_{1}\!\left(\!\mu S^{*}\!\left(\!2 - s - \frac{1}{s}\right)\! + \!\sum\limits_{n = 1}^{5} \!\beta_{n} S^{*}I_{n}^{*}\! \left(\!2 - \frac{s \, i_{n}}{i_{1}} - \frac{1}{s} - i_{1} + i_{n}\!\right)\!\!\!\right) \\ \displaystyle&+ c_{2}\!\left(\!\theta_{1} I_{1}^{*}\!\left(\!1 - \frac{i_{1}}{i_{2}} - i_{2} + i_{1}\!\right)\!\!\right) \\ \displaystyle&+ c_{3}\!\left(\!\theta_{2} I_{2}^{*}\!\left(\!1 - \frac{i_{2}}{i_{3}} - i_{3} + i_{2}\right)\! + \xi I_{5}^{*}\!\left(\!1 - \frac{i_{5}}{i_{3}} - i_{3} + i_{5}\right)\!\!\right) \\ \displaystyle&+ c_{4}\!\left(\theta_{3} I_{3}^{*}\!\left(\!1 - \frac{i_{3}}{i_{4}} - i_{4} + i_{3}\right)\!\!\right) \\ \displaystyle&+ c_{5}\!\left(\!\delta_{2} I_{2}^{*}\!\left(\!\!1 - \frac{i_{2}}{i_{5}} - i_{5} + i_{2}\!\right)\!\!+\!\delta_{3} I_{3}^{*}\!\left(\!\!1 - \frac{i_{3}}{i_{5}} - i_{5} + i_{3}\!\right)\!\!+\!\delta_{4} I_{4}^{*}\!\left(\!\!1 - \frac{i_{4}}{i_{5}} - i_{5} + i_{4}\!\right)\!\!\right)\!. \end{split} \label{DerivLyapunovFunction2}$

(28)

Such a function is non-positive in the interior of $\Omega$ by choosing the Lyapunov coefficients $c_1, \ldots, c_5$ according to the method proposed in [18] (see the Appendix for details), namely

$\begin{array}{l} \displaystyle c_{1} = \Theta_{1} \Theta_{3}\!\!\left(\!\Theta_{2}\!\sum\limits_{n = 2}^{4}\!\Delta_{n}\! + \!\Xi_{5} \Delta_{2}\!\right), \\ \displaystyle c_{2} = \!\sum\limits_{n = 2}^{5}\!B_{n} \Theta_{3}\!\!\left(\!\Theta_{2}\!\sum\limits_{n = 2}^{4}\!\Delta_{n}\! + \!\Xi_{5} \Delta_{2}\!\right), \\ \displaystyle c_{3} = \Theta_{1} \Theta_{3}\!\!\left(\!B_{3}\!\sum\limits_{n = 2}^{4}\!\Delta_{n}\! + \!B_{4}\!\sum\limits_{n = 2}^{4}\!\Delta_{n}\! + \!B_{5}\!\sum\limits_{n = 3}^{4}\!\Delta_{n}\right), \\ \displaystyle c_{4} = \Theta_{1} \!\!\left(\!B_{3} \Xi_{5} \Delta_{4} + \!B_{4}\!\!\left(\!\Theta_{2}\!\sum\limits_{n = 2}^{4}\!\Delta_{n}\! + \!\Xi_{5} \!\bigg(\!\Delta_{2} + \Delta_{4}\!\bigg)\!\!\right)\! + \!B_{5} \Delta_{4} (\Theta_{2} + \Xi_{5})\right), \\ \displaystyle c_{5} = \Theta_{1} \Theta_{3}\!\!\left(\!\Xi_{5}\!\sum\limits_{n = 3}^{5}\!B_{n}\! + \!\Theta_{2} B_{5}\!\right), \end{array} \label{LyapunovCoefficients}$

(29)

where

$\begin{array}{l} \Theta_{1} = \theta_{1} I_{1}^{*}, \Theta_{2} = \theta_{2} I_{2}^{*}, \Theta_{3} = \theta_{3} I_{3}^{*}, \\ \Delta_{2} = \delta_{2} I_{2}^{*}, \Delta_{3} = \delta_{3} I_{3}^{*}, \Delta_{4} = \delta_{4} I_{4}^{*}, \\ \, \Xi_{5} = \xi I_{5}^{*}, \, \, B_{n} = \beta_{n} S^{*} I_{n}^{*}, \, n = 2, \ldots, 5. \end{array} \label{LyapunovCoefficients2}$

(30)

Therefore, we obtain $\dot{L}(S, I_{1}, \dots, I_{5}) \leq 0$, with $\dot{L}(S, I_{1}, \ldots, I_{5}) = 0$ only in $\{E_{en}\}$ (see Appendix). This proves the asymptotic stability of $E_{en}$ in the interior of $\Omega$. Since on the boundary of $\Omega$ out of the disease-free equilibrium the vector field of our system points inwards, the basin of attraction of $E_{en}$ can be extended to $\Omega \setminus \{E_{df}\}$.

In order to illustrate the role of antiviral treatment, it may be of interest to compute the reproduction number $R$ for different treatment scenarios. With $\Lambda = \bar{\Lambda}$, $\Phi = \Phi(t=2013)$, $\alpha=0$, and setting the values of the other parameters as in Table 1, $R$ is equal to 2.004. In the absence of treatment ($\delta_n=0$, $n=2, 3, 4$), $R$ decreases to 0.678. So we notice that, in general, treatment is not necessarily accompanied by a lower value of $R$, since treatment implies a longer average life-time of infected individuals and then may imply a greater per capita number of infections. With $\xi=0$ and keeping the other parameter values of Table 1 (i.e. in case of life-long therapy), $R$ increases to 2.796.

---

5.2. Existence and stability of positive equilibria for $\mathbf{\alpha > 0}$

In this section we report the equilibrium analysis of model 4 when $\alpha > 0$, i.e. in the presence of infected among immigrants.

For the generic equilibrium point, the following algebraic system holds

$\begin{array}{l} \displaystyle \Lambda + (1 - \alpha) \Phi = \mu S^{*} + \left(\sum\limits_{n = 1}^{5}\beta_{n} I_{n}^{*}\right)S^{*}, \\ \displaystyle Q_{1}I_{1}^{*} = \left(\sum\limits_{n = 1}^{5}\beta_{n} I_{n}^{*}\right)S^{*} + a_{1} \alpha \Phi, \\ Q_{2} I_{2}^{*} = a_{2} \alpha \Phi + \theta_{1} I_{1}^{*}, \\ Q_{3} I_{3}^{*} = a_{3} \alpha \Phi + \theta_{2} I_{2}^{*} + \xi I_{5}^{*}, \\ Q_{4} I_{4}^{*} = \theta_{3} I_{3}^{*}, \\ Q_{5} I_{5}^{*} = \delta_{2} I_{2}^{*} + \delta_{3} I_{3}^{*} + \delta_{4} I_{4}^{*}, \end{array} \label{EqulibriumSystem_Phi}$

(31)

where the quantities $Q_{n}$, $n = 1, \ldots, 5$, are given by 11. From the last four equations in 31 we obtain:

$I_{n}^{*} \!=\! q_{n} + m_{n} I_{1}^{*}, n = 2, \ldots, 5, \label{X_Y_Z_T_EquilibriumRealtions_Phi}$

(32)

where the coefficients $m_{n}$ are explicitly given in 13, while the terms $q_{n}$ are defined as follows:

$\begin{array}{l} q_{2} \!=\! \displaystyle \frac{a_{2} \alpha \Phi}{Q_{2}}, \\ q_{3} \!=\! \displaystyle \frac{a_{3} \alpha \Phi}{Q_{3}\!-\!\delta_{3} (\xi/Q_{5})\!-\!\theta_{3} (\delta_{4}/Q_{4})\!(\xi/Q_{5})}\!\!+\!\!\displaystyle \frac{\theta_{2} + \delta_{2} (\xi/Q_{5})}{Q_{3}\!-\!\delta_{3} (\xi/Q_{5})\!-\!\theta_{3} (\delta_{4}/Q_{4})\!(\xi/Q_{5})} q_{2}, \\ q_{4} \!=\! \displaystyle \frac{\theta_{3}}{Q_{4}} q_{3}, \\ q_{5} \!=\! \displaystyle \frac{\delta_{2}}{Q_{5}} q_{2} + \frac{\delta_{3}}{Q_{5}} q_{3} + \frac{\delta_{4}}{Q_{5}} q_{4}. \end{array} \label{q_n_quantities}$

(33)

Theorem 5.3. When $\alpha > 0$ and $\Lambda + (1 - \alpha) \Phi > 0$, a unique positive equilibrium exists in $\Omega$ and it is globally asymptotically stable.

Proof. Recalling the definition of $\eta$ given in Section 5.1 and introducing the notation

$\gamma = \beta_{2} q_{2} + \beta_{3} q_{3} + \beta_{4} q_{4} + \beta_{5} q_{5},$

from the first two equations of 31 and from equations 32 we obtain

$S^{*} \!=\! \frac{\Lambda + (1 - \alpha) \Phi}{\gamma + \eta I_{1}^{*} + \mu}, \label{S_EquilibriumRelations}$

(34)

and the quadratic equation in the variable $I_{1}^{*}$

$a {I_{1}^{*}}^{2} + b I_{1}^{*} + c = 0, \label{PE_equation}$

where

$\begin{array}{l} a = - Q_{1} \eta, \\ b = (\Lambda + (1 + a_{1} \alpha - \alpha) \Phi) \eta - Q_{1} (\gamma + \mu), \\ c = a_{1} \alpha \Phi \mu + (\Lambda + (1 + a_{1} \alpha - \alpha) \Phi) \gamma. \end{array}$

From 5 we have $c > 0$, and then, since $a < 0$, there exists a unique positive solution for $I_{1}^{*}$, that is

$\begin{array}{l} I_{1}^{*} = \zeta + \sqrt{\zeta^{2} + \nu}, \\ \displaystyle \zeta = \frac{(\Lambda + (1 + a_{1} \alpha - \alpha) \Phi) \eta - Q_{1} (\gamma + \mu)}{2 Q_{1} \eta}, \\ \displaystyle \nu = \frac{a_{1} \alpha \Phi \mu + (\Lambda + (1 + a_{1} \alpha - \alpha) \Phi) \gamma}{Q_{1} \eta}. \end{array} \label{PE}$

(35)

Note that if a positive equilibrium exists, Eq. 34 implies that 5 must hold. Thus, 5 is a necessary and sufficient condition for the existence of the unique positive equilibrium defined by Eqs. 34, 35, and 32.

The stability of such an endemic equilibrium can be proved considering the same Lyapunov function defined by 25. Indeed, from the derivative of $L$ given by 27, from system 4 with $\alpha >0$, and from the equilibrium equations 31, we obtain

$\dot{L}(S, I_{1}, \ldots, I_{5}) = \dot{L}^{\prime}(S, I_{1}, \ldots, I_{5}) + \dot{L}^{\prime\prime}(S, I_{1}, \ldots, I_{5}), \label{DerivLyapunovFunction_Phi}$

(36)

where

$\dot{L}^{\prime}(S, I_{1}, \ldots, I_{5})\! = \! \sum\limits_{n = 1}^{3} \! c_{n} a_{n} \alpha \Phi \left(2 - i_{n} - \frac{1}{i_{n}} \right), \label{DerivLyapunovFunction_1term}$

(37)

while the term $\dot{L}^{\prime\prime}(S, I_{1}, \ldots, I_{5})$ equals the r.h.s. of 28. Since in the proof of Theorem 5.2 it has been demonstrated the non-positivity of the function 28, we have $\dot{L}^{\prime\prime}(S, I_{1}, \ldots, I_{5}) \leq 0$ for the same choice of the coefficients $c_{n}$ given in 29, 30. Conversely, it is easy to see that for any (positive) choice of $c_{n}$, it results $\dot{L}^{\prime}(S, I_{1}, \dots, I_{5}) \leq 0$ since it is $(2 - i_{n} - 1/i_{n}) \leq 0$ for the arithmetic-mean/geometric-mean inequality. Moreover, it is $\dot{L}(S, I_{1}, \ldots, I_{5}) = 0$ only in the set $\{(S^*, I_{1}^*, \ldots, I_{5}^*)\}$. Such a property can be easily verified by simultaneously equating to zero all the polynomials in the brackets of Eq. 37, as well as of Eqs. 42, 44.

Note that if the weaker condition $\Lambda + (1 - \alpha) \Phi \ge 0$ were required instead of 5, when $\Lambda + (1 - \alpha) \Phi =0$ a unique equilibrium would still exist with $S^{*} = 0$ and $I_1^{*}, \dots, I_5^{*} > 0$, namely with the whole population infected.

---

6. Concluding remarks

The present paper presents a mathematical model developed to assess the impact of early treatment strategies in Italy. If the treatment is initiated regardless of CD4 cell counts at diagnosis, a significant impact over AIDS incidence and mortality is predicted, provided that the trend of reduced resistance to drug, observed in the decade $2003 \div 2012$, is maintained in time. Moreover, the model forecasts a significant decrease in the number of newly infected individuals, thus suggesting, in agreement with other theoretical studies (see e.g [21]), that benefits can derive from the implementation of a 'test and treat' strategy for the control of HIV infection. The model also predicts that these benefits would be enhanced if early treatments were accompanied by a decrease of the infectivity of treated subjects.

It must be recalled, however, that the knowledge about the HIV prevalence among immigrants and the infectivity of treated subjects is rather uncertain up to now. So, a better assessment of the value of these parameters would improve the quantitative reliability of the model. Conversely, more subtle phenomena like the possible dependence of some model parameters on the virus strain, as well as the possible correlation between the transmission routes and the drug resistance mutations of the transmitted virus [30], are expected to be of minor importance owing to the large prevalence of wild type HIV infections in Europe (more than 90$\%$ at 2006 [39]). Moreover, we want to stress that the proposed model was intentionally kept simple, postponing possible extensions to future investigations.

The formal analysis of the model has allowed to establish the existence of equilibria and their global stability, classifying them on the basis of the prevalence of infected among immigrants. Only in the case in which the fraction of infected among immigrants is zero, the disease-free equilibrium exists and a threshold parameter for its stability can be identified in the reproduction number. The value of the reproduction number, moreover, has a non-trivial relationship with the parameters characterizing the dynamics of the compartment related to patients under treatment.

However, the knowledge of equilibria and their stability is scarcely useful for predicting the HIV epidemic over reasonable time horizons. Numerical simulations have shown indeed a very slow convergence to the equilibrium point ($>$ 100 years), since the time from infection to death is quite long, even in untreated subjects ($\sum_{n=1}^{4} \theta_{n}^{-1} \simeq 11$ years, with the values of Table 1).

It may be questioned whether the disease progression in patients for whom therapy is interrupted or is no longer effective can be equated with the progression experienced by naive infected individuals. Thus, it would be worthy to design and study a little more complex model in which different disease progressions for naive and previously treated patients are considered. Finally, it cannot be guaranteed that the trend of reduced drug resistance will continue over the next decade. Therefore, a model of the evolution in Italy of the resistance to different classes of drugs is still needed and it is the main goal of our future work.

Acknowledgments. The present work was partially supported by the PRIN "Study of determinants of secondary and primary resistance to antiretrovirals to support control strategies of HIV transmission". The authors wish to thank the anonymous referees for their useful suggestions and comments.

Appendix. It can be easily verified that the basic hypotheses ($H_{1}$)-($H_{8}$) assumed in [18] are satisfied by model 4 with $\alpha = 0$.

Additional properties are required to prove the global asymptotic stability of the disease-free equilibrium and the global asymptotic stability of the unique endemic equilibrium. In this section, it is shown that all these additional properties are satisfied by our model.

With this aim we introduce the following notation:

• $\chi(S) = \Lambda + \Phi - \mu S$, growth rate of model 4 with $\alpha = 0$;

• $g_{n}(S, I_{n}) = \beta_{n} S I_{n}$, $n = 1, \ldots, 5$, bilinear incidence term for infections occurring by contact of susceptibles with infected in stage $n$;

• $\psi_{n}(I_{n}) = Q_{n} I_{n}$, $n = 1, \ldots, 5$, removal rate from compartments $I_{n}$, including natural death, disease progression, possibly treatment initiation/interruption or death due to the disease (see the equations of $Q_{n}$ in 11);

• $\omega_{j, n}(I_{n}) = \sigma_{j, n} I_{n}$, $1 \leq j, n \leq 5$, transfer rate from compartment $I_{n}$ to compartment $I_{j}$, where $\sigma_{2, 1} = \theta_{1}$, $\sigma_{3, 2} = \theta_{2}$, $\sigma_{4, 3} = \theta_{3}$, $\sigma_{3, 5} = \xi$ and $\sigma_{5, n} = \delta_{n}$, $n = 2, 3, 4$, while $\sigma_{j, n} = 0$ for all the remaining pairs of indexes.

We also denote by $\tilde{S}$ the number of susceptible individuals at the disease-free equilibrium (that is the value of $S$ at which the growth rate vanishes, i.e. $\chi(\tilde{S}) = 0$) and with $S^{*}$, $I_{n}^{*}$, $n = 1, \ldots, 5$, the values of the state variables at the endemic equilibrium. We also define the function $\varphi(S) = S$. On the basis of the above notation, we rewrite hypotheses ($A_{1}$), ($A_{2}$) and ($B_{1}$)-($B_{4}$) required in [18] showing that they are satisfied by our rate functions and by the chosen $\varphi(S)$:

(A1): For $n = 1, \ldots, 5$, we have

$\begin{split} &g_{n}(S, I_{n}) - g_{n}(\tilde{S}, I_{n}) = \beta_{n}(S - \tilde{S})I_{n} \leq 0, 0 \leq S \leq \tilde{S}, I_{n} \geq 0, \\ &g_{n}(S, I_{n}) - g_{n}(\tilde{S}, I_{n}) = \beta_{n}(S - \tilde{S})I_{n} = 0, I_{n} > 0, \Longrightarrow S = \tilde{S}, \end{split}$

$\displaystyle 0 < \frac{g_{n}(\tilde{S}, I_{n})}{\psi(I_{n})} = \frac{\beta_{n} \tilde{S}}{Q_{n}} < \infty, I_{n} > 0;$

(A2): For $j, n \in \{1, \ldots, 5\}$, we have

$\displaystyle 0 \leq \frac{\omega_{j, n}(I_{n})}{\psi(I_{n})} = \frac{\sigma_{j, n}}{Q_{n}} < \infty, I_{n} > 0;$

(B1): For $S \neq S^{*}$,

$(\chi(S) - \chi(S^{*}))(\varphi(S) - \varphi(S^{*})) = - \mu(S - S^{*})^{2} < 0;$

(B2): For $I_{n} > 0$, $n = 1, \ldots, 5$, and $0 \leq S \leq \tilde{S}$,

$\displaystyle \left(\frac{g_{n}(S, I_{n})}{\varphi(S)} - \frac{g_{n}(S^{*}, I_{n}^{*})}{\varphi(S^{*})}\right)\!\left(\frac{g_{n}(S, I_{n})}{\varphi(S) \psi_{n}(I_{n})} - \frac{g_{n}(S^{*}, I_{n}^{*})}{\varphi(S^{*}) \psi_{n}(I_{n}^{*})}\right)\!\! \!=\! \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \displaystyle \beta_{n} (I_{n} - I_{n}^{*}) \left(\frac{\beta_{n}}{Q_{n}} - \frac{\beta_{n}}{Q_{n}}\right)\! \!=\! 0;$

(B3): For $I_{n} > 0$, $n = 1, \ldots, 5$, and for each pair $j, n \in \{1, \ldots, 5\}$,

$\displaystyle (\omega_{j, n}(I_{n}) - \omega_{j, n}(I_{n})^{*})\!\left(\frac{\omega_{j, n}(I_{n})}{\psi_{n}(I_{n})} - \frac{\omega_{j, n}(I_{n})^{*})}{\psi_{n}(I_{n}^{*})}\right)\!\! \!=\! \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \displaystyle \sigma_{j, n} (I_{n} - I_{n}^{*}) \left(\frac{\sigma_{j, n}}{Q_{n}} - \frac{\sigma_{j, n}}{Q_{n}}\right)\! \!=\! 0;$

(B4): All the functions $g_{n}(S, I_{n})$, $\psi_{n}(I_{n})$, $\omega_{j, n}(I_{n})$, with $j, n = 1, \ldots, 5$, and $\sigma_{j, n} \neq 0$ are strictly monotone with respect to $I_{n}$.

Let us finally define the weight matrix $M$ of the infection-transfer graph $G$ of model 4 (see Fig. 12), as $M = (m_{j, n})$, where $m_{1, n} = \omega_{1, n}(I_{n}^{*}) + g_{n}(S^{*}, I_{n}^{*})$, $n = 1, \ldots, 5$, while $m_{j, n} = \omega_{j, n}(I_{n}^{*})$, $n = 1, \ldots, 5$, for $j = 2, \ldots, 5$. Substituting the expressions of the functions $g_{n}(S^{*}, I_{n}^{*})$, $\omega_{j, n}(I_{n}^{*})$, $j, n = 1, \ldots, 5$, we have

$M = \left[ \begin{array}{ccccc} \beta_{1}S^{*}I_{1}^{*}&\beta_{2}S^{*}I_{2}^{*}&\beta_{3}S^{*}I_{3}^{*}&\beta_{4}S^{*}I_{4}^{*}&\beta_{5}S^{*}I_{5}^{*} \\ \theta_{1} I_{1}^{*}&0&0&0&0\\ 0&\theta_{2} I_{2}^{*}&0&0&\xi I_{5}^{*}\\ 0&0&\theta_{3} I_{3}^{*}&0&0\\ 0&\delta_{2} I_{2}^{*}&\delta_{3} I_{3}^{*}&\delta_{4} I_{4}^{*}&0 \end{array} \right]. \label{WeightMatrix}$

(38)

As the the weighted graph ($M, G$) of the model 4 is strongly connected, the matrix $M$ results to be irreducible.

On the basis of the properties given above, Theorems 4.1 and 5.1 in [18] hold in our case, guaranteeing that the disease-free equilibrium is globally asymptotically stable for $R \leq 1$ and unstable for $R > 1$, while the endemic equilibrium is unique and also asymptotically stable when $R > 1$ ($R$ is the reproduction number defined in Section 5.1).

However, for the completeness of our study, we report here some details of the proof of Theorem 5.2 given in Section 5.1. As suggested in [18], we consider the following general Lyapunov function

$L(S, I_{1}, \ldots, I_{5}) = c_{1} \int_{S^{*}}^{S} \frac{\varphi(\tau) - \varphi(S^{*})}{\varphi(\tau)} d\tau + \sum\limits_{n = 1}^{5} c_{n} \int_{I_{n}^{*}}^{I_{n}} \frac{\psi(\tau) - \psi(I_{n}^{*})}{\psi(\tau)} d\tau. \label{GeneralLyapunovFunction}$

(39)

Recalling that in our case the removal rates $\psi_{n}(I_{n})$ are equal to $Q_{n} I_{n}$, $n = 1, \ldots, 5$, for the particular choice $\varphi(S)=S$ the general function 39 becomes the function given in 25. Then, we choose the coefficient $c_{n}$ as the cofactor of the $n$-th diagonal entry of the algebraic Laplacian matrix of $M$ defined as

$\begin{split} \mathcal{L}(M) =&\, {\rm diag}\left(\sum\limits_{n = 1}^{5} m_{1, n}, \ldots, \sum\limits_{n = 1}^{5} m_{5, n}\right) - M \\ = & \left[ \begin{array}{ccccc} \displaystyle \sum\limits_{n = 2}^{5} \beta_{n}S^{*}I_{n}^{*}&-\beta_{2}S^{*}I_{2}^{*}&-\beta_{3}S^{*}I_{3}^{*}&-\beta_{4}S^{*}I_{4}^{*}&-\beta_{5}S^{*}I_{5}^{*} \\ -\theta_{1} I_{1}^{*}&\theta_{1} I_{1}^{*}&0&0&0\\ 0&-\theta_{2} I_{2}^{*}&\theta_{2} I_{2}^{*}+\xi I_{5}^{*}&0&-\xi I_{5}^{*}\\ 0&0&-\theta_{3} I_{3}^{*}&\theta_{3} I_{3}^{*}&0\\ 0&-\delta_{2} I_{2}^{*}&-\delta_{3} I_{3}^{*}&-\delta_{4} I_{4}^{*}&\displaystyle \sum\limits_{n = 2}^{4} \delta_{n} I_{n}^{*} \end{array} \right]. \end{split}$

(40)

It can be easily verified that the five cofactors of $\mathcal{L}(M)$ are given by 29, 30. By inferring the relation $c_{2} = (\sum_{n = 2}^{5}\!B_{n}/\Theta_{1}) c_{1}$ from equations 29, the function 28 can be rewritten as

$\dot{L}(S, I_{1}, \ldots, I_{5}) = T^{\prime}(S, I_{1}, \ldots, I_{5}) + T^{\prime\prime}(S, I_{1}, \ldots, I_{5}), \label{DerivLyapunovFunction3}$

(41)

where

$T^{\prime}(S, I_{1}, \ldots, I_{5}) = \displaystyle c_{1} \left(\!(\mu S^{*} + B_{1})\!\left(\!2 - s - \frac{1}{s}\right)\! + B_{2}\!\left(\!3 - \frac{s i_{2}}{i_{1}} - \frac{1}{s} - \frac{i_{1}}{i_{2}}\right)\right), \label{DerivLyapunovFunction3_Term1}$

(42)

and

$\begin{split} &T^{\prime\prime}(S, I_{1}, \ldots, I_{5}) \\ \displaystyle =&\, c_{1}\!\sum\limits_{n = 3}^{5}\!B_{n}\!\left(\!3 - \frac{s i_{n}}{i_{1}} - \frac{1}{s} - \frac{i_{1}}{i_{2}} - i_{2} + i_{n}\!\right) \\ \displaystyle&+ c_{3}\!\left(\!\theta_{2} I_{2}^{*}\!\left(\!1 - \frac{i_{2}}{i_{3}} - i_{3} + i_{2}\right)\! + \xi I_{5}^{*}\!\left(\!1 - \frac{i_{5}}{i_{3}} - i_{3} + i_{5}\!\right)\!\!\right) \\ \displaystyle&+ c_{4}\!\left(\theta_{3} I_{3}^{*}\!\left(\!1 - \frac{i_{3}}{i_{4}} - i_{4} + i_{3}\!\right)\!\!\right) \\ \displaystyle&+ c_{5}\!\left(\!\delta_{2} I_{2}^{*}\!\left(\!\!1 - \frac{i_{2}}{i_{5}} - i_{5} + i_{2}\!\right)\!\!+\!\delta_{3} I_{3}^{*}\!\left(\!\!1 - \frac{i_{3}}{i_{5}} - i_{5} + i_{3}\!\right)\!\!+\!\delta_{4} I_{4}^{*}\!\left(\!\!1 - \frac{i_{4}}{i_{5}} - i_{5} + i_{4}\!\right)\!\!\right)\!. \end{split} \label{DerivLyapunovFunction3_Term2}$

(43)

Since from the arithmetic-mean/geometric-mean inequality we have $(2 - s - 1/s) \leq 0$ and $(3 - s i_{2}/i_{1} - 1/s - i_{1}/i_{2}) \leq 0$, the function 42 is non-positive. In order to prove the non-positiveness of the second term of the function 41, we directly substitute in Eq. 43 the expressions of $c_{1}$, $c_{3}$, $c_{4}$ and $c_{5}$ given in 29. Then, collecting the terms with respect to $B_{n}$, $n = 3, 4, 5$, we obtain

$\begin{align*} &T^{\prime\prime}(S, I_{1}, \ldots, I_{5}) \\ \displaystyle =&\, B_{3}\!\!\left(\!\!\Theta_{1}\Theta_{3}\!\!\left(\!\!\Theta_{2}\!\!\sum\limits_{n = 2}^{4}\!\Delta_{n}\!\!\left(\!\!4 - \frac{s i_{3}}{i_{1}} - \frac{1}{s} - \frac{i_{1}}{i_{2}} - \frac{i_{2}}{i_{3}} \!\right)\right.\right. \\ \displaystyle&+ \Xi_{5}\Delta_{2}\!\!\left(\!\!5 \!-\! \frac{s i_{3}}{i_{1}} - \frac{1}{s} - \frac{i_{1}}{i_{2}} - \frac{i_{5}}{i_{3}} - \frac{i_{2}}{i_{5}}\!\right) \\ \displaystyle&+\!\left.\left.\Xi_{5}\Delta_{3}\!\left(\!2 - \frac{i_{5}}{i_{3}} - \frac{i_{3}}{i_{5}}\!\right)\! + \!\Xi_{5}\Delta_{4}\!\left(\!3 - \frac{i_{5}}{i_{3}} - \frac{i_{3}}{i_{4}} - \frac{i_{4}}{i_{5}}\!\right)\!\right)\!\right)\\ \displaystyle&+ B_{4}\!\!\left(\!\!\Theta_{1}\Theta_{3}\!\!\left(\!\!\Theta_{2}\!\!\sum\limits_{n = 2}^{4}\!\Delta_{n}\!\!\left(\!\!5 - \frac{s i_{4}}{i_{1}} - \frac{1}{s} - \frac{i_{1}}{i_{2}} - \frac{i_{2}}{i_{3}} - \frac{i_{3}}{i_{4}}\!\right)\right.\right. \\ \displaystyle&+ \Xi_{5}\Delta_{2}\!\!\left(\!\!6 \!-\! \frac{s i_{4}}{i_{1}} - \frac{1}{s} - \frac{i_{1}}{i_{2}} - \frac{i_{5}}{i_{3}} - \frac{i_{3}}{i_{4}} - \frac{i_{2}}{i_{5}}\!\right) \\ \displaystyle&+\!\left.\left.\Xi_{5}\Delta_{3}\!\left(\!2 - \frac{i_{5}}{i_{3}} - \frac{i_{3}}{i_{5}}\!\right)\! + \!\Xi_{5}\Delta_{4}\!\left(\!3 - \frac{i_{5}}{i_{3}} - \frac{i_{3}}{i_{4}} - \frac{i_{4}}{i_{5}}\!\right)\!\right)\!\right)\\ \displaystyle&+ B_{5}\!\!\left(\!\!\Theta_{1}\Theta_{3}\!\!\left(\!\!\Theta_{2}\Delta_{2}\!\left(\!4 - \frac{s i_{5}}{i_{1}} - \frac{1}{s} - \frac{i_{1}}{i_{2}} - \frac{i_{2}}{i_{5}}\!\right)\right.\right. \\ \displaystyle&+ \Theta_{2}\Delta_{3}\!\left(\!5 - \frac{s i_{5}}{i_{1}} - \frac{1}{s} - \frac{i_{1}}{i_{2}} - \frac{i_{2}}{i_{3}} - \frac{i_{3}}{i_{5}}\!\right) \\ \displaystyle&+ \Theta_{2}\Delta_{4}\!\left(\!6 - \frac{s i_{5}}{i_{1}} - \frac{1}{s} - \frac{i_{1}}{i_{2}} - \frac{i_{2}}{i_{3}} - \frac{i_{3}}{i_{4}} - \frac{i_{4}}{i_{5}}\!\right) \\ \displaystyle&+ \Xi_{5}\Delta_{2}\!\left(\!4 - \frac{s i_{5}}{i_{1}} - \frac{1}{s} - \frac{i_{1}}{i_{2}} - \frac{i_{2}}{i_{5}}\!\right)\! + \! \end{align*}$

$\begin{split} \displaystyle&+\!\left.\Xi_{5}\Delta_{3}\!\left(\!2 - \frac{i_{5}}{i_{3}} - \frac{i_{3}}{i_{5}}\!\right)\! + \!\Xi_{5}\Delta_{4}\!\left(\!3 - \frac{i_{5}}{i_{3}} - \frac{i_{3}}{i_{4}} - \frac{i_{4}}{i_{5}}\!\right)\!\right). \end{split} \label{DerivLyapunovFunction3_Term2Final}$

(44)

From 44 it is easy to deduce that $T^{\prime\prime}(S, I_{1}, \ldots, I_{5}) \leq 0$ since, exploiting again the arithmetic-mean/geometric-mean inequality, all the polynomials in the brackets result to be actually non-positive. With this last statement we can conclude that $\dot{L}(S, I_{1}, \ldots, I_{5}) \leq 0$.

Finally, in order to prove the asymptotic stability of $E_{en}$, we need to verify that $\dot{L}(S, I_{1}, \ldots, I_{5}) = 0$ if and only if $(S, I_{1}, \ldots, I_{5}) \equiv (S^{*}, I_{1}^{*}, \ldots, I_{5}^{*})$. Indeed, if $\dot{L}(S, I_{1}, \ldots, I_{5}) = 0$ it follows that all the polynomials in the brackets of equations 42 and 44 are equal to zero. Then, it is easy to obtain that each point satisfying $\dot{L}(S, I_{1}, \ldots, I_{5}) = 0$ is such that

$S = S^{*}, i_{n} = \lambda, \ \ n = 1, \ldots, 5, \label{PropertiesOfPointsDotL=0}$

(45)

where $\lambda$ is a positive constant. In other words, each solution satisfying $\dot{L}(S, I_{1}, \ldots,$ $I_{5}) = 0$ is a time-invariant solution identified by $S(t) = S^{*}$, $I_{n}(t) = \lambda I_{n}^{*}$, $n = 1, \ldots, 5$. Consequently, from the equation of $\dot{S}$ in 4 we obtain

$\displaystyle 0 = \Lambda - \mu S^{*} - \lambda \sum\limits_{n = 1}^{5} \beta_{n} I_{n}^{*} S^{*}. \label{Time-InvariantEq_S}$

(46)

Since the same equation holds for the endemic equilibrium with $\lambda = 1$, we deduce that Eq. 46 can be verified only for $\lambda = 1$, because of its linearity with respect to $\lambda$. So, we can conclude $\dot{L}(S, I_{1}, \ldots, I_{5}) = 0 \, \, \Leftrightarrow \, \, (S, I_{1}, \ldots, I_{5}) \equiv (S^{*}, I_{1}^{*}, \ldots, I_{5}^{*})$.

---