Influence of the topology on the dynamics of a complex network of HIV/AIDS epidemic models

In this paper, we propose an original complex network model for an epidemic problem in an heterogeneous geographical area. The complex network is constructed by coupling nonidentical instances of a HIV/AIDS epidemiological model for which a disease-free equilibrium and an endemic equilibrium can coexist. After proving the existence of a positively invariant region for the solutions of the complex network problem, we investigate the effect of the coupling on the dynamics of the network, and establish the existence of a unique disease-free equilibrium for the whole network, which is globally asymptotically stable. We prove the existence of an optimal topology that minimizes the level of infected individuals, and apply the theoretical results to the case of the Cape Verde archipelago.

separated zones. We focus on the situation where the set of vertices of the graph G is split into at least two subsets: the first subset being coupled with instances of the HIV model for which the basic reproduction number satisfies R 0 < 1, thus admitting a unique equilibrium which is a DFE; and the second subset being coupled with instances of the HIV model for which R 0 > 1, thus presenting the coexistence of a DFE and an EE. We address the following questions. Does the coupling between the two subsets of vertices create new equilibrium points? Is it possible to eliminate the endemic equilibriums with a suitable disposition of the couplings? If not, is it possible to minimize the propagation of the epidemic in the network by searching an optimal topology?
Those questions are of great interest, and still have not been deeply studied. Recently, it has been proved, in the particular case of a behavioral model [5], that oriented chains play an important role for reaching an expected equilibrium. Here, the equilibrium we aim to favor in the network is the DFE. The main goal of this paper is to find a topology for which the DFEs of a given subset of the network drive the whole network to a global DFE.
This paper is organized as follows. In Section §2. , we recall the equations and stability results of the HIV/AIDS model proposed in [26] and analyzed in [27] that is considered in our study, and we improve the previous model by constructing a novel complex network model. In Section §3. , we show that the complex network model is well-posed, by proving the existence of a positively invariant region that guarantees the non-negativity of the solutions, together with their boundedness and global existence. In Section §4. , we explore the influence of the coupling on the equilibrium points of the network, by computing the global basic reproduction number for particular small networks, and prove under reasonable assumptions that the network admits a unique DFE which is globally asymptotically stable. A case study is considered in Section §5. , where the HIV/AIDS epidemic in the Cape Verde archipelago is analyzed and relevant numerical simulations are presented. Moreover, the existence of an optimal topology minimizing the level of infection in the network is investigated. We end our paper with Section §6. with conclusions and discussion of possible future works. §2. Problem statement In this section, we recall the HIV/AIDS compartmental model, given by a system of differential equations firstly proposed in [26], which describes the transmission dynamics of HIV in a homogeneous population with variable size. In Subsection 2.2. we explicit the construction of a complex network determined with nonidentical instances of the previous HIV/AIDS model.

HIV/AIDS model
We consider a human population affected by a HIV/AIDS epidemics, from [26]. The total population is divided into four mutually-exclusive compartments: susceptible individuals (S); HIV-infected individuals with no clinical symptoms of AIDS (the virus is living or developing in the individuals but without producing symptoms or only mild ones) but able to transmit HIV to other individuals (I); HIV-infected individuals under ART treatment (the so called chronic stage) with a viral load remaining low (C); and HIV-infected individuals with AIDS clinical symptoms (A). The total population at time t, denoted by N (t), is given by The SICA model is given by a system of four ordinary differential equations that can be written as follows (see [27]): The parameters Λ, µ, β, η C , η A , φ, ρ, α, ω and d take fixed values and their significance is detailed in Table 1. We recall that system (1) admits a disease-free equilibrium (DFE) given by We introduce the basic reproduction number R 0 , which represents the expected average number of new HIV infections produced by a single HIV-infected individual when in contact with a completely susceptible population, given by where Theorem 1. [27] The disease free equilibrium Σ 0 given by (2) is globally asymptotically stable for R 0 < 1.
The HIV/AIDS model (1) can be rewritten in the following waẏ Nonidentical instances of system (4) can be coupled with the vertices of a graph in order to give rise to a complex network, as we are going to see in the coming subsection.

Construction of the complex network
Let us consider a graph G = (V , E ) made of a finite set V of n vertices, where n denotes an integer greater than 2, and a finite set E of m edges, where m denotes a positive integer. This graph models the geographical zone which is affected by the epidemics. We assume that V can be split into at least two subsets of vertices V 1 and V 2 . We couple the vertices of V 1 with an instance of system (4) for which R 0 < 1, and the vertices of V 2 with an instance of system (4) for which R 0 > 1. The complex network is determined by the following non-linear and autonomous differential system: where X = (x 1 , . . . , x n ) T ∈ R 4 n , HX = (Hx 1 , . . . , Hx n ) T ∈ R 4 n , P = (p 1 , . . . , p n ) ∈ R 10 n , and F determines the internal dynamic of each vertex: Furthermore, L is the matrix of connectivity, which is defined as follows. For each edge (k, j) ∈ E , k = j, we have L j,k > 0. If (k, j) / ∈ E , k = j, we set L j,k = 0. The diagonal coefficients satisfy Finally, H is the matrix of the coupling strengths and it is given by with non negative coefficients ε S , ε I , ε C and ε A .
In this complex network model, we consider that an edge (k, j) ∈ E , k = j, models a connection between two vertices k and j, which corresponds to human displacements from vertex k towards vertex j. Moreover, the parameter ε S models the rate of susceptible individuals on vertex k which migrate towards vertex j. The parameters ε I , ε C and ε A are defined analogously for the compartments I, C and A respectively. This implies that our model can take into account the situation where a part of the population is not concerned with the migrations. Additionally, each connection (k, j) is weighted by a positive coefficient L j,k , which means that the displacements can be different in some places of the network. It is worth emphasizing that the migrations are assumed to be instantaneous, and that individuals are not subject to an evolution from one compartment to another during migration from one node to another.
Next we explicit the equations which describe the state of vertex j ∈ {1, . . . , n}: where the time dependence is omitted, in order to lighten the notations. The coupling terms can be divided into fluxes exiting from vertex j and fluxes entering in vertex j, that is for all j ∈ {1, . . . , n}. §3.

Positively invariant region
In this section, we prove that the complex network problem (5) is well posed, and admits a positively invariant region. This property is obtained as a consequence of the conservation of the couplings terms in system (5), which corresponds to the fact that the matrix of connectivity L is a zero column sum matrix.

Preliminary results
We begin recalling two preliminary results. Since their proofs are well-known, we omit it.

Lemma 2. Let us consider the Cauchy problem
where f and g are two continuous functions defined on R, and t 0 ∈ R. We assume that g(t) ≥ 0 for all t ∈ R, and that ζ(

Lemma 3. Let ζ be a continuous function defined on
, and satisfying the differential inequalitẏ Then we have

Non-negativity of the solutions of the complex network problem
The next theorem guarantees the non-negativity of the solutions of the complex network problem (5), which is an obvious property to be satisfied for population dynamics models. Since the proof uses classical techniques [11,17], we only give the main steps.
Theorem 3. For any initial condition X 0 ∈ (R + ) 4n , the Cauchy problem where F , P , L and H are defined as above (see Section 2.

), admits a unique solution defined on
Proof. Let us consider an initial condition X 0 ∈ (R + ) 4n . We denote by X(t, X 0 ) the solution of the Cauchy problem (8), defined on [0, T ] with T > 0. for each j ∈ {1, . . . , n}, where the coefficient γ j corresponds to the fluxes exiting from node j, and is given by Let us denote byX(t, X 0 ) the solution of the auxiliary problem (??), stemming from the same initial condition X 0 ∈ (R + ) 4n , defined on [0,T ].
Applying Lemma 2, we easily prove that the components ofX(t, X 0 ) are non-negative. This implies thatX(t, X 0 ) is also a solution of the Cauchy problem (8) on [0,T ]. By uniqueness, we haveX(t, X 0 ) = X(t, X 0 ) for all t ∈ [0, T ] ∩ [0,T ]. Finally, it is seen that T =T , which achieves the proof.

Boundedness of the solutions of the complex network problem
Let us introduce the minimum mortality rate µ 0 defined by and the compact region The total population in the complex network, defined by since the matrix of connectivity L is a zero column sum matrix. Applying Lemma 3 leads to thus we obtain the following theorem.

Remark 1.
It is easily seen that the positively invariant region Ω for the complex network problem (5) satisfies where In this section, we explore the effect of the couplings on the dynamics of the complex network (5). We use symbolic computational methods, in the case of small networks. Furthermore, we prove the existence of a unique disease-free equilibrium which is globally asymptotically stable.

Asymmetric two-nodes network
Let us consider a two-nodes networks, with one vertex (1) for which R 0 < 1, another vertex (2) for which R 0 > 1, and a directed connection from vertex (1) towards vertex (2) (see Figure 1). In that case, the matrix of connectivity is given by and therefore the equations of the network read where we omit the dependence in t in order to lighten our notations.  (1) is associated with an instance of system (1) for which the basic reproduction number R 0 satisfies R 0 < 1, whereas the red node (2) is coupled with an instance of system (1) for which R 0 > 1.
Roughly speaking, the coupling coefficient L 2,1 acts on vertex (1) as if the mortality rate µ 1 increases, which changes the value of the basic reproduction number on vertex (1). Thus the following question arises. How does R 0 vary when L 2,1 increases? Proposition 1 below partly answers this question.
First, we easily prove that system (10) admits a disease-free equilibrium point Σ 0 given by Following the method proposed in [8] for computing the basic reproduction number, we have the following result. (10) is

Proposition 1. The basic reproduction number of model
where and with and Proof. Let F i (t) be the rate at which new infections appear in the i-th compartment and V + i (t) be the "individuals" transfer rate into the i-th compartment in all other ways. Similarly, let V − i (t) denote the "individuals" transfer rate out of the i-th compartment, for which Therefore, we take The Jacobian matrices F of F(t) and V of V(t) are given by Evaluating the matrices F and V at the disease-free equilibrium Σ 0 given by (11), we find The eigenvalues of the matrix F 0 V −1 0 are given by: The basic reproduction number is given by the dominant eigenvalue of the matrix F 0 V −1 0 , that is, R 0 takes the value given by (12).

Remark 2.
We emphasize that R 0,1 and R 0,2 correspond to the basic reproduction numbers of nodes (1) and (2) respectively, in absence of coupling (that is ε S = ε I = ε C = ε A = 0). Thus, the expression R 0 = max(R 0,1 , R 0,2 ) implies that if R 0,1 > 1 or R 0,2 > 1, then R 0 > 1. In other words, the node admitting a basic reproduction number R 0 > 1 drives the other node to a global Endemic Equilibrium (EE). It is a work in progress to generalize this pattern to more general topologies (e.g. chain networks). However, one should not conclude for the general case that the good solution is to "cut" the couplings (there may exist an optimal coupling topology which globally tempers the level of infected individuals, as we are going to show in Section §5. below).
At the opposite, one can find parameters values for which R 0,1 is a decreasing function of L 2,1 , but also other parameters values for which R 0,1 is an increasing function of L 2,1 in a neighborhood of 0. Figure 2(b) presents an example for which R 0,1 admits a maximum with respect to L 2,1 ; this example has been obtained for the following parameters values: In parallel, the coupling strengths ε S , ε I , ε C and ε A , stored in matrix H (see section 2.2. ), are also observed to play an important role (see Figure 3). For the same set of parameters as above, and a frozen coefficient L 2,1 = 0.15, we have computed the values of the basic reproduction numbers R 0,1 and R 0,2 with respect to a variation of ε S , ε I , ε C and ε A . It seems that R 0,1 and R 0,2 are robust to a variation of the coupling strengths ε I and ε A , whereas a variation of ε S or ε C can induce an important variation in R 0,1 and R 0,2 , which can imply a change in the dynamics of both nodes in the complex network (10). Moreover, the coupling strengths ε S and ε C seem to play antagonistic roles, since an increase of ε S provokes an increase of R 0,2 , whereas an increase of ε C provokes an increase of R 0,1 . Figure 4). However, the output is unreadable, even with relevant simplifications of the parameters of the system.

Remark 3. After tedious symbolical computations, it is possible to obtain the expressions of the basic reproduction numbers in the case of a symmetric two-nodes network (see
In the same manner, it is possible to obtain the complete symbolical expression of the disease-free equilibrium Σ 0 of a three-nodes chain (see Figure 5): Similarly, the global basic reproduction number of a three-nodes chain reveals that the couplings coefficients L 2,1 and L 3,2 affect the dynamics of each node of the network, and are likely to produce undesirable phenomenon. But its nebulous expression seems to forbid any relevant interpretation. However, we are going to see in the next subsection, that the existence of a unique stable disease-free equilibrium for the network is guaranteed under reasonable assumption.

Disease-Free Equilibrium of the complex network
Here, our aim is to overcome the computational difficulties met in the previous subsections. Thus we establish in the general case that the complex network (5) admits a unique stable disease-free equilibrium under reasonable assumption.
Theorem 5. The complex network (5) admits a unique disease-free equilibrium Σ 0 , which is globally asymptotically stable in the region Ω defined by (9), provided for all i ∈ {1, . . . , n}, where N i and D i are defined by Proof. The equilibrium points of the network problem (5) are the solutions of the systeṁ We determine the disease-free equilibria by assuming that I j = C j = A j = 0 for all j ∈ {1, . . . , n}. Then we directly obtainİ The latter system is a linear system which can be written with B = B 1 − ε S L, L being the matrix of connectivity defined as in Section 2.2. , and B 1 is a diagonal matrix storing the mortality rates, that is B 1 = diag {µ 1 , . . . , µ n }. L being a zero column-sum matrix, it follows that B is a strictly diagonally dominant matrix. By virtue of Levy-Desplanques Theorem [12], B is an invertible matrix. Hence, system (18) admits a unique solution, which corresponds to the unique disease-free equilibrium Σ 0 of the network problem (5). Next, we introduce the Lyapunov functional V defined by where V i is the Lyapunov function introduced in [27] (proof of Theorem 1), given by where the coefficients k 1,i , k 2,i and k 3,i are determined by (17). We compute the orbital derivativeV of the Lyapunov functional V along a solution X starting in Ω : which guarantees thatV ≤ 0, since we assume that Λ0 µ0 Ni Di < 1 for all i ∈ {1, . . . , n}. Finally, it is seen thatV = 0 if and only if I i = C i = A i = 0 for all i ∈ {1, . . . , n}. The conclusion follows from LaSalle invariance principle [16].

Remark 4.
Since we have Λ i ≤ Λ 0 and µ i ≥ µ 0 for all i ∈ {1, . . . , n}, assumption (16) implies that As it is relevant to introduce again . . , n}, it is seen that assumption (16) is a sufficient condition for the existence of a unique stable disease-free equilibrium in the network, which requires that every node in the network has a "small" basic reproduction number R 0,i . If only one node violates this condition, then the network is likely to exhibit undesirable equilibrium states. In other words, Theorem 5 generalizes the pattern discovered with a two-nodes network in Proposition 1.

§5. A case study: Cape Verde archipelago
In this section, we study the case of Cape Verde archipelago, which has been affected by HIV/AIDS epidemics for several decades. Our aim is to determine a topology which could temper the spreading of the epidemics.

Geographical background
Cape Verde is an archipelago of 10 volcanic islands, located in the Atlantic Ocean, at about 570 kilometers from the Northwest African coast. Since 1 of those 10 islands has no inhabitants, we propose to model this archipelago with a 9 nodes network (see Figure 6 below). We assume that the network is divided into 3 groups of nodes: group 1 is composed with nodes 1, 2, 3, 4, 5, group 2 with nodes 6, 7, 8, and group 3 with single node 9, corresponding to Santiago island which is the most important island in the archipelago, with the greatest number of HIV infected inhabitants. The parameters values are given in Table 3. In absence of coupling, it is relevant to compute the basic reproduction number R 0 for each group: R 0 0.914 for group 1, R 0 1.371 for group 2 and R 0 7.312 for group 3.

Remark 5.
The value of the basic reproduction number for group 3 implies that assumption (16) of Theorem 5 may not be fulfilled, which could lead to the emergence of undesirable equilibrium states, with a persistence of the infection within the population for instance. Thus it appears crucial to limit the spreading of the infection at a reasonable level, by finding a suitable topology of the network.
The coupling strengths are fixed as follows: ε S = 0.02, ε I = 0.01, ε C = 0.01, ε A = 0.01, in the case of weak coupling, or in the case of strong coupling. The initial condition X 0 partially corresponds to official data: approximate values of the total population N j (0) for each node (1 ≤ j ≤ 9) in 2015 can be found in [24], as well as approximate values of infected individuals I j (0). The values of C j (0) and A j (0) have been assumed, so that the corresponding subpopulations are in proportionality with the total population.

Randomly generated topologies
The numerical integration on a finite time interval [0, T ] of the complex network (5) modeling Cape Verde archipelago has been performed using the python language, in a GNU/LINUX environment. For each set of parameters, let us introduce the final level of infected individuals, given by   In absence of coupling (see Figure 6(a)), we obtain L f 9112.77 with T = 200, whereas the complete graph topology (see Figure 6(b)) leads to L f 9161.02. Since the couplings are likely to produce emerging equilibria, we propose to explore the possible topologies for the complex network modeling Cape Verde archipelago. The set of possible topologies being finite, there obviously exists an optimal topology minimizing the level of infection L f . Thus our goal is to determine a near-optimal topology. However, it is easily seen that a 9 nodes network can admit at most 72 edges, assuming that there are no loops nor parallel edges. The total number of possible topologies is given by the sum of binomial coefficients thus it is not reasonable to explore the total set of topologies. We propose to investigate a sample of randomly generated topologies, by choosing a random number of edges 1 ≤ |E | ≤ 72, and a random subset of |E | edges. We have computed the final level L f of infected individuals for a sample of 1400 randomly generated topologies. The result is depicted in Figure 7, where each red cross has coordinates  We observe that the final level of infected individuals L f varies a lot with respect to the number |E | of edges. It seems that a dense topology, with a number of edges neighbor to the maximal number 72, corresponding to the complete graph topology, produces a high level of infection. Meanwhile, a weakly dense topology is not a warranty for a low final level of infection. However, this random simulation has detected an optimal topology (marked with a green circle in Figure 7) for which the final level of infection is lesser than the benchmark L f 9113 obtained for an empty topology. Furthermore, we observe that the two clouds of points obtained for weak or strong coupling roughly admit similar shapes. In other words, the topology seems to be more important than the coupling strength.

Weakly dense topologies
The random simulation presented in the previous section seems to exclude dense topologies. The question of how to select a weakly dense topology, in order to temper the final level of infected individuals L f remains delicate. Finally, we present the times series corresponding to two weakly dense topologies.
The time series of the corresponding complex network are shown in Figure 8. The final level of infected individuals for that optimal topology is L f 9085.09.
The time series of the corresponding complex network are shown in Figure 9. The final level of infected individuals for that second weakly dense topology is L f 9087.50.

Remark 6.
The numerical results presented in this section can help finding a favorable situation for limiting the level of infection in Cape Verde archipelago. Indeed, the results of the simulation of randomly generated topologies appear to exclude dense topologies, which means that one should avoid important human migrations from one island to another. In the mean time, weakly dense topologies (c) and (d) presented in Figure 6 seem to favor migrations stemming from islands admitting a small basic reproduction number. Nevertheless, those interpretations should be prudently nuanced, since they are the result of a mathematical model whose scope is necessarily limited. In this work, we presented the analysis of a complex network of dynamical systems for the study of the spread of HIV/AIDS epidemics. Built with nonidentical instances of a compartmental model for which a disease-free equilibrium and an endemic equilibrium can coexist, this complex network exhibits a positively invariant region and presents a unique disease-free equilibrium which is globally asymptotically stable, under the assumption that each node composing the network admits a small basic reproduction number. However, emerging equilibria are likely to appear if this assumption is not fulfilled, and we proposed a numerical strategy in order to detect a near-optimal topology for which the level of infection is minimized. This method has been applied to the case of the Cape Verde archipelago, and we exhibited a near-optimal topology which seems to be robust with respect to a variation of the coupling strength. However, it seems delicate to identify the characteristic features of such a near-optimal topology, since weakly dense topologies can produce a high level of infection as well as limit the infection at a low level. In a future work, we aim to deepen this subtle question, which could lead to establishing a necessary and sufficient condition of synchronization in the network. In parallel, we propose to improve our model by applying an optimal control process, in order to reach a global disease-free equilibrium, in spite of the risk that a small group of nodes in the network could admit a high basic reproduction number. This control process could be introduced at a double scale, with control actions exerted into the dynamics of each node, and simultaneously control actions exerted along the connections of the network.
As a final perspective, we also intend to study of the effect of introducing delays in the migrations supported by the connections of the network, since it is likely to reveal new dynamics which might be hidden at this stage.