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

1.   Introduction

In the last decades, mathematical models have played a very important role on the study of the analysis of the spread of infectious diseases. One of those infectious diseases, is the acquired immunodeficiency syndrome (AIDS), which is caused by the infection with the human immunodeficiency virus (HIV). HIV continues to be a major global public health issue, having claimed more than 35 million lives so far. In 2017,940 000 people died from HIV-related causes globally. There were approximately 36.9 million people living with HIV, at the end of 2017, with 1.8 million people becoming newly infected in 2017 globally ^[29]. HIV is spread all over the world, however the World Health Organization African Region is the most affected region, with 25.7 million people living with HIV in 2017. The African region also accounts for over two thirds of the global total of new HIV infections ^[29]. In this context, numerous mathematical models have been proposed for HIV/AIDS transmission and other infectious diseases, but although they can incorporate important informations about the characteristics of epidemic outbreaks, many of them do not take into account the heterogeneity of the geographical landscape (see e.g. ^{[14,15,22,23,25,26]}). However, this geographical heterogeneity, which seems to represent a key factor in understanding the spreading of infectious diseases, can be studied through the complex networks approach, which combines dynamical systems with graph theory, in order to propose refined mathematical models.

Usually, a complex network is built by considering a graph, given by a set of vertices and a set of edges, and by coupling each vertex with an instance of a given dynamical system, which can be determined by a set of differential equations. Recent works have been devoted to complex networks for studying real-world applications, like neural networks ^[6,7,9,28], behavioral networks ^[4,5] or biological networks ^[18,19]. Other works have also been dedicated to the study of epidemiological problems and their relationship with social networks ^[13,20,21]. The effect of the topology of the graph, which is determined by the disposal of its edges, on the dynamics of the network and the possible synchronization state are widely studied ^[1,2,3,10], but it is still unknown for complex networks of nonidentical systems if one can find some topology that favors a particular dynamics.

In this paper, we propose an original model of complex network of nonidentical dynamical systems, in order to analyze the spread of epidemics within an heterogeneous geographical zone. To this end, we consider a recent HIV/AIDS model, proposed in ^[26], given by a system of ordinary differential equations, for which a basic reproduction number $R_0$ can be determined. It has been proved in ^[27] that this refined model admits a disease-free equilibrium (DFE), which is globally asymptotically stable if $R_0 < 1$, and an endemic equilibrium (EE), which is globally asymptotically stable if $R_0 > 1$. Thus we construct a complex network by coupling nonidentical instances of the HIV/AIDS system proposed in ^[26], that is, multiple instances of the epidemiological system with distinct parameters values. We construct the complex network so that it can take into account the situation where a part of the population is not concerned with the migrations. Moreover, we assume that the displacements can be different in some places of the network, that migrations are instantaneous, and that individuals are not subject to an evolution from one compartment to another during migration from one node to another. Up to our knowledge, these assumptions are a novelty on the mathematical modeling of the spread of HIV/AIDS epidemics. We consider a graph $\mathscr{G} = (\mathscr{V}, \, \mathscr{E})$, where the set of vertices $\mathscr{V}$ models the zones of high population concentration, and the set of edges $\mathscr{E}$ corresponds to human displacements among those separated zones. We focus on the situation where the set of vertices of the graph $\mathscr{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.

2.1.   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 $\Lambda$, $\mu$, $\beta$, $\eta_C$, $\eta_A$, $\phi$, $\rho$, $\alpha$, $\omega$ and $d$ take fixed values and their significance is detailed in Table 1.

We recall that system (2.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 $\Sigma_0$ given by (2.2) is globally asymptotically stable for $R_0 < 1$.

We also recall that model (2.1) has a unique endemic equilibrium $\Sigma_+ = (S^*, I^*, C^*, A^*)$, which is globally asymptotically stable, whenever $R_0 > 1$.

Lemma 1. ^[27] The model (2.1) has a unique endemic equilibrium $\Sigma_+ = (S^*, I^*, C^*, A^*)$ whenever $R_0 > 1$, which is given by

Theorem 2. ^[27] The endemic equilibrium $\Sigma_+$ given by (1) is globally asymptotically stable for $R_0 > 1$.

The HIV/AIDS model (2.1) can be rewritten in the following way

where $x = (S, \, I, \, C, \, A)^T$, $p = (\Lambda, \, \beta, \, \eta_C, \, \eta_A, \, \mu, \, \rho, \, \phi, \, \omega, \, \alpha, \, d)$ and $f(x, \, p)$ is defined by:

Nonidentical instances of system (2.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.

2.2.   Construction of the complex network

Let us consider a graph $\mathscr{G} = (\mathscr{V}, \, \mathscr{E})$ made of a finite set $\mathscr{V}$ of $n$ vertices, where $n$ denotes an integer greater than $2$, and a finite set $\mathscr{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 $\mathscr{V}$ can be split into at least two subsets of vertices $\mathscr{V}_1$ and $\mathscr{V}_2$. We couple the vertices of $\mathscr{V}_1$ with an instance of system (2.4) for which $R_0 < 1$, and the vertices of $\mathscr{V}_2$ with an instance of system (2.4) for which $R_0 > 1$. The complex network is determined by the following non-linear and autonomous differential system:

where

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) \in \mathscr{E}$, $k \neq j$, we have $L_{j, k} > 0$. If $(k, \, j) \notin \mathscr{E}$, $k \neq 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 $\varepsilon_S$, $\varepsilon_I$, $\varepsilon_C$ and $\varepsilon_A$.

In this complex network model, we consider that an edge $(k, \, j) \in \mathscr{E}$, $k \neq j$, models a connection between two vertices $k$ and $j$, which corresponds to human displacements from vertex $k$ towards vertex $j$. Moreover, the parameter $\varepsilon_S$ models the rate of susceptible individuals on vertex $k$ which migrate towards vertex $j$. The parameters $\varepsilon_I$, $\varepsilon_C$ and $\varepsilon_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 \in \lbrace 1, \, \dots, \, n \rbrace$:

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 \in \lbrace 1, \, \dots, \, n \rbrace$.

3.   Positively invariant region

In this section, we prove that the complex network problem (2.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 (2.5), which corresponds to the fact that the matrix of connectivity $L$ is a zero column sum matrix.

3.1.   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 $\mathbb{R}$, and $t_0 \in \mathbb{R}$. We assume that $g(t) \geq 0$ for all $t \in \mathbb{R}$, and that $\zeta(t_0) \geq 0$. Let $\zeta(t)$ be a solution of (3.1) defined on $[t_0, \, t_0+\tau]$ with $\tau > 0$, such that $\zeta(t_0) \geq 0$.

Then we have $\zeta(t) \geq 0$, for all $t \in [t_0, \, t_0+\tau]$.

Lemma 3. Let $\zeta$ be a continuous function defined on $[0, \, T]$, with $T > 0$, continuously differentiable on $]0, \, T]$, and satisfying the differential inequality

with two positive coefficients $\delta_1$, $\delta_2$.

Then we have

3.2.   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 (2.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 \in (\mathbb{R}^+)^{4n}$, the Cauchy problem

where $F$, $P$, $L$ and $H$ are defined as above (see Section 2.2), admits a unique solution defined on $[0, \, T]$ with $T > 0$, whose components are non-negative on $[0, \, T]$.

Proof. Let us consider an initial condition $X_0 \in (\mathbb{R}^+)^{4n}$. We denote by $X(t, \, X_0)$ the solution of the Cauchy problem (3.2), defined on $[0, \, T]$ with $T > 0$.

We introduce an auxiliary problem defined by

for each $j \in \lbrace 1, \, \dots, \, n \rbrace$, where the coefficient $\gamma_j$ corresponds to the fluxes exiting from node $j$, and is given by

Let us denote by $\tilde{X}(t, \, X_0)$ the solution of the auxiliary problem (3.3), stemming from the same initial condition $X_0 \in (\mathbb{R}^+)^{4n}$, defined on $[0, \, \tilde{T}]$.

Applying Lemma 2, we easily prove that the components of $\tilde{X}(t, \, X_0)$ are non-negative. This implies that $\tilde{X}(t, \, X_0)$ is also a solution of the Cauchy problem (3.2) on $[0, \, \tilde{T}]$. By uniqueness, we have $\tilde{X}(t, \, X_0) = X(t, \, X_0)$ for all $t \in [0, \, T] \cap [0, \, \tilde{T}]$. Finally, it is seen that $T = \tilde{T}$, which achieves the proof.

3.3.   Boundedness of the solutions of the complex network problem

Let us introduce the minimum mortality rate $\mu_0$ defined by

the positive coefficient $\Lambda_0$ defined by

and the compact region

The total population in the complex network, defined by

satisfies

since the matrix of connectivity $L$ is a zero column sum matrix. Applying Lemma 3 leads to

thus we obtain the following theorem.

Theorem 4. The region $\Omega$ defined by (3.4) is positively invariant under the flow induced by the complex network (2.5).

Remark 1. It is easily seen that the positively invariant region $\Omega$ for the complex network problem (2.5) satisfies

where $\Omega_i = \left\lbrace (x_j)_{1\leq j \leq 4} \in (\mathbb{R}^+)^{4}\; ;\; \sum\limits_{j = 1}^{4} x_j \leq \frac{\Lambda_i}{\mu_i} \right\rbrace$ corresponds to the positively invariant region of the node $(i)$ in absence of coupling. Roughly speaking, the couplings can enlarge the phase space of the flow induced by the network problem.

4.   Stability analysis of the complex network

In this section, we explore the effect of the couplings on the dynamics of the complex network (2.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.

4.1.   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.

Roughly speaking, the coupling coefficient $L_{2, 1}$ acts on vertex $(1)$ as if the mortality rate $\mu_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 (4.1) admits a disease-free equilibrium point $\Sigma^0$ given by

Following the method proposed in ^[8] for computing the basic reproduction number, we have the following result.

Proposition 1. The basic reproduction number of model (4.1) is

where

and

with

and

Proof. Let $\mathcal{F}_i(t)$ be the rate at which new infections appear in the $i$-th compartment and $\mathcal{V}_i^+(t)$ be the "individuals" transfer rate into the $i$-th compartment in all other ways. Similarly, let $\mathcal{V}_i^-(t)$ denote the "individuals" transfer rate out of the $i$-th compartment, for which

Therefore, we take

and

The Jacobian matrices $F$ of $\mathcal{F}(t)$ and $V$ of $\mathcal{V}(t)$ are given by

where

and

Evaluating the matrices $F$ and $V$ at the disease-free equilibrium $\Sigma^0$ given by (4.2), we find

with

The eigenvalues of the matrix $F_0 V_0^{-1}$ are given by:

The basic reproduction number is given by the dominant eigenvalue of the matrix $F_0 V_0^{-1}$, that is, $R_0$ takes the value given by (4.3).

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 $\varepsilon_S = \varepsilon_I = \varepsilon_C = \varepsilon_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).

We easily prove that $R_{0, 2}$ is an increasing function of $L_{2, 1}$. Indeed, we have

with $k = \frac{\mathcal{N}_2}{\mathcal{D}_2 \mu_2}$, $d_1 = \varepsilon_S (\Lambda_1 + \Lambda_2)$, $d_2 = \Lambda_2 \mu_1$, $d_3 = \varepsilon_S$ and $d_4 = \mu_1$. Since $d_1 d_4 - d_2 d_3 = \mu_1 \varepsilon_S \Lambda_1 > 0$, we can conclude that $R_{0, 2}$ is an increasing function of $L_{2, 1}$. Figure 2(a) illustrates this increasing shape of $R_{0, 2}$ with respect to $L_{2, 1}$ for the following parameters values:

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 $\varepsilon_S$, $\varepsilon_I$, $\varepsilon_C$ and $\varepsilon_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 $\varepsilon_S$, $\varepsilon_I$, $\varepsilon_C$ and $\varepsilon_A$. It seems that $R_{0, 1}$ and $R_{0, 2}$ are robust to a variation of the coupling strengths $\varepsilon_I$ and $\varepsilon_A$, whereas a variation of $\varepsilon_S$ or $\varepsilon_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 (4.1). Moreover, the coupling strengths $\varepsilon_S$ and $\varepsilon_C$ seem to play antagonistic roles, since an increase of $\varepsilon_S$ provokes an increase of $R_{0, 2}$, whereas an increase of $\varepsilon_C$ provokes an increase of $R_{0, 1}$.

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 Figure 4). However, the output is unreadable, even with relevant simplifications of the parameters of the system.

In the same manner, it is possible to obtain the complete symbolical expression of the disease-free equilibrium $\Sigma^0$ of a three-nodes chain (see Figure 5):

where

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.

4.2.   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 (2.5) admits a unique stable disease-free equilibrium under reasonable assumption.

Theorem 5. The complex network (2.5) admits a unique disease-free equilibrium $\Sigma_0$, which is globally asymptotically stable in the region $\Omega$ defined by (3.4), provided

for all $i \in \lbrace 1, \, \dots, \, n \rbrace$, where $\mathcal{N}_i$ and $\mathcal{D}_i$ are defined by

with

Proof. The equilibrium points of the network problem (2.5) are the solutions of the system

We determine the disease-free equilibria by assuming that $I_j = C_j = A_j = 0$ for all $j \in \lbrace 1, \, \dots, \, n \rbrace$. Then we directly obtain

and simultaneously

The latter system is a linear system which can be written

with $B = B_1 - \varepsilon_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 = \text{diag}\left\{{\mu_1, \, \dots, \, \mu_n}\right\}$. $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 (4.9) admits a unique solution, which corresponds to the unique disease-free equilibrium $\Sigma_0$ of the network problem (2.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 (4.8). We compute the orbital derivative $\dot{V}$ of the Lyapunov functional $V$ along a solution $X$ starting in $\Omega$ :

which guarantees that $\dot{V} \leq 0$, since we assume that $\frac{\Lambda_0}{\mu_0} \frac{\mathcal{N}_i}{\mathcal{D}_i} < 1$ for all $i \in \lbrace 1, \, \dots, \, n \rbrace$.

Finally, it is seen that $\dot{V} = 0$ if and only if $I_i = C_i = A_i = 0$ for all $i \in \lbrace 1, \, \dots, \, n \rbrace$. The conclusion follows from LaSalle invariance principle ^[16].

Remark 4. Since we have $\Lambda_i \leq \Lambda_0$ and $\mu_i \geq \mu_0$ for all $i \in \lbrace 1, \, \dots, \, n \rbrace$, assumption (4.7) implies that

As it is relevant to introduce again $R_{0, i} = \dfrac{\Lambda_i}{\mu_i} \dfrac{\mathcal{N}_i}{\mathcal{D}_i}$ for each $i \in \lbrace 1, \, \dots, \, n \rbrace$, it is seen that assumption (4.7) 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.

5.1.   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 (Table 2) (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 \simeq 0.914$ for group $1$, $R_0 \simeq 1.371$ for group $2$ and $R_0 \simeq 7.312$ for group $3$.

Remark 5. The value of the basic reproduction number for group $3$ implies that assumption (4.7) 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:

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 \leq j \leq 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.

5.2.   Randomly generated topologies

The numerical integration on a finite time interval $[0, \, T]$ of the complex network (2.5) modeling Cape Verde archipelago has been performed using the $\textsf{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 \simeq 9112.77$ with $T = 200$, whereas the complete graph topology (see Figure 6(b)) leads to $L_f \simeq 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 \leq \left\vert{\mathscr{E}}\right\vert \leq 72$, and a random subset of $\left\vert{\mathscr{E}}\right\vert$ 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 $\left(L_f, \, \left\vert{\mathscr{E}}\right\vert\right)$. The green dotted vertical line corresponds to the level of infected individuals for an empty topology.

We observe that the final level of infected individuals $L_f$ varies a lot with respect to the number $\left\vert{\mathscr{E}}\right\vert$ 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 \simeq 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.

5.3.   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 first weakly dense topology we aim to analyze is a near-optimal topology detected by the random simulation (see Figure 6(c)); it admits a set of $14$ edges, given by

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 \simeq 9085.09$.

The second weakly dense topology (see Figure 6(d)) we focus on is another near-optimal topology; it admits a set of $9$ edges is given by

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 \simeq 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.

6.   Conclusion

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.

Acknowledgments

The authors are very grateful to the anonymous reviewers whose comments greatly improved the presentation of the paper.

This research was partially supported by the Portuguese Foundation for Science and Technology (FCT) within projects UID/MAT/04106/2019 (CIDMA) and PTDC/EEI-AUT/2933/2014 (TOCCATTA), co-funded by FEDER funds through COMPETE2020 – Programa Operacional Competitividade e Internacionalização (POCI) and by national funds (FCT). Silva is also supported by national funds (OE), through FCT, I.P., in the scope of the framework contract foreseen in the numbers 4, 5 and 6 of the article 23, of the Decree-Law 57/2016, of August 29, changed by Law 57/2017, of July 19.

Conflict of interest

The authors declare there is no conflict of interest in this paper.