A network model for the propagation of Hepatitis C with HIV co-infection

TL;DR

This paper models the spread of Hepatitis C with HIV co-infection across interconnected urban networks, capturing both local interactions and inter-city travel to understand epidemic dynamics.

Contribution

It introduces a novel network-of-networks model that encodes city-level interactions and travel, providing insights into epidemic spread across urban areas.

Findings

Propagation speed depends on initial infection location.

Network structure influences epidemic extent.

Travel between cities accelerates spread.

Abstract

We define and examine a model of epidemic propagation for a virus such as Hepatitis C (with HIV co-infection) on a network of networks, namely the network of French urban areas. One network level is that of the individual interactions inside each urban area. The second level is that of the areas themselves, linked by individuals travelling between these areas and potentially helping the epidemic spread from one city to another. We choose to encode the second level of the network as extra, special nodes in the first level. We observe that such an encoding leads to sensible results in terms of the extent and speed of propagation of an epidemic, depending on its source point.

Tables

Table 1. Table 1: Annual origin-destination matrix.

	Destination
	$m_p=0.1\%$				$m_p=1\%$
	$A$	$B$	$C$		$A$	$B$	$C$
$A$	0	6	4		0	62	39
$B$	5	0	3		53	0	27
$C$	3	2	0		28	22	0

Table 2. Table 2: Parameters used in the model to simulate the spread of HCV .

Parameter	Notation	Value	Source
Risk of transmission per contact
HCV, %	$\pi$	1	[6]
HIV, %	${\pi }^{\ast }$	0.63	[19]
HCV and HIV, %	${\pi }^{†}$	$6.3\times {10}^{-5}$	–
Semi-annual HCV spontaneous clearance
HCV mono-infected, %	$p$	26	[23]
HCV/HIV co-infected, %	$p^{\ast }$	10	[24]
Annual risk of death
HCV/HIV uninfected, ‰	$\mu$	10.4	[22]
HCV mono-infected, ‰	$\delta$	22.7	[22]
HIV mono-infected, ‰	${\mu }^{\ast }$	44.3	[22]
HCV/HIV co-infected, ‰	${\delta }^{\ast }$	54.9	[22]
Network characteristics
Mean degree	$⟨k⟩$	3.0	[9]
Scale parameter	$\alpha$	1.8	[9]
Exponential cutoff	$\kappa$	40	–
Anti-HCV treatment characteristics
Recovery rate, %	$\nu$	93	[20]
Proportion of treated patients, %	$q$	5.2	[21]
Other parameters
Per-partner sharing frequency, per month	$\omega$	5.33	[15]
Duration of the acute phase, years	$1/\sigma$	1/2	[23]
Unit time, years	–	1	–

Urban area	Population		Urban area	Population
Paris	12,341,418		Saint-Brieuc	171,721
Lyon	2,214,068		Béziers	165,498
Marseille - Aix-en-Provence	1,727,070		Montbéliard	162,326
Toulouse	1,270,760		Niort	152,721
Lille†	1,166,452		Vannes	150,860
Bordeaux	1,158,431		Chartres	145,735
Nice	1,004,914		Bourges	139,968
Nantes	897,713		Thionville	135,627
Strasbourg†	768,868		Châlon-sur-Saône	133,557
Rennes	690,467		Boulogne-sur-Mer	133,062
Grenoble	679,863		Maubeuge	129,931
Rouen	658,285		Arras	128,784
Toulon	611,237		Colmar	127,625
Montpellier	569,956		Blois	127,053
Douai - Lens	540,981		Calais	126,266
Avignon	515,536		Quimper	125,487
Saint-Etienne	512,830		Beauvais	125,095
Tours	483,744		Bourg-en-Bresse	122,806
Clermont-Ferrand	469,922		Laval	121,399
Nancy	434,479		La Roche-sur-Yon	117,965
Orléans	423,123		Creil	117,654
Angers	403,633		Cherbourg-Octeville	116,517
Caen	403,765		Tarbes	115,557
Metz	389,700		Belfort	114,077
Dijon	377,590		Alès	113,769
Béthune	368,633		Vienne	112,334
Valenciennes†	367,094		Agen	111,663
Le Mans	344,893		Saint-Quentin	111,474
Reims	317,611		Evreux	111,449
Brest	314,844		Roanne	107,209
Perpignan	309,962		Charleville-Mézières	106,835
Amiens	293,671		Montauban	105,654
Genève - Annemasse†	292,180		Cholet	104,917
Le Havre	290,890		Périgueux	102,417
Bayonne†	288,359		Sarrebruck - Forbach†	101,806
Mulhouse	284,739		Nevers	101,586
Limoges	282,971		Brive-la-Gaillarde	101,435
Nîmes	259,348		Ajaccio	100,643
Dunkerque	257,773		Mâcon	99,873
Poitiers	255,831		Carcassonne	97,801
Besançon	246,841		Albi	97,667
Pau	240,857		Compiègne	97,502
Annecy	221,111		Bastia	93,971
Chambéry	217,356		Epinal	93,891
Lorient	215,591		Fréjus	93,562
Saint-Nazaire	213,083		Bâle - Saint-Louis†	93,018
La Rochelle	207,211		Châteauroux	92,723
Troyes	191,505		Auxerre	92,307
Angoulême	180,593		Sète	91,101
Valence	175,636		Cluses	90,872

Equations

$$\displaystyle\eqalign{\frac{\mathrm{d}s(t)}{\mathrm{d}t}=\Lambda-\beta s(t)\big{(}a(t)+c(t)\big{)}+p\sigma a(t)+\gamma c(t)-\mu s(t),\\ \frac{\mathrm{d}a(t)}{\mathrm{d}t}=\beta s(t)\big{(}a(t)+c(t)\big{)}-(\sigma+\mu)a(t),\\ \frac{\mathrm{d}c(t)}{\mathrm{d}t}=(1-p)\sigma a(t)-(\gamma+\delta+\mu)c(t).}$$

$$s(t)=\sum _{k=1}^{N-1}{s}_k(t),a(t)=\sum _{k=1}^{N-1}{a}_k(t),c(t)=\sum _{k=1}^{N-1}{c}_k(t).$$

$${\theta }_k(t)=\sum _{k^{\prime }=1}^{N-1}{\rho }_N(k^{\prime }|k)\left(a_{k^{\prime }}(t)+c_{k^{\prime }}(t)\right)$$

$$\frac{\mathrm{d}{s}_k(t)}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{a}_k(t)}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{c}_k(t)}{\mathrm{d}t}$$

$$\rho (k)=\frac{k^{-\alpha }{e}^{-k/\kappa }}{{\mathrm{Li}}_{\alpha }\left(e^{-1/\kappa }\right)}$$

$$\displaystyle\theta_{k}(t)=\displaystyle\sum_{k^{\prime}=1}^{N-1}\rho_{N}(k^{\prime}|k)\Big{(}a_{k^{\prime}}(t)+c_{k^{\prime}}(t)+a^{*}_{k^{\prime}}(t)+c^{*}_{k^{\prime}}(t)\Big{)},$$

$$\displaystyle\theta^{*}_{k}(t)=\displaystyle\sum_{k^{\prime}=1}^{N-1}\rho_{N}(k^{\prime}|k)\Big{(}s^{*}_{k^{\prime}}(t)+a^{*}_{k^{\prime}}(t)+c^{*}_{k^{\prime}}(t)\Big{)},$$

$$\displaystyle\theta^{**}_{k}(t)=\displaystyle\sum_{k^{\prime}=1}^{N-1}\rho_{N}(k^{\prime}|k)\Big{(}a^{*}_{k^{\prime}}(t)+c^{*}_{k^{\prime}}(t)\Big{)}.$$

$$\frac{\mathrm{d}{s}_k(t)}{\mathrm{d}t}=\Lambda -k\beta s_k(t){\theta }_k(t)-k{\beta }^{\ast }{s}_k(t){\theta }_k^{\ast }(t)-k{\beta }^{†}{s}_k(t){\theta }_k^{\ast \ast }(t)$$

$$+p\sigma a_k(t)+\gamma c_k(t)-\mu s_k(t),$$

$$\frac{\mathrm{d}{s}_k^{\ast }(t)}{\mathrm{d}t}=-k\beta s_k^{\ast }(t){\theta }_k(t)+k{\beta }^{\ast }{s}_k(t){\theta }_k^{\ast }(t)+p^{\ast }\sigma a_k^{\ast }(t)+{\gamma }^{\ast }{c}_k^{\ast }(t)-{\mu }^{\ast }{s}_k^{\ast }(t),$$

$$\frac{\mathrm{d}{a}_k(t)}{\mathrm{d}t}=k\beta s_k(t){\theta }_k(t)-(\sigma +\mu )a_k(t)-k{\beta }^{\ast }{a}_k(t){\theta }_k^{\ast }(t),$$

$$\frac{\mathrm{d}{a}_k^{\ast }(t)}{\mathrm{d}t}=k\beta s_k^{\ast }(t){\theta }_k(t)+k{\beta }^{†}{s}_k(t){\theta }_k^{\ast \ast }(t)+k{\beta }^{\ast }{a}_k(t){\theta }_k^{\ast }(t)-(\sigma +{\mu }^{\ast })a_k^{\ast }(t),$$

$$\frac{\mathrm{d}{c}_k(t)}{\mathrm{d}t}=(1-p)\sigma a_k(t)-k{\beta }^{\ast }{c}_k(t){\theta }_k^{\ast }(t)-(\gamma +\delta +\mu )c_k(t),$$

$$\frac{\mathrm{d}{c}_k^{\ast }(t)}{\mathrm{d}t}=(1-p^{\ast })\sigma a_k^{\ast }(t)+k{\beta }^{\ast }{c}_k(t){\theta }_k^{\ast }(t)-({\gamma }^{\ast }+{\delta }^{\ast }+{\mu }^{\ast })c_k^{\ast }(t).$$

$$\frac{\mathrm{d}{s}_{k,A}(t)}{\mathrm{d}t}=-k\beta s_{k,A}(t){\theta }_{k,A}(t)-k{\beta }^{\ast }{s}_{k,A}(t){\theta }_{k,A}^{\ast }(t)-k{\beta }^{†}{s}_{k,A}(t){\theta }_{k,A}^{\ast \ast }(t)-\mu s_{k,A}(t)$$

$$\displaystyle-s_{k,A}(t)\left(\frac{\tau_{B,A}}{n_{A}+\tau_{B,A}}\right)\cdot\bigg{[}\beta\Big{(}a_{B}(t)+c_{B}(t)+a^{*}_{B}(t)+c^{*}_{B}(t)\Big{)}\bigg{.}$$

$$\displaystyle+\bigg{.}\beta^{*}\Big{(}s^{*}_{B}(t)+a^{*}_{B}(t)+c^{*}_{B}(t)\Big{)}+\beta^{\dagger}\Big{(}a^{*}_{B}(t)+c^{*}_{B}(t)\Big{)}\bigg{]}+p\sigma a_{k,A}(t)+\gamma c_{k,A}(t)$$

$$n_A=s_A(t)+a_A(t)+c_A(t)+s_A^{\ast }(t)+a_A^{\ast }(t)+c_A^{\ast }(t).$$

$$\frac{\mathrm{d}{s}_{k,A}(t)}{\mathrm{d}t}=-k\beta s_{k,A}(t){\theta }_{k,A}(t)-k{\beta }^{\ast }{s}_{k,A}(t){\theta }_{k,A}^{\ast }(t)-k{\beta }^{†}{s}_{k,A}(t){\theta }_{k,A}^{\ast \ast }(t)-\mu s_{k,A}(t)$$

$$\displaystyle-\displaystyle s_{k,A}(t)\sum_{\Phi\in\{B,C\}}\left(\frac{\tau_{\Phi,A}}{n_{A}+\sum_{\Phi\in\{B,C\}}\tau_{\Phi,A}}\right)\cdot\bigg{[}\beta\Big{(}a_{\Phi}(t)+c_{\Phi}(t)+a^{*}_{\Phi}(t)+c^{*}_{\Phi}(t)\Big{)}\bigg{.}$$

$$\displaystyle+\bigg{.}\beta^{*}\Big{(}s^{*}_{\Phi}(t)+a^{*}_{\Phi}(t)+c^{*}_{\Phi}(t)\Big{)}+\beta^{\dagger}\Big{(}a^{*}_{\Phi}(t)+c^{*}_{\Phi}(t)\Big{)}\bigg{]}+p\sigma a_{k,A}(t)+\gamma c_{k,A}(t).$$

$${\tau }_{A,B}=\frac{m_p{n}_A{n}_B}{n_B+n_C}$$

$$\frac{\mathrm{d}{s}_{k,\cdot }(t)}{\mathrm{d}t}=\Lambda -k\beta s_{k,\cdot }(t){\theta }_{k,\cdot }(t)-k{\beta }^{\ast }{s}_{k,\cdot }(t){\theta }_{k,\cdot }^{\ast }(t)-k{\beta }^{†}{s}_{k,\cdot }(t){\theta }_{k,\cdot }^{\ast \ast }(t)-\mu s_{k,\cdot }(t)$$

$$\displaystyle-\displaystyle s_{k,\cdot}(t)\sum_{\Phi}\left(\frac{\tau_{\Phi,\cdot}}{n_{\cdot}+\sum_{\Phi}\tau_{\Phi,\cdot}}\right)\cdot\bigg{[}\beta\Big{(}a_{\Phi}(t)+c_{\Phi}(t)+a^{*}_{\Phi}(t)+c^{*}_{\Phi}(t)\Big{)}\bigg{.}$$

$$\displaystyle+\bigg{.}\beta^{*}\Big{(}s^{*}_{\Phi}(t)+a^{*}_{\Phi}(t)+c^{*}_{\Phi}(t)\Big{)}+\beta^{\dagger}\Big{(}a^{*}_{\Phi}(t)+c^{*}_{\Phi}(t)\Big{)}\bigg{]}+p\sigma a_{k,\cdot}(t)+\gamma c_{k,\cdot}(t),$$

$$\frac{\mathrm{d}{s}_{k,\cdot }^{\ast }(t)}{\mathrm{d}t}=-k\beta s_{k,\cdot }^{\ast }(t){\theta }_{k,\cdot }(t)+k{\beta }^{\ast }{s}_{k,\cdot }(t){\theta }_{k,\cdot }^{\ast }(t)+p^{\ast }\sigma a_{k,\cdot }^{\ast }(t)+{\gamma }^{\ast }{c}_{k,\cdot }^{\ast }(t)-{\mu }^{\ast }{s}_{k,\cdot }^{\ast }(t)$$

$$\displaystyle+\displaystyle\sum_{\Phi}\left(\frac{\tau_{\Phi,\cdot}}{n_{\cdot}+\sum_{\Phi}\tau_{\Phi,\cdot}}\right)\cdot\bigg{[}\beta^{*}s_{k,\cdot}(t)\Big{(}s^{*}_{\Phi}(t)+a^{*}_{\Phi}(t)+c^{*}_{\Phi}(t)\Big{)}\bigg{.}$$

$$\displaystyle\bigg{.}-\beta s^{*}_{k,\cdot}(t)\Big{(}a_{\Phi}(t)+c_{\Phi}(t)+a^{*}_{\Phi}(t)+c^{*}_{\Phi}(t)\Big{)}\bigg{]},$$

$$\frac{\mathrm{d}{a}_{k,\cdot }(t)}{\mathrm{d}t}=k\beta s_{k,\cdot }(t){\theta }_{k,\cdot }(t)-(\sigma +\mu )a_{k,\cdot }(t)-k{\beta }^{\ast }{a}_{k,\cdot }(t){\theta }_{k,\cdot }^{\ast }(t)$$

$$\displaystyle+\displaystyle\beta s_{k,\cdot}(t)\sum_{\Phi}\left(\frac{\tau_{\Phi,\cdot}}{n_{\cdot}+\sum_{\Phi}\tau_{\Phi,\cdot}}\right)\cdot\Big{(}a_{\Phi}(t)+c_{\Phi}(t)+a^{*}_{\Phi}(t)+c^{*}_{\Phi}(t)\Big{)},$$

$$\frac{\mathrm{d}{a}_{k,\cdot }^{\ast }(t)}{\mathrm{d}t}=k\beta s_{k,\cdot }^{\ast }(t){\theta }_{k,\cdot }(t)+k{\beta }^{†}{s}_{k,\cdot }(t){\theta }_{k,\cdot }^{\ast \ast }(t)+k{\beta }^{\ast }{a}_{k,\cdot }(t){\theta }_{k,\cdot }^{\ast }(t)-(\sigma +{\mu }^{\ast })a_{k,\cdot }^{\ast }(t)$$

$$\displaystyle+\displaystyle\sum_{\Phi}\left(\frac{\tau_{\Phi,\cdot}}{n_{\cdot}+\sum_{\Phi}\tau_{\Phi,\cdot}}\right)\cdot\bigg{[}\beta s^{*}_{k,\cdot}(t)\Big{(}a_{\Phi}(t)+c_{\Phi}(t)+a^{*}_{\Phi}(t)+c^{*}_{\Phi}(t)\Big{)}\bigg{.}$$

$$\displaystyle+\bigg{.}\beta^{*}a_{k,\cdot}(t)\Big{(}s^{*}_{\Phi}(t)+a^{*}_{\Phi}(t)+c^{*}_{\Phi}(t)\Big{)}+\beta^{\dagger}s_{k,\cdot}(t)\Big{(}a^{*}_{\Phi}(t)+c^{*}_{\Phi}(t)\Big{)}\bigg{]},$$

$$\frac{\mathrm{d}{c}_{k,\cdot }(t)}{\mathrm{d}t}=(1-p)\sigma a_{k,\cdot }(t)-k{\beta }^{\ast }{c}_{k,\cdot }(t){\theta }_{k,\cdot }^{\ast }(t)-(\gamma +\delta +\mu )c_{k,\cdot }(t)$$

$$\displaystyle-\displaystyle\beta^{*}(a_{k,\cdot}(t)+c_{k,\cdot}(t))\sum_{\Phi}\left(\frac{\tau_{\Phi,\cdot}}{n_{\cdot}+\sum_{\Phi}\tau_{\Phi,\cdot}}\right)\cdot\Big{(}s^{*}_{\Phi}(t)+a^{*}_{\Phi}(t)+c^{*}_{\Phi}(t)\Big{)},$$

$$\frac{\mathrm{d}{c}_{k,\cdot }^{\ast }(t)}{\mathrm{d}t}=(1-p^{\ast })\sigma a_{k,\cdot }^{\ast }(t)+k{\beta }^{\ast }{c}_{k,\cdot }(t){\theta }_{k,\cdot }^{\ast }(t)-({\gamma }^{\ast }+{\delta }^{\ast }+{\mu }^{\ast })c_{k,\cdot }^{\ast }(t)$$

$$\displaystyle+\displaystyle\beta^{*}c_{k,\cdot}(t)\sum_{\Phi}\left(\frac{\tau_{\Phi,\cdot}}{n_{\cdot}+\sum_{\Phi}\tau_{\Phi,\cdot}}\right)\cdot\Big{(}s^{*}_{\Phi}(t)+a^{*}_{\Phi}(t)+c^{*}_{\Phi}(t)\Big{)}.$$

Full text

Keywords: Networks, Epidemic modelling

1 Introduction

Modelling the propagation of epidemics is a long-standing endeavour in the medical and mathematical sciences. In the past decade or so, multiple developments in the science of complex systems in general, and networks in particular, have led to new models, along with the possibility to test their predictions on datasets, the size of which keeps increasing (see for instance [1, 2] and references therein).

Particularly challenging are the study and modeling of epidemics involving specific, sometimes marginal, segments of the general population. For instance, viruses and diseases that are particularly prevalent among people who inject drugs (PWID). Things become even more involved if one tries to take into account correlations between interrelated epidemics, such as Hepatitis C virus (HCV) and AIDS (HIV) among PWID populations [3].

In this paper, we introduce a simple model of propagation for HCV, with potential HIV co-infection, on a network of networks, representing PWID populations in distinct but interconnected urban areas. This model is intended as a first step towards a more comprehensive model, and to serve as a basis for comparison with public health statistical data when it becomes available.

In our model, the PWID population of each urban area is a network of individuals, and the fact that people may travel or even move houses from one city to another makes for a second level of links, between cities themselves. We choose to represent this second level via the introduction of special nodes at the first level, that will represent at the individual level in each city, contacts with individuals from other cities. This means that if cities $A$ and $B$ are linked (in the sense that a number of individuals travel from one of them to the other), then we add a node in the network representing individuals in city $A$. The degree and other characteristics of this special node in city $A$ will depend on the number of individuals from city $B$ travelling to city $A$ and on the size of the infected population in city $B$. This makes for a simple treatment of the links between cities. We also adapt the standard susceptible-infected-recovered (SIR) equations to reflect the specifics of the HCV epidemics with possible HIV co-infection, as well as to take into account the special type of network with which we are working. The construction of this simple model is detailed in Section 2 and 3.

In Section 4, we examine, via numerical simulations, how our model works, first on a toy-country with just three cities and then on the network of France’s 100 largest urban areas. We find that simulations run for different values of the parameters lead to sensible stylized facts in terms of the extent and speed of propagation of an HCV epidemics in France.

2 Epidemic models

2.1 Initial model of the spread of Hepatitis C

Let us consider the compartmental model describing the infection dynamics of Hepatitis C virus (HCV) displayed in Figure 1. This model is based on a susceptible-infected-susceptible (SIS) model which is a derivative of the classic susceptible-infected-recovered (SIR) model introduced in [4].

We denote by $S$ the number of susceptible individuals, by $A$ the number of individuals in the acute phase of the infection and by $C$ the number of individuals in the chronic phase of the infection. Assuming that the number of individuals in the total population is constant over time, which is not a strong assumption for large populations, we write $S(t)+A(t)+C(t)=N$. In particular, we assume that the number of people who died from HCV or from some other cause equals the number of people joining the population as new susceptibles.

From now on, we work in terms of proportions rather than absolute numbers. The proportion of susceptible, acute-HCV infected and chronic-HCV infected at time $t$ are denoted by $s(t),a(t)$ and $c(t)$ respectively. Writing $\beta$ for the transmission rate, $\gamma$ for the recovery rate, $\delta$ for the death rate, $\mu$ for the natural mortality rate, $\Lambda$ for the joining rate, $1/\sigma$ for the duration of the acute phase and $p$ for the spontaneous clearance, i.e. the probability that a given individual clears the infection, one can write the following set of differential equations:

[TABLE]

Since HCV is an asymptomatic disease, very limited individuals are diagnosed during the acute phase of the infection. Thus we assume that only individuals in the chronic phase have access to a treatment, and only people in the chronic phase can recover.

Note that the transmission rate $\beta$ may be written as $\beta=\omega\times\pi$, where $\omega$ is the contact rate (the number of contacts per unit time) and $\pi$ the transmission risk (the risk of infection for a given contact). Moreover, as not every patient is actually treated, $\gamma$ may be written as $\gamma=\nu\times q$ where $\nu$ is the recovery rate per patient and $q$ the proportion of treated patients.

The model defined in the set of equations (2.4), also called mass action model, offers a simple and convenient manner of describing the dynamics of an epidemic in a large population. However, it fails to incorporate population heterogeneities and individual interactions that may greatly affect the virus’s rate of transmission.

Regarding HCV, it has been stated that the primary route of transmission in the developed world is intravenous drug use (IDU) [5]. As people who inject drugs (PWID) tend to share drug injecting equipment with a limited number of partners (see for example the social network described in [6]), using a simple mass action model might by questionable.

Hence, assuming that individuals with the same number of sharing partners present the same risk in the infectious process, we opt for the heterogeneous mean field approach, first introduced in [7]. Each individual is assigned a degree (a number of partners) $k$, where $1\leqslant k\leqslant N-1$, and constitutes an element of the network (also referred to as a node, in graph theory). As individuals with [math] partner do not contribute to the spread of the epidemic, only individuals with at least one partner will be considered. Figure 2 offers a visual representation of a static network of sharing partners infected by HCV. For every individual dying of HCV, it is assumed that a new susceptible individual enters the network, substituting for them.

Let us denote by $\rho$ the degree distribution, i.e. $\rho(k)$ is the proportion of individuals with degree $k$. Writing $s_{k}(t)$ (resp. $a_{k}(t)$ and $c_{k}(t)$) for the proportion of susceptible (resp. acute-HCV infected and chronic-HCV infected) individuals with degree $k$ at time $t$, one has:

[TABLE]

Referring now to [7], the proportion of infected partners for an individual with degree $k$ is given by:

[TABLE]

where $\rho_{N}(k^{\prime}|k)$ is the proportion of partners with degree $k^{\prime}$ among all partners of an individual with degree $k$. Under the assumption that all nodes of a given degree are statistically equivalent, one finds that $\rho_{N}(k^{\prime}|k)=\left(k^{\prime}\rho(k^{\prime})\right)/\langle k\rangle$, where $\langle k\rangle$ is the mean degree of the network [8]. As $\rho_{N}(k^{\prime}|k)$ then does not depend on $k$, we shall simply write $\rho_{N}(k^{\prime})$. Thus (2.4) becomes:

[TABLE]

To determine which degree distribution should be used to model the drug equipment sharing pattern between PWID, we refer to [9] and [10]. The authors study the social network of PWID in the Bushwick neighborhood in Brooklyn, New York, and analyse its degree distribution. This analysis showed that the resulting degree distribution could be assimilated to a power law distribution. Typical of scale-free networks, power law distributions reflect a preferential attachment dynamics [11], which is consistent with the fact that new PWID network members tend to join a particular injector network via someone who is already an influential part of it.

Henceforth, we assume that the social network of PWID admits a scale-free topology. In particular, one can write $\rho(k)\sim k^{-\alpha}$ where $\alpha$ is called scale parameter. However, according to [12], it appears that power laws with exponential cutoff tend to fit natural networks better than pure power laws, as almost surely there are no individuals with an extremely large number of partners. Thus, based on this result and for an arbitrary cutoff $\kappa$, one can write:

[TABLE]

where $\mathrm{Li}_{\alpha}$ is the polylogarithm of order $\alpha$.

2.2 Incorporating contact with HIV-infected individuals

The model developed in the previous subsection does not take into account the fact that among PWID, people infected by HCV are often co-infected by HIV (human immunodeficiency virus), which greatly increases the risk of death [3].

Therefore, we introduce three new compartments representing individuals infected only by HIV ($S^{*}$), individuals in the acute phase of the HCV infection and infected by HIV ($A^{*}$) and individuals in the chronic phase of the HCV infection and infected by HIV ($C^{*}$). We now have $S(t)+A(t)+C(t)+S^{*}(t)+A^{*}(t)+C^{*}(t)=N$.

We also define $\gamma^{*}$ and $\delta^{*}$, the recovery and death rates of HCV for people already infected by HIV. As well as $\beta^{*}=\omega\times\pi^{*}$ the transmission rate of HIV, $\beta^{\dagger}=\omega\times\pi^{\dagger}$ the transmission rate of both HCV and HIV and $p^{*}$ the spontaneous clearance of HCV for people infected by HIV. Finally, we write $\mu^{*}$ for the mortality rate of PWID infected by HIV. A visual representation of the model is given in Figure 3. As previously, proportions of the different quantities are denoted with small letters.

We now introduce the probabilities that an individual with degree $k$ has an infected partner whichever the infection. Hence, we have the probability that an individual with degree $k$ has an HCV infected partner:

[TABLE]

the probability that an individual with degree $k$ has an HIV infected partner:

[TABLE]

and the probability that an individual with degree $k$ has an HCV/HIV co-infected partner:

[TABLE]

From those equations and based on the model of Figure 3, one easily rewrites (2.6):

[TABLE]

Such a framework allows for the study of an epidemic on a given network. However, when dealing with a large population, it might be wiser to deal with several subnetworks connected with each other. Indeed, when dealing with an epidemic, say at a country scale, computing a subnetwork for each city and allowing human movements between them might be more accurate than computing a single giant network. Based on this hypothesis, we propose to study, in the next section, the spread of the epidemic of HCV defined by (2.2) across several subnetworks.

3 Epidemic propagation on a metapopulation model

Metapopulation models, first described in [13], consist in modeling the interactions of several spatially separated populations. Those models have been extensively used in the field of epidemiology to describe the spread of infectious diseases at a large scale by dividing the population into subpopulations corresponding to different households, cities, etc. Depending on human movements, a subpopulation, or patch, containing infected individuals is allowed to interact with other patches spreading the infection to neighboring subpopulations. We refer the reader for instance to [1] and [2] for an illustration of the role of the global aviation network in the spread of the SARS infection.

Here we adapt this type of model to study the spread of HCV on a network of cities, assuming that the population of each city can be considered as a patch of a global metapopulation model. Each subpopulation admits a scale-free structure as defined in Section 2.1 and the spread of the virus is governed by the set of differential equations (2.2). To keep things simple while making them concrete, we work in this section on networks of two and three cities.

Two-city network. Let us consider two cities $A$ and $B$ with injecting drug populations $N_{A}$ and $N_{B}$ respectively. Assuming that these populations are organized in scale-free networks, we denote respectively by $\langle k_{A}\rangle$ and $\langle k_{B}\rangle$ their mean degree and by $\rho_{A}(k)$ and $\rho_{B}(k)$ their degree distributions. We update the notations in (2.2) by adding subscripts corresponding to cities, e.g. $s_{A}(t)$ and $s_{B}(t)$ denote the proportions of susceptible individuals in cities $A$ and $B$ respectively.

To model the effect of $B$ on the spread of the disease in $A$, we propose to add a node in the network of $A$ which is interacting with every node of $A$. In the same way, we add an extra node into the network of $B$ to model the effect of $A$ on $B$. These extra nodes allow the population of $A$ (resp. $B$) to connect with the new individuals moving from $B$ to $A$ (resp. from $A$ to $B$). A visual representation of the process is given in Figure 4.

The extra node in the network of $A$ contains the main characteristics of $B$, i.e. the proportion of infected individuals in $B$ and $\tau_{B,A}$ the number of people travelling from $B$ to $A$. Similarly, the extra node in the network of $B$ contains the main characteristics of $A$. Hence, by focusing on the first equation of (2.2), one can write:

[TABLE]

where $n_{A}$ is the population still alive in city A. In particular, $n_{A}$ is defined as:

[TABLE]

Three-city network. Let us introduce a third city $C$ with population $N_{C}$. As above, we assume that $C$ is organized in a scale-free network with mean degree $\langle k_{C}\rangle$ and degree distribution $\rho_{C}(k)$. To model the interactions between $A$ and $C$ (resp. $B$ and $C$), we add another node in the graph of $A$ (resp. $B$) containing the main characteristics of $C$. See Figure 5 for a graphical visualization of the process in the particular case of $A$.

Adding the new node corresponding to city $C$, equation (3) becomes:

[TABLE]

Full details of the complete set of differential equations are available in Appendix A.1.

4 Simulations

4.1 Parameter values

Let us first consider three abstract cities $A$, $B$ and $C$ which PWID communities of 10000, 8000 and 5000 injecting drug users respectively. Cities with larger populations tend to be more attractive, so we use a constrained gravitational model [14] to determine the proportion of a city’s travelling population driven to a particular city among all other cities:

[TABLE]

where $m_{p}$ is the moving population percentage. Values of the different population movements were rounded up and reported in Table 1. Those percentages were kept low to account for the fact that the population of a given city is far greater than the moving population.

As far as network parameters are concerned, very limited data is available for their estimation. Hence, we use the values obtained by [9], which were $\alpha=1.8$ for the scale parameter and $\langle k\rangle=3.0$ for the mean degree. Based on this, we initiate 1000 simulations of the degree distribution with $\alpha=1.8$ and different values of the exponential cut-off until reaching the value of 3.0 for the mean degree. The resulting cutoff was $\kappa=40$. Figure 6 illustrates the degree distribution of a network of N=10,000.

As illustrated on Figure 6, degrees may attain high values ($>10^{2}$). Since the model described in A.1 is degree based, a high degree can lead to a very large amount of differential equations. Hence, for the sake of practicality, we decide to split the population into four arbitrary groups: people with degree 1, people whose degree lies between 2 and 5, people whose degree lies between 6 and 10 and people whose degree is greater than 10.

The values of the transmission risks of HCV and HIV are difficult to estimate. Several studies reported different values for the per contact probability of HCV transmission among PWID, ranging from 0.5% to 10% (see [6, 15, 16, 17, 18]). As for the transmission of HCV, we refer to [6] who estimated a transmission risk of 1% in a network of PWID in Australia. Considering the transmission of HIV, we refer to [19] who reported 63 transmissions for 10,000 exposures for the US epidemic. These transmission risks represent the probabilities of getting infected per act of sharing needles. Assuming that these probabilities are independent, the transmission risk $\pi^{\dagger}$ of both HCV and HIV is given by $\pi^{\dagger}=\pi\times\pi^{*}$.

As for the frequency of needle/syringe sharing, difficulties arise when one tries to obtain accurate estimations. Referring to [15] we choose the value of 16 sharing acts per month. However, as this frequency is shared among the different partners of a given individual, we assume a per-partner sharing frequency of $16/3=5.33$, where 3 is the mean number of partners.

The values of the transmission risks introduced above reflect somehow the true HCV and HIV transmission risks for PWID. However, for illustrative purposes, we propose to temporary choose greater values for those risks by multiplying the transmission risks $\pi$, $\pi^{*}$ and $\pi^{\dagger}$ by some constant $\lambda$ to observe the impact of migrating people on the spread of the epidemic.

Considering the treatment of HCV, we decided to focus on the recent direct-acting antivirals (DAAs) which leads to high recovery rates. Indeed, based on a recent ANRS press release, the recovery rate of the new DAAs has been estimated up to $93\%$ for HCV/HIV co-infected patients (see [20]). Similar results were observed for HCV mono-infected patients. Due to the fact that HCV is an asymptomatic disease, a large part of HCV infected patients are not aware of their serostatus. Globally, this involves a very low proportion of infected patients accessing to health care. Indeed, although France has one of the highest treatment rates in Europe, this proportion hardly reaches 5.2% (see [21]).

As for the mortality rates of PWID infected by HCV, HIV or HCV and HIV, we refer to the work of [22] who calculated all-cause and cause-specific crude mortality rates (per 1000 person-years) and standardized mortality ratios for PWID in the Netherlands. A summary of the parameters used in the model are given in Table 2.

Since we want to quantify the impact of the interactions between cities, we assume that the epidemic starts in one of the three cities and observe the impact on the two others. We initiate the epidemic with a relatively low number of infected individuals in city $A$, that is $s_{A}(0)=0.985$, $c_{A}(0)=0.01$ (HCV mono-infected individuals in chronic phase) and $c^{*}_{A}(0)=0.005$ (HIV/HCV co-infected individuals in chronic phase). The values chosen for the initiation of the model are purely arbitrary.

Due to the random nature of the network framework, we will focus on the average of each quantity over 1000 simulation runs.

4.2 Results on an abstract three-city network

Figures 7, 8 and 9 display the evolution of the susceptible, the HCV mono-infected and HIV/HCV co-infected populations in the whole system (cities $A$, $B$ and $C$). The different curves of the three figures represent the population dynamics when no migration is assumed (blue markers), 0.1% of the population is migrating (red markers) and 1% of the population is migrating (green markers). Each quantity displayed in Figures 7, 8 and 9 is the sum of the corresponding quantity over the four degree groups.

In Figure 7, one sees that the susceptible population gradually depletes with the percentage of migrating population. As more and more people move from city $A$, the epidemic spreads faster and has a stronger impact on cities $B$ and $C$. As expected, the spread of the epidemic is even stronger when coefficient $\alpha$ increases i.e. when the scale parameter of the degree distribution increases and individuals with larger degrees are present. One can see that the blue curves reach a minimum value between $60$% and $70$%. This phenomenon is due to the fact that the susceptible populations in cities $B$ and $C$ stay intact (no initial infected individual nor migration from infected cities).

In Figure 8, one observes the same dynamics as in Figure 7. The proportion of HCV mono-infected individuals in the population is greatly impacted by the percentage of migrating people and by the value of $\alpha$. Similarly, the blue curves here correspond to the proportion of HCV mono-infected individuals in the city $A$. This explains why that proportion is never higher than 33.3%.

Lastly, Figure 9 exhibits the same evolutions as the two previous figures. One can particularly observe on Figure 9 that the evolution of the HIV/HCV co-infected population is greater than the evolution of the HCV mono-infected population of Figure 8. This is due to the fact that, contrary to HCV, no treatment for HIV is considered in the model. Thus, each HCV mono-infected individual tends to get also infected by HIV as time goes by.

In the next section we study the spread of the infection on the network of France’s $100$ largest urban centers. In particular, we examine the influence of the starting point of the epidemic.

4.3 Results on the network of France’s 100 largest urban areas

We generalize the metapopulation model of Section 4.1 by studying the spread of HCV across a whole country: France. As we want our predictions to be relatively accurate without considering every single city, we choose to focus on the first hundred biggest urban areas in France (metropolitan France only) for a population of 43,187,838 individuals in 2012. Indeed, in [25] three main drivers of mobility for PWID were identified: legal problems, entering drug treatment programs, and drug tourism. We assume that the great majority of PWID tend to move towards and between the main urban areas.

Then, the population of these urban areas is multiplied by the prevalence of PWID in France estimated at $5.9$ ‰  by the OFDT (Observatoire français des drogues et des toxicomanies, the French monitoring centre for drugs & drug addiction) in 2006 [26]. Since this prevalence has been estimated for the population aged $15$ to $64$ ($63,8$% of the total French population in 2012, i.e 40,433,870 individuals), we multiply the population of the French urban areas by $40433870/43187838=0.94$ for both populations to match. The original populations retained for the simulations are reported in Appendix A.2.

Regarding the origin-destination matrix, no data is available for the annual migration of PWID. Hence, we refer to a French database of residential mobility during a five-year period published on the website of INSEE (Institut national de la statistique et des études économiques, the French office for national statistics). Values are divided by 5 and rounded to obtain an annual origin-destination matrix. As we focus on PWID, values are multiplied by 5.9 ‰.

To study the impact of migrating populations and of the size of the initial city infected on the spread of the infection, we shall consider two different starting points. We first choose Paris as a starting point since this city has the highest level of interactions with nearly all other French cities. To select the second starting point, we need a city which possesses a far lower level of interactions with all cities other than Paris. However, to prevent the snowball effect, it is better to also focus on a city which has very limited interactions with Paris. Hence, we choose Forbach, a small urban area in the North-Eastern part of France, as the starting point in a second set of simulations.

Referring to Section 4.1, we set the proportion of susceptible individuals to $98.5$% in the city where the epidemic starts. The proportions of HCV mono-infected (chronic phase) and HIV mono-infected patients are set to $1$% and $0.5$% respectively in that same city. All other cities are initiated at $100$% of susceptible individuals. Similarly, we assume that each urban area is organized according to a scale-free network, the degree distribution of which follows a power law with exponential cutoff and a multiplicative constant $\lambda$ of $10$. Parameters considered for this model are reported in Table 2.

Results. Evolution of the different mean quantities are represented in Figure 10 where the two different scenarios are exhibited in plain and dashed curves. One notices that, as expected, the epidemic spreads much more easily to the whole country when starting from Paris rather than Forbach. Indeed, as Paris has many more interactions with other urban areas than Forbach, the epidemic spreads faster. Particularly, one sees that the proportion of susceptible individuals depletes rapidly between the third and the eighth year when the epidemic starts in Paris while this proportion depletes between the seventh and the fourteenth year when the epidemic starts in Forbach. Similarly, one oberves that the proportion of HCV mono-infected patients increases over a few years and then decreases due to the superinfection with HIV.

Figures 11 and 12 give a comparison of the spread of the infection in the country. As intuitively expected, the epidemic spreads further and faster when starting in Paris. Indeed, ten years after the beginning of the epidemic, the high majority of urban areas admits a proportion of HCV infected individuals (mono-infection or co-infection) superior to $40$% when the epidemic outbreak takes place in Paris. For comparison, when the epidemic breaks out in Forbach, not all urban areas are affected by the infection even after ten years. Finally, one notes that the last epidemic-affected city is Cluses in both cases. This can be explained by the fact that this urban area presents the lowest level of interaction with other French cities (it only interacts significantly with Genève-Annemasse).

5 Conclusion

In this paper, we have introduced as a metapopulation model a compartmental model for the propagation of HCV on a network of networks, representing PWID populations in distinct but interconnected French cities.

Tested on a three-city network, it exhibits, as expected, an epidemic spread correlated with the amount of interactions between the cities: the epidemic grows faster and stronger when the amount of interactions between cities is high. Applied on a network of $100$ interconnected urban areas (corresponding to the $100$ lqrgest urban areas in France), the model yields very sensible stylized facts concerning the epidemic spread, in particular when chosing as starting points two different cities differing from each other by the amount of interactions they have with other urban areas. In particular, sensible differences are observed in terms of prevalence and propagation.

Our model combines two important features. First, it takes into account heterogeneities in the PWID population, which plays a crucial role in the transmission process. Secondly, it uses a metapopulation model to describe the infectious dynamics at a country scale.

Indeed, the use of a metapopulation model allows to divide the general PWID population into several PWID subpopulations and thus to better handle population heterogeneities at the subpopulation level with a contact network. Also, the model is relatively easy to handle, thanks to the use of extra-nodes working, in each network coorsedpong to a city, as proxys for individuals coming from other cities.

Having obtained sensible outcomes with this model, our next steps will be to obtain reliable data in order to calibrate and test it further. This is particularly challenging when dealing with PWID population. Indeed, whether it be needle sharing behavior for the estimation of HCV transmission or the degree of each individual for the establishment of a reasonable contact network, the lack of data forced us to make strong assumptions. Moreover, no cessation of treatment was taken into account despite the fact that the rate of non-compliance among PWID may be significant. Regarding the population moving from city to city, our only source was a French database of residential mobility which may misestimate the actual migrating dynamics of the PWID population.

It will also be of particular interest to explore part of the phase space arising when the values chosen for certain crucial parameters (such as the transmission rate) are allowed to vary across ranges away from their current reported or assumed typical value. This may reveal some unexpected transitions and critical values in the phase space corresponding to some of the parameters.

This study was funded by a CIFRE grant to A.N. from the Bristol-Myers Squibb company, subsidized by the ANRT (CIFRE convention number: 2012/0879).