Repositorio Institucional de Documentos

Modelos metapoblacionales para la difusi´on de epidemias Trabajo de fin de grado Alberto Aleta Casas Director Yamir Moreno Vega Departamento de F´ısica Te´orica Universidad de Zaragoza Zaragoza, Junio 2014 According to the evidence put forward in the preceding pages the space-time events in the body of a living being which correspond to the activity of its mind, to its selfconscious or any other actions, are (considering also their complex structure and the accepted statistical explanation of physicochemistry) if not strictly deterministic at any rate statistico-deterministic. What is life? Erwin Schr¨odinger, 1944 ´ Indice 4 ´ Indice 1 Introducci´on 6 2 Fundamentos de la teor´ıa de redes 7 2.1 Lamatrizdeadyacencia...................................... 7 2.2 Grado de un v´ertice y distribuci´on de grado . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Distribuci´on power-law ....................................... 8 2.4 Redes scale-free ........................................... 9 3 Difusi´on de epidemias 13 3.1 Modeloscompartimentales..................................... 13 3.2 Homogeneous mixing ........................................ 14 3.3 Difusi´on en redes scale-free .................................... 17 4 Modelos metapoblacionales 21 4.1 La difusi´on en metapoblaciones: procesos reacci´on-difusi´on . . . . . . . . . . . . . . . . . . 21 4.2 Global invasion threshold ..................................... 23 4.3 Aplicaciones: la pandemia de gripe A/H1N1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 5 Estudio de la persistencia en metapoblaciones 29 5.1 Modelo ............................................... 29 5.2 Resultados.............................................. 30 6 Conclusiones 33 Introducci´on 6 1 Introducci´on Un sistema complejo es aquel que muestra comportamientos que son consecuencia de la interacci´on entre los elementos que lo constituyen y que no pueden explicarse del an´alisis de dichos elementos por separado. A estos comportamientos se les denomina comportamientos emergentes o colectivos[1]. Encontramos sistemas complejos en multitud de disciplinas como la f´ısica, la biolog´ıa o la econom´ıa. Sin embargo, detr´as de cada uno de ellos siempre existe una red que define las interacciones entre sus componentes, la cual es fundamental para comprender el funcionamiento de estos sistemas. Un proceso que resulta de particular inter´es en el estudio de los sistemas complejos es el fen´omeno de la difusi´on en redes complejas, capaz de explicar desde la propagaci´on de una enfermedad hasta la de un virus de ordenador o la de un rumor. La introducci´on de la f´ısica estad´ıstica, la teor´ıa de transiciones de fase y la de fen´omenos cr´ıticos ha resultado de gran utilidad para describir el comportamiento macrosc´opico de los brotes epid´emicos. En particular, la aproximaci´on de campo medio ha resultado un gran ´exito al conseguir reducir enormemente los grados de libertad del sistema[2]. En los ´ultimos a˜nos el inter´es por el estudio de las epidemias se ha incrementado enormemente. Vivimos en una sociedad cada d´ıa m´as grande y conectada, lo que es el caldo de cultivo perfecto para una epidemia. As´ı, mientras que la peste negra tard´o 4 a˜nos en pasar de Francia a Suecia en el siglo XIV, en el a˜no 2009 la gripe A se propag´o de M´exico a Espa˜na en apenas un mes. Los modelos cl´asicos de difusi´on de epidemias se centran en estudiar la propagaci´on dentro de poblaciones aisladas. Sin embargo, hoy en d´ıa se hace necesario extender estos modelos para abarcar m´ultiples poblaciones separadas espacialmente que no son independientes entre s´ı. Es en este contexto en el que surgen los modelos metapoblacionales. El presente trabajo tiene como objetivo ´ultimo desarrollar un modelo computacional que nos permita estudiar la persistencia de una epidemia en una metapoblaci´on. Para conseguirlo ser´a necesario analizar y comprender c´omo se puede modelar la difusi´on de una epidemia en una poblaci´on y el papel que juega la estructura de la red subyacente en el proceso. Por ello, el texto se organiza como sigue: i) Comenzaremos con una breve introducci´on a la teor´ıa de redes en la que definiremos el concepto de grado y su influencia en la estructura de una red. ii) A continuaci´on estableceremos los principios b´asicos que gobiernan la difusi´on de las epidemias. Formularemos las hipot´esis de “homogeneous mixing” y de campo medio que formar´an una parte fundamental de los modelos metapoblacionales. iii) Integraremos todos estos conceptos en la construcci´on de un modelo metapoblacional y determinaremos algunas de sus propiedades m´as importantes. iv) Finalizaremos aplicando todos los conocimientos adquiridos en el desarrollo de un modelo computacional. Estudiaremos la persistencia de una epidemia en una metapoblaci´on en funci´on del tiempo promedio de inmunizaci´on de los individuos, de su movilidad y de la estructura de la red. Los resultados que obtengamos podr´an servir como punto de partida para futuros desarrollos te´oricos. Fundamentos de la teor´ıa de redes 7 2 Fundamentos de la teor´ıa de redes Un grafo es un conjunto de v´ertices unidos por enlaces (edges). Esta denominaci´on se emplea al hablar de la representaci´on matem´atica, mientras que para sistemas reales se prefiere la palabra red[3]. As´ı, en f´ısica, tenemos redes formadas por sites ybonds; en sociolog´ıa redes de actores yties; y en los sistemas complejos redes de nodos unidos por links. A pesar de la aparente simplicidad de este esquema est´a claro que existen numerosas diferencias entre una red cristalogr´afica, una red de transporte o una red metab´olica. A continuaci´on veremos algunas de las propiedades m´as importantes de las redes, en especial aquellas relacionadas con el tipo de redes que vamos a emplear para estudiar la difusi´on de epidemias. 2.1 La matriz de adyacencia Existen muchas formas de representar una red matem´aticamente, pero la m´as c´omoda para realizar c´alculos es emplear una matriz de adyacencia. Figura 1: Red de ´areas y departamentos de f´ısica en la Universidad de Zaragoza (distribuci´on espacial). Se define como matriz de adyacencia a la matriz A cuyos elementos Aij valen 1 si existe un enlace entre los v´ertices iyjy 0 en cualquier otro caso. Por ejemplo, la matriz de adyacencia correspondiente a la red mostrada en la figura 1 ser´a A=         011000 100100 100100 011011 000101 000110         En ocasiones es ´util asociar a los enlaces un n´umero real en lugar de un 0 o un 1. A estas redes se les denomina redes pesadas. Un ejemplo lo podemos encontrar en el modelo tigh-binding de la f´ısica del estado s´olido, donde los enlaces entre sites se representan por un par´ametro tcuyo valor depende de la integral de solapamiento entre los ´orbitales de sites vecinos[4, p.143]. 2.2 Grado de un v´ertice y distribuci´on de grado El grado de un v´ertice es igual al n´umero de enlaces que posee. As´ı, en una red de Nv´ertices el grado kide un nodo ise puede expresar en t´erminos de la matriz de adyacencia como ki= N X j=1 Aij Sea Lel n´umero total de enlaces. Como cada enlace est´a presente en dos v´ertices tendremos que 2L= N X i=1 ki⇒L=1 2X ij Aij Con estas relaciones el grado medio hkide un v´ertice se puede expresar como hki=1 N N X i=1 ki=2L N Fundamentos de la teor´ıa de redes 8 Un tipo especial de redes son aquellas cuyos v´ertices tienen todos el mismo grado. A estas redes se les conoce como redes regulares y son las que nos solemos encontrar al estudiar el modelo tight-binding mencionado anteriormente (periodic lattice). Se conoce como distribuci´on de grado a la distribuci´on de probabilidad de un grado en la red. Definimos como P(k) a la fracci´on de v´ertices de una red que poseen grado k. As´ı, P(k) representa la probabilidad de que un v´ertice aleatorio de la red sea de grado k. Aunque la distribuci´on de grado no determina completamente una red, s´ı que es una caracter´ıstica fundamental de su estructura y nos permite conocer muchas de sus propiedades. 2.3 Distribuci´on power-law El modelo m´as sencillo en el que podemos pensar para generar una red es el de Erd˝os-Renyi (ER)[5]. En este modelo cada par de nodos se une con la misma probabilidad p, de forma que se obtiene una red aleatoria. La distribuci´on de grado de una red aleatoria finita sigue una distribuci´on binomial P(k) = N−1 kpk(1 −p)(N−1)−k De forma que cuando N→ ∞ la distribuci´on tiende a una distribuci´on de Poisson P(k) = e−hkihkik k! Veamos qu´e implica este resultado. El n´umero medio de amigos que un usuario estadounidense tiene en Facebook es de 350[6] siendo el n´umero total de usuarios de ese pa´ıs ∼2·108. As´ı, si las relaciones sociales siguiesen el patr´on de una red aleatoria, la probabilidad de que un usuario tuviese 100 amigos ser´ıa de ∼5·10−48 y la de tuviese 600 amigos ∼4·10−26. En otras palabras, en una red aleatoria la mayor´ıa de los v´ertices tienen aproximadamente el mismo grado. Esto contradice claramente nuestra experiencia, pues es comprensible suponer que una persona famosa tendr´a bastantes m´as amigos que la media. Ante este resultado surge una pregunta, ¿es este modelo relevante para sistemas reales? La respuesta es no[7]. A pesar de ello se mantuvo vigente hasta 1999 debido a la carencia de los datos necesarios para comprobar experimentalmente su validez. No obstante, a d´ıa de hoy se sigue empleando a nivel acad´emico debido a su sencillez. El punto de inflexi´on se produjo en 1999. En este a˜no Barab´asi et al. realizaron un mapa de Internet donde los nodos eran las p´aginas web y los links los enlaces entre unas y otras. El resultado fue que el 80% de los nodos ten´ıan 4 links o menos, mientras que un 0.01% de ellos enlazaban con m´as de 1000 p´aginas, todo lo contrario a lo que uno encuentra en una red aleatoria[8]. Tras repetir el an´alisis con otras redes llegaron a la conclusi´on de que la probabilidad de que un v´ertice de una red real interaccione con otros k v´ertices (i.e. que tenga grado k) decae siguiendo una ley de potencias o power-law[9], P(k)∼k−α A pesar de su simplicidad, este tipo de distribuci´on tiene algunas propiedades muy interesantes. Comencemos estudiando el valor de la constante de normalizaci´on. Sabemos que toda distribuci´on de probabilidad tiene que cumplir que ∞ X k=0 P(k) = 1 Difusi´on de epidemias 15 n´umeros no est´an determinados un´ıvocamente. Si la enfermedad volviese a extenderse entre la poblaci´on, aunque parti´esemos de las mismas condiciones iniciales, los resultados ser´ıan ligeramente diferentes. Sin embargo, por ahora supondremos que estamos ante un proceso deterministac. Supongamos que en cada paso de tiempo cada infectado entra en contacto con una persona. La probabilidad de que ´esta pertenezca al estado susceptible ser´a S(t)/N donde Nes el n´umero total de personas. As´ı, por unidad de tiempo, la probabilidad de que un infectado contagie a uno sano ser´a βS(t)/N. Como el n´umero total de enfermos es I(t), en cada paso de tiempo se estar´an generando βI(t)S(t)/N nuevas infecciones. De forma similar, µI(t) personas pasar´an al estado inmune y γR(t) volver´an a ser susceptibles. Si escribimos esto en forma de ecuaci´on, S(t+δt) = S(t)−βS(t)I(t) Nδt +γR(t)δt ⇒dS(t) dt =S(t+δt)−S(t) δt =−βS(t)I(t) N+γR(t) Procediendo de forma an´aloga con los otros dos estados y normalizando (s≡S/N,i≡I/N yr≡R/N) se obtiene el sistema de ecuaciones diferenciales que determina la evoluci´on de una epidemia en el modelo SIRS bajo la hip´otesis de homogeneous mixing, ds(t) dt =−βs(t)i(t) + γr(t) (3.1) di(t) dt =βs(t)i(t)−µi(t) (3.2) dr(t) dt =µi(t)−γr(t) (3.3) Realmente una de las ecuaciones es redundante ya que S+I+R=N⇒s+i+r= 1, por lo que el sistema se puede resumir en dos ecuaciones diferenciales. A pesar de ello, el modelo SIRS no se puede resolver anal´ıticamente en el caso general. La predicci´on m´as significativa de estas ecuaciones es la existencia de un epidemic threshold o umbral epid´emico[21]. Es decir, existe una cierta combinaci´on de las condiciones iniciales que determina si es posible o no que se produzca una epidemia. Ve´amoslo. Tomemos como condiciones iniciales s(0) = s0>0, i(0) = i0>0 y r(0) = 0. As´ı, en t= 0 la ecuaci´on (3.2) queda di dtt=0 =i0(βs0−µ) Para que se pueda producir una epidemia ser´a necesario que, como m´ınimo, inicialmente el n´umero de infectados se incremente, es decir di dtt=0 >0⇒βs0−µ > 0 ⇒s0>µ β(3.4) Pero s≤1⇒s0≤1. Teniendo en cuenta esta condici´on, 1 ≥s0>µ β, por lo que λ≡β µ>1 (3.5) Esto nos permite definir el umbral epid´emico como λc≡1 de forma que si λ>λcla propagaci´on de la epidemia es posible, mientras que si λ<λcel brote inicial desaparece antes de infectar a una cantidad no despreciable de la poblaci´on. cPodemos considerar que los resultados que vamos a obtener son los que conseguir´ıamos si repti´esemos el proceso muchas veces con las mismas condiciones iniciales y despu´es hicieramos el promedio[1]. Difusi´on de epidemias 16 Volviendo a la ecuaci´on (3.4) definimos R0≡s0β µ=s0λ(3.6) de forma que si R0>1 la epidemia es posible y en caso contrario desaparece. En epidemiolog´ıa cl´asica a esta cantidad se le denomina basic reproductive number y da cuenta del n´umero de infecciones secundarias que se producen cuando se introduce un individuo infectado en una poblaci´on donde todos son susceptibles. Esta condici´on es m´as restrictiva que la del umbral epid´emico, aunque como normalmente s0≈1 podemos considerar que R0≈λ[22]. Deteng´amonos un momento en esta expresi´on. Los par´ametros βyµdependen principalmente de la enfermedad, por lo que ser´a dif´ıcil modificarlos. Sin embargo, si conseguimos reducir el valor de s0podemos llegar a evitar que se produzca la epidemia, cosa que podemos realizar mediante la vacunaci´on. Veamos un par de ejemplos reales. La viruela tiene aproximadamente un umbral epid´emico λ= 5. Si queremos eliminar la enfermedad, necesitaremos que R0<1⇒s0<1/λ = 0.2. As´ı, si conseguimos vacunar al 80% de la poblaci´on seremos capaces de eliminar la enfermedad como, de hecho, se logr´o en los a˜nos 70. La rubeola, en cambio, tiene un umbral epid´emico de λ= 18, lo que implica que necesitamos vacunar al ∼94% de la poblaci´on, algo dif´ıcil pero a simple vista no imposible. El problema de ´esta y otras enfermedades con valores tan altos del umbral epid´emico es que en toda poblaci´on hay una serie de individuos que no pueden recibir la vacuna: al´ergicos, enfermos, personas con el sistema inmune debilitado, ancianos... Y en ocasiones son los suficientes como para que no se pueda reducir el valor de s0hasta los niveles que previenen la epidemia [18]. 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 r =0.0 =0.1 =0.3 =0.6 (a) Inicialmente s0= 0.99, i0= 0.01 y r0= 0. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 3 4 5 6 7 8 9 10 r λ γ=0.1 γ=0.3 γ=0.6 (b) El 80% de la poblaci´on est´a vacunada. Figura 6: Soluci´on num´erica de las ecuaciones (3.1), (3.2) y (3.3). Se representa la fracci´on de poblaci´on recuperada de la enfermedad en estado estacionario para varios valores de λyγmanteniendo µ= 0.2. La condici´on necesaria para que r6= 0 es que i(t→ ∞)6= 0 si γ6= 0, mientras que si γ= 0 basta con que se haya producido un brote epid´emico. Para finalizar con el an´alisis de la hip´otesis de homogeneous mixing podemos comprobar la validez de estas predicciones resolviendo num´ericamente el modelo, figura 6. La expresi´on del umbral epid´emico que hemos obtenido depende ´unicamente de βyµ, por lo que deber´ıamos obtener λc= 1 independientemente del valor de γsiempre que s0≈1. En la figura 6a se aprecia claramente la validez de este resultado. Cabe resaltar que en el caso en el que γ= 0 se recupera el modelo SIR, de ah´ı que exista tanta diferencia entre dicha curva y el resto, aunque vemos que el umbral es el mismo que en el modelo SIRS. Si consider´asemos el modelo SIS tambi´en obtendr´ıamos el mismo resultado, ya que el an´alisis se ha efectuado sobre la ecuaci´on (3.2) y ´esta no depende de la presencia de la fase inmune. En la figura 6b se puede comprobar c´omo cambian los resultados al partir de una situaci´on en la que solo el ∼20% de la poblaci´on es susceptible de contraer la enfermedad. Bajo estas condiciones ya no es Difusi´on de epidemias 17 v´alida la aproximaci´on R0≈λ, por lo que λcdeber´a ser distinto de uno. Seg´un la ecuaci´on (3.6) para que R0>1 ser´a necesario que λ > 5, es decir, λc= 5. Vemos claramente como los resultados num´ericos est´an de acuerdo con esta condici´on. 3.3 Difusi´on en redes scale-free La suposici´on de que es posible entrar en contacto con cualquier persona de una poblaci´on contradice claramente nuestra experiencia diaria. En el mundo real, la mayor´ıa de las personas tienen un c´ırculo de conocidos con el que se relaciona normalmente: familia, amigos, compa˜neros de trabajo... Matem´aticamente podemos representar al conjunto de contactos potenciales de una persona mediante una red. Como veremos a continuaci´on, la estructura de dicha red puede tener un gran efecto sobre la forma en la que se difunde la epidemia. La diferencia principal respecto al caso de homogeneous mixing es que un nodo solo podr´a infectarse si uno de sus vecinos est´a infectado. Debido a la heterogeneidad de la red esto implica que hay que considerar la evoluci´on de cada nodo por separado. Para poder analizar este sistema anal´ıticamente vamos a plantear una hip´otesis de campo medio: los nodos que poseen el mismo grado se comportan de la misma manera[23][2]. Por esta hip´otesis, si el estado de un nodo viene dado por itendremos que i=jsi el grado de los nodos iyjes el mismo. Denotando por Dkal conjunto de nodos de grado k, esto es equivalente a decir que i, j ∈Dk. En otras palabras, si la probabilidad de que el nodo ipertenezca al estado susceptible es si(t), la probabilidad de que los nodos de grado kpertenezcan al estado susceptible ser´a sk(t) = 1 NkX i∈Dk si(t) = si(t) por lo que en vez de plantear las ecuaciones de evoluci´on de cada nodo, basta con plantear una por grado. Denotemos por Θ(t) a la probabilidad de que un link de un nodo apunte a un nodo infectado. Un nodo de grado ktendr´a klinks, de forma que la probabilidad de que est´e en contacto con uno infectado ser´a kΘ(t). As´ı, tendremos que dsk(t) dt =−βsk(t)kΘ(t) + γrk(t) (3.7) dik(t) dt =βsk(t)kΘ(t)−µik(t) (3.8) drk(t) dt =µik(t)−γrk(t) (3.9) donde sk+ik+rk= 1. Falta determinar una expresi´on para Θ(t). La probabilidad de que un link apunte a un nodo con k links ser´a proporcional a kP(k). Por tanto, la probabilidad de que un link est´e unido a un nodo infectado ser´a Θ(t) = PkkP(k)ρk(t) PssP(s)=PkkP(k)ρk(t) hki(3.10) En esta aproximaci´on estamos considerando que la probabilidad de que un link se una a un nodo infectado es independiente del nodo de partida. Existen modelos m´as detallados en los que se tiene en cuenta la probabilidad condicional P(k|k0) de que un nodo de grado k0se una con uno de grado k. Sin embargo, esta aproximaci´on es suficiente para obtener las condiciones que determinan si se va a producir o no una epidemia[24]. Difusi´on de epidemias 18 Para determinar el punto cr´ıtico que diferencia el estado epid´emico del no epid´emico vamos a suponer que r≈0. Esto es cierto en el modelo SIS para todo tiempo t, mientras que en el SIR y en el SIRS lo es solo al inicio del proceso. No obstante, los resultados ser´an v´alidos para todos ellos ya que si no se produce un brote al inicio del proceso es imposible que se produzca despu´es. Bajo estas suposiciones la condici´on necesaria para que se pueda producir un brote es dik(t) dt t=0 >0⇒β(1 −ik(t))kΘ(t)> µik(t) para alg´un valor de k. Despejando ik, ik<βkΘ µ+βkΘ Y sustituyendo esta ecuaci´on en la definici´on de Θ (3.10), Θ>1 hkiX k k2P(k)βΘ µ+βkθ ≡f(Θ) A la vista de la figura 7 para que exista alg´un valor de Θ que cumpla la condici´on es necesario que 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Θ Θ f(Θ) f(Θ)2 Figura 7: Representaci´on gr´afica de y= Θ e y= f(Θ) con diferentes valores de kβ. df(Θ) dΘΘ=0 >1 Es decir, 1 hkiX k k2P(k)·β µ>1 Recordando que λ=β/µ llegamos a que λ > hki hk2i Por lo que el umbral epid´emico para la difusi´on en una red es λc=hki hk2i(3.11) Adem´as, de forma similar al caso anterior, podemos definir el basic reproductive number como R0=β µ hk2i hki=λhk2i hki(3.12) Ambas expresiones son muy similares a las obtenidas suponiendo homogeneous mixing, (3.5) y (3.6). Sin embargo, en esta ocasi´on tenemos t´erminos adicionales relacionados con la estructura de la red, como era de esperar. Veamos qu´e implican estos resultados. En el caso de una red aleatoria se puede demostrar[7] que hk2i=hki1−hki N Difusi´on de epidemias 19 En el l´ımite N→ ∞, por tanto, hk2i=hkipor lo que el punto cr´ıtico queda λrand =1 hki(3.13) De manera que cuanto m´as conectada est´e la red menor ser´a el punto cr´ıtico y m´as f´acil ser´a que se produzca una epidemia. No obstante, el resultado realmente interesante surge al considerar una red scale-free que, adem´as, son las redes que nos encontramos en el mundo real. Recordemos que en la secci´on 2.3 hab´ıamos determinado que en una red scale-free hk2i=∞cuando 2< α ≤3. Adem´as, en la 2.4 ve´ıamos como la mayor´ıa de las redes reales cumplen esta condici´on. Por ejemplo, la red de routers tiene un exponente α= 2,1d. Aplicando este resultado a la ecuaci´on (3.11) obtenemos algo sorprendente, λc=hki hk2i= 0 Es decir, en una red scale-free con un exponente tal que 2 < α ≤3 no existe umbral epid´emico. Por tanto, existe una probabilidad no nula de que se produzca una epidemia para cualquier valor de βyµ (siempre que β > 0). Otra manera de verlo es mediante el basic reproductive number. Por la ecuaci´on (3.12) si hk2i=∞entonces R0>1. Pensemos de nuevo en la vacunaci´on. En este caso no nos basta con reducir la poblaci´on inicial susceptible ya que el que el umbral sea cero depende de la estructura de la red. Es en este punto donde cobra sentido la discusi´on sobre la importancia de los hubs de la p´agina 10. Ve´ıamos que, aunque una red scale-free era muy resistente a la p´erdida (inmunizaci´on) de nodos aleatorios, era muy d´ebil frente a ataques dirigidos (vacunaci´on). De esta forma, lo que era una desventaja podemos convertirlo en algo positivo. Si conseguimos localizar a los grandes hubs de la red y los inmunizamos, seremos capaces de hacer que la red deje de ser conexa y detendremos la difusi´on de la epidemia. Adem´as, el eliminar a los nodos de mayor grado har´a que la diferencia entre los de menor y mayor grado no sea tan grande, incluso es posible que hk2itome un valor finito. En consecuencia, se reestablecer´a un umbral finito y ser´a posible detener la epidemia aunque la red siga siendo conexa. Para finalizar, vamos a comprobar si estas predicciones se ajustan a las simulaciones. En este caso vamos a emplear simulaciones de Monte Carlo, lo que nos va a permitir introducir el factor de aleatoriedad inherente a este tipo de procesos (recordemos que los par´ametros β,µyγson probabilidades). Estos efectos son especialmente importantes al inicio de la epidemia, cuando el n´umero de infectados es todav´ıa muy peque˜no[26]. En primer lugar empleamos una red aleatoria con hki= 4.04. Si la predicci´on de la ecuaci´on (3.13) es correcta, deber´ıamos obtener λc≈0.25. Como se puede comprobar en la figura 8a el punto cr´ıtico es precisamente ∼0.25e. En la red scale-free, a pesar de tener un grado medio similar, los resultados son bastante diferentes. De acuerdo a la ecuaci´on (3.11) idealmente deber´ıamos obtener el punto cr´ıtico en el 0, pero al usar una red finita obtenemos un resultado ligeramente diferente. La red empleada posee hk2i= 66.44, lejos de un valor infinito. No obstante, el valor es lo suficientemente grande como para introducir una gran diferencia respecto al caso anterior. As´ı, teniendo en cuenta este dCabe destacar que en la pr´actica la difusi´on de virus entre ordenadores y la de epidemias entre personas se puede tratar con los mismos modelos. La primera vez que se propuso aplicar un modelo epid´emico a una red scale-free fue al tratar de explicar la persistencia de los virus en los ordenadores[25]. eEl que la transici´on entre un estado y otro no sea abrupta se denomina efecto de tama˜no finito y es una consecuencia del tama˜no finito de la red[1]. Difusi´on de epidemias 20 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 r λ α=0.1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 r λ α=0.1 α=0.3 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 r λ α=0.1 α=0.3 α=0.5 -0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 0.22 0.23 0.24 0.25 0.26 0.27 0.28 λc α=0.1 -0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 0.22 0.23 0.24 0.25 0.26 0.27 0.28 λc α=0.1 α=0.3 -0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 0.22 0.23 0.24 0.25 0.26 0.27 0.28 λc α=0.1 α=0.3 α=0.5 (a) Modelo SIRS sobre una red aleatoria 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 r λ α=0.1 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 r λ α=0.1 α=0.3 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 r λ α=0.1 α=0.3 α=0.5 -0.01 -0.005 0 0.005 0.01 0.015 0.055 0.06 0.065 0.07 0.075 0.08 λc α=0.1 -0.01 -0.005 0 0.005 0.01 0.015 0.055 0.06 0.065 0.07 0.075 0.08 λc α=0.1 α=0.3 -0.01 -0.005 0 0.005 0.01 0.015 0.055 0.06 0.065 0.07 0.075 0.08 λc α=0.1 α=0.3 α=0.5 (b) Modelo SIRS sobre una red scale-free Figura 8: Simulaci´on MC de la difusi´on de una epidemia seg´un el modelo SIRS sobre redes de ∼10.000 nodos. La red aleatoria posee hki= 4.04, mientras que la scale-free tiene hki= 4.59 y hk2i= 66.44 par´ametro deber´ıamos obtener λc= 0.69. Podemos ver en la figura 8b como claramente este resultado concuerda con las simulaciones. A pesar de la simplicidad del modelo, estos resultados concuerdan cualitativamente con la epidemiolog´ıa de una gran cantidad de pat´ogenos. No obstante, para realizar predicciones cuantitativas que puedan emplearse para evaluar qu´e pol´ıticas aplicar ante un brote y analizar el coste-beneficio de ellas es necesario refinar este modelo. Es en este contexto en el que surgen los modelos metapoblacionales[26]. Modelos metapoblacionales 21 4 Modelos metapoblacionales Como sucede en muchos otros campos, a la hora de proponer un modelo debemos tener en cuenta que cuanto m´as sencillo sea, peor precisi´on tendr´an las predicciones y viceversa. Por otra parte, el agregar m´as par´ametros al modelo, adem´as de dificultar su manejo, posee un problema inherente: ser´a necesario contar con muchos m´as datos para poder determinarlos. Si no podemos estimar estos par´ametros con gran precisi´on no nos servir´a de nada complicar el modelo ya que sus predicciones pueden no tener menos incertidumbre que la que nos dar´ıan modelos m´as sencillos pero con par´ametros bien determinados. Figura 9: Diferentes estructuras empleadas en el estudio de la difusi´on de epidemias. Los colores representan el estado de cada nodo. De izquierda a derecha: el modelo m´as b´asico es considerar homogeneous mixing; tambi´en podemos tener en cuenta la estructura social (edad, g´enero...); para modelar la heterogeneidad de las relaciones sociales se emplean redes de contacto; en los modelos multi-escala (entre los que se incluyen los metapoblacionales) se tiene en cuenta la interacci´on entre diferentes subpoblaciones (como podr´ıan ser las ciudades) pero dentro de ellas se emplea homogeneous mixing; en los modelos agent-based se recrean los movimientos e interacciones de cada una de las personas individualmente con gran precisi´on[17]. En la figura 9 se representan esquem´aticamente los diferentes modelos que podemos emplear. Hasta ahora nos hemos limitado a los tres primeros ya que, como vamos a ver a continuaci´on, son ingredientes clave para el desarrollo de los siguientes, los modelos metapoblacionales. 4.1 La difusi´on en metapoblaciones: procesos reacci´on-difusi´on Una metapoblaci´on es una poblaci´on subdividida en comunidades discretas unidas por conexiones que representan los posibles movimientos de los individuos. Dentro de cada comunidad la din´amica de transmisi´on de las enfermedades se describe mediante los modelos compartimentales est´andar[27]. Matem´aticamente, la difusi´on en redes metapoblacionales se estudia como un proceso de reacci´ondifusi´on (RD) sobre una red compleja. As´ı, consideramos que los nodos son ciudades y que las personas se pueden mover (difundir) de una a otra a trav´es de los links. En las secci´on 2 hemos presentado todos los conceptos de la teor´ıa de redes que vamos a necesitar, pero no hemos hablado todav´ıa sobre los procesos RD. Los modelos RD se desarrollaron para explicar las reacciones qu´ımicas, aunque hoy en d´ıa tienen aplicaciones en muchos otros campos. En f´ısica, por ejemplo, podemos encontrar este tipo de procesos al estudiar la din´amica no-lineal del transporte en semiconductores[28]. A nivel microsc´opico, los procesos RD consisten en un conjunto de part´ıculas que se difunde en el espacio y que est´an sujetas a una serie de reacciones que dependen del problema que se est´e tratando (figura 10). Mientras que en los modelos RD fermi´onicos se considera que hay un l´ımite en el n´umero de part´ıculas que puede haber en cada nodo (por el principio de exclusi´on), en los modelos bos´onicos no existe tal restricci´on[29]. Para modelar la difusi´on de una epidemia vamos a considerar un proceso RD bos´onico en el cual las part´ıculas ser´an las personas. ´ Estas se dividir´an en tres categor´ıas (susceptible, infectado o inmune) pudi- Modelos metapoblacionales 22 endo pasar de una a otra mediante una “reacci´on”. Cabe destacar que el modelo que vamos a desarrollar es determinista, por lo que no es del todo v´alido para estudiar el inicio de la epidemia (en dicho punto las fluctuaciones estoc´asticas no son despreciables)[30]. Figura 10: Esquema de un proceso RD aplicado en una red. Consideremos, por simplicidad, que podemos describir la epidemia mediante el modelo SIR. Las ecuaciones de las reacciones ser´an: I+Sβ/N −−−→ 2I Iµ −−−→ R Reacciones que se producir´an en cada poblaci´on j. Bajo la hip´otesis de homogeneous mixing (secci´on 3.2) la tasa de infecci´on en cada nodo jla podemos denotar por βΓjdonde Γj=IjSj Nj Bajo la aproximaci´on de campo medio (secci´on 3.3) tendremos que Ik=1 VkXIj, Sk=1 VkXSj siendo Vkel n´umero de nodos de grado k. En esta aproximaci´on la tasa de infecci´on ser´a βΓk=βIkSk/Nk. Podemos escribir la variaci´on del n´umero de infectados en cada subpoblaci´on de grado kmediante una ecuaci´on maestra, Ik(t+ ∆t)−Ik(t) = W+ k−W− k donde W± krepresenta el n´umero de infectados que entran/abandonan los nodos de grado ktanto por contagio como por difusi´on. Si suponemos que el ritmo de difusi´on depende del grado del nodo, pk, tendremos que W− k=pkIk+ (1 −pk)µIk El primer t´ermino representa el n´umero de individuos que abandonan el nodo y el segundo el n´umero que se cura de los que se quedan. De forma similar, W+ k= (1 −pk)βΓk+kX k0 P(k0|k)dk0k[(1 −µ)Ik0+βΓk0] En este caso el primer t´ermino da cuenta del n´umero de gente que se queda en su nodo y se infecta y el segundo el n´umero de infectados que llegan de otros nodos. En efecto, dk0kes la fracci´on de infectados que van de un nodo de grado k0a uno de grado k. El n´umero de infectados del nodo de grado k0ser´a [(1 −µ)Ik0+βΓk0] (los que no se recuperan m´as los nuevos infectados). El otro t´ermino del sumatorio, P(k0|k), es la probabilidad de que haya un link entre un nodo de grado ky otro de grado k0. As´ı, al sumar sobre todos los grados, obtendremos el n´umero total de infectados que pueden llegar a trav´es de un link al nodo de grado k. Por ´ultimo, multiplicamos esta cantidad por el n´umero de links del nodo, k. Introduciendo estas definiciones llegamos a la ecuaci´on diferencial que determina la evoluci´on del n´umero de infectados presentes en los nodos de grado k, ∂tIk=−pkIk+ (1 −pk)[−µIk+βΓk] + kX k0 P(k0|k)dk0k[(1 −µ)Ik0+βΓk0] De forma an´aloga se pueden obtener las expresiones para la evoluci´on de SkyRk. Con estas ecuaciones quedar´ıa determinada la evoluci´on del sistema dadas las condiciones iniciales adecuadas. Modelos metapoblacionales 23 4.2 Global invasion threshold En los modelos metapoblacionales, al igual que en los estudiados anteriormente, uno de los principales objetivos es hallar qu´e condiciones determinan que se pueda producir una epidemia. Es posible expresar estas condiciones mediante un ´unico par´ametro, R∗, que da cuenta del n´umero medio de subpoblaciones a las que una subpoblaci´on infectada puede transmitir la enfermedad mediante la movilidad de individuos infectados durante la duraci´on de la epidemia. Este par´ametro se conoce como global invasion threshold[31]. As´ı pues, si R∗<1 la epidemia permanece en la subpoblaci´on inicial hasta que desaparece, mientras que si R∗>1 se puede propagar espacialmente en el sistema y alcanzar una dimensi´on global. Como vemos, este par´ametro es una extensi´on del basic reproductive number de las secciones 3.2 y 3.3. Las ecuaciones que acabamos de deducir no nos van a servir para calcular R∗ya que no tienen en cuenta ni las fluctuaciones inherentes al proceso de difusi´on ni la probabilidad de extinci´on del brote epid´emico[30]. Por esta raz´on, vamos a tener que desarrollar otro modelo en el que se tenga en cuenta la aleatoriedad presente en el proceso. Figura 11: Esquema del modelo SIR sobre una metapoblaci´on[30]. Siguiendo con el modelo SIR, supongamos que la epidemia empieza en una subpoblaci´on de grado ky Nkpersonas en la que R0>1, de forma que localmente es posible que se produzca un brote. Denotemos por αNkal n´umero total de individuos infectados que se generan en esta subpoblaci´on. Cada persona infectada permanecer´a de media en dicho estado durante un tiempo 1/µ, durante el cual podr´a viajar a una subpoblaci´on vecina de grado k0con una probabilidad dkk0. As´ı, en primera aproximaci´on, podemos considerar que el n´umero de nuevos infectados que llegan a una subpoblaci´on de grado k0provenientes de la subpoblaci´on original es λkk0=dkk0 αNk µ Sea D0 kel n´umero de subpoblaciones “infectadas” de grado kal principio del proceso (generaci´on 0). Cada una de ellas podr´a comunicar la infecci´on a sus vecinas, determinando el valor de D1 k. Es decir, el n´umero de subpoblaciones infectadas en la siguiente generaci´on. Al inicio del proceso el n´umero de poblaciones en las que hay un brote ser´a lo suficientemente peque˜no como para poder considerar un proceso tipo ´arbol de forma que Dn kpodremos relacionarlo con Dn−1 k. As´ı, Dn k=X k0 Dn−1 k0(k0−1)λk0k(R0−1)P(k|k0) 1−Dn−1 k Vk!(4.1) Analicemos detenidamente esta expresi´on, - Cada subpoblaci´on de grado k0podr´a infectar a k0−1 subpoblaciones. El factor −1 da cuenta del hecho de que no se puede traspasar la epidemia a la subpoblaci´on que le ha infectado. Modelos metapoblacionales 24 - Como hay Dn−1 k0subpoblaciones de grado k0infectadas, el n´umero total de subpoblaciones a las que puede pasar la epidemia desde una de grado k0ser´a Dn−1 k0(k0−1). - A continuaci´on multiplicamos por la probabilidad de tengan un link con un vecino de grado k,P(k|k0) - Suponiendo que R0−11 podemos aproximar la probabilidad de que se produzca un brote en una subpoblaci´on de grado kpor λkk0(R0−1) - Por ´ultimo, debemos comprobar que la subpoblaci´on de grado kno est´a infectada todav´ıa, (1 − Dn−1 k/Vk) Sumando las contribuciones de todas las subpoblaciones obtenemos el n´umero de subpoblaciones de grado kinfectadas en la siguiente generaci´on. Al estar al comienzo de la epidemia la probabilidad de que haya subpoblaciones infectadas es muy peque˜na y podremos considerar (1 −Dn−1 k/Vk)≈1. El siguiente paso es definir c´omo se produce la difusi´on. Supongamos que en cada generaci´on la fracci´on de individuos que abandona cada subpoblaci´on es p. L´ogicamente, la fracci´on de personas que viaja por cada link ser´a p/k donde kes el n´umero de links del nodo (difusi´on homog´enea). As´ı, λk0k=p k0 αNk0 µ Si la red no posee correlaciones entonces P(k|k0) = kP(k)/hki[1]. Adem´as, en difusi´on homog´enea se cumple que Nk0=¯ Nk0/hki, donde ¯ Nes el tama˜no medio de las subpoblaciones. Introduciendo estas expresiones en (4.1), Dn k= (R0−1)pα ¯ NkP (k) µhki2X k0 Dn−1 k0k0(k0−1) Definiendo Θn≡Pk0Dn k0k0(k0−1) podemos reescribir esta expresi´on como Θn= (R0−1)hk2i−hki hki2 p¯ Nα µΘn−1 Para que el n´umero de subpoblaciones infectadas se incremente ser´a necesario que R∗≡(R0−1)hk2i−hki hki2 p¯ Nα µ>1 (4.2) Cualitativamente este resultado es muy similar al obtenido al estudiar la difusi´on en redes scale-free (3.12). Cuanto mayor sea la heterogeneidad de la red (mayor hk2i) m´as f´acil ser´a que se cumpla la condici´on R∗>1. Es m´as, en el caso de una red scale-free ideal R∗→ ∞, el mismo resultado que en la difusi´on en una poblaci´on. Resulta interesante escribir la condici´on en funci´on de p. Si R0−11 la constante αse puede aproximar en el modelo SIR por[21] α≈2(R0−1) R2 0 De forma que p¯ N≥hki2 hk2i−hki µR2 0 2(R0−1)2(4.3) La interpretaci´on de este resultado es inmediata. Si existe una gran movilidad (pmuy elevado) ser´a f´acil que se cumpla la condici´on y se producir´a una epidemia a escala global. An´alogamente, un valor Estudio de la persistencia en metapoblaciones 31 As´ı, podemos considerar que en cualquier punto (p,γ−1) que est´e por encima de la l´ınea roja la epidemia se ha extinguido. De la misma forma, en cualquier punto que est´e por debajo de la verde la epidemia permanece siempre. En el l´ımite en el que p→0 el sistema equivaldr´ıa a un mont´on de poblaciones desconectadas entre s´ı. Vemos como este comportamiento se mantiene as´ı hasta que se alcanza un cierto valor cr´ıtico de la movilidad en el que empieza a haber una fracci´on apreciable de intercambio de individuos entre las poblaciones. Este intercambio hace que, aunque el periodo de inmunidad sea muy elevado, haya un suministro continuo de individuos susceptibles entre las subpoblaciones, lo que permite a la epidemia persistir. Sin embargo, se alcanza un valor m´aximo a partir del cual la tendencia se invierte. Esto se puede achacar a que el intercambio es tan grande que en pr´acticamente todas las subpoblaciones habr´a un outbreak. As´ı, la epidemia se propagar´a muy rapidamente y descender´a el n´umero de susceptibles en todas las poblaciones, por lo que un valor peque˜no de γ−1ser´a suficiente para eliminar la epidemia. A valores muy elevados de la movilidad los efectos de la red se deber´ıan ver disminuidos, recuperando en el l´ımite p→1 los resultados de una ´unica poblaci´on con 108individuos. 0 100 200 300 400 500 600 700 800 900 1000 1100 1e-06 1e-05 0.0001 0.001 0.01 0.1 1 -1 p SF ER Figura 19: Persistencia frente a duraci´on de la inmunizaci´on y movilidad en las redes SF y ER de 1000 nodos. En efecto, tal y como se ve en la figura 19 a movilidades bajas y altas recuperamos el mismo resultado en la red SF y en la ER. Adem´as, la m´axima inmunidad es mayor en la red SF que en la ER, por lo que tiene que estar relacionado con la estructura de la red. Ser´a m´as f´acil que la epidemia persista en un hub ya que al tener muchos vecinos tendr´a un gran intercambio de individuos. As´ı, por una parte recibir´a una gran cantidad de susceptibles y por otra propagar´a la epidemia a muchas subpoblaciones. Se recuperan resultados an´alogos al comparar las redes de 100 y de 10000. Por ´ultimo, comparamos los resultados para las redes SF, figura 20. La poblaci´on total en todas ellas es de 108, por lo que en las m´as grandes la poblaci´on de cada nodo ser´a menor. As´ı, cuando la epidemia est´a localizada en unas pocas subpoblaciones (movilidad baja) la persistencia es menor en las redes grandes. Al tener menor poblaci´on se reducir´a antes el n´umero de susceptibles, de forma que ser´a necesario un valor peque˜no de γ−1para que la epidemia pueda mantenerse. Estudio de la persistencia en metapoblaciones 32 0 200 400 600 800 1000 1200 1400 1600 1800 2000 1e-06 1e-05 0.0001 0.001 0.01 0.1 1 -1 p N=10000 N=1000 N=100 Figura 20: Persistencia frente a duraci´on de la inmunizaci´on y movilidad en las redes SF. Para movilidades elevadas, como coment´abamos anteriormente, se pierden los efectos de la red ya que se tiende a la situaci´on de una ´unica poblaci´on con 108individuos. Estos resultados son muy interesantes ya que abren nuevos mecanismos para combatir una epidemia. As´ı, aunque a simple vista pudiera parecer contradictorio, podr´ıamos aumentar las probabilidades de erradicar una epidemia incrementando el intercambio de personas sea cual sea la red. Esta medida, l´ogicamente, ser´ıa mucho m´as rentable en t´erminos econ´omicos y log´ısticos que restringir el tr´afico[38]. Como coment´abamos anteriormente estos resultados abren la puerta a nuevos estudios te´oricos. En concreto, resultar´ıa interesante determinar los tres puntos cr´ıticos de la movilidad en funci´on de la estructura de la red. Conclusiones 33 6 Conclusiones A lo largo del trabajo hemos podido comprobar la presencia de los comportamientos emergentes que mencion´abamos en la introducci´on y su estrecha relaci´on con la estructura de las redes complejas: - En el apartado de introducci´on a las redes hemos visto que la mayor´ıa de ellas siguen un comportamiento libre de escala. - Al aplicar los modelos epidemiol´ogicos a una red hemos podido ver el primer comportamiento colectivo consecuencia de la interacci´on entre los individuos: el umbral epid´emico desaparece en las redes libres de escala. - En los modelos metapoblacionales vuelve a aparecer un ejemplo de comportamiento colectivo: el global invasion threshold depende expl´ıcitamente de la red y de la intensidad de la interacci´on (la movilidad) entre las subpoblaciones. - En el estudio de la persistencia hemos comprobado la gran dependencia que existe con la movilidad y la estructura de la red. Podemos concluir, por tanto, que la difusi´on de una epidemia es un ejemplo claro de lo que hemos denominado sistemas complejos. As´ı, si queremos ser capaces de predecir e incluso detener la propagaci´on de una epidemia no nos bastar´a con estudiar las caracter´ısticas de la enfermedad, sino que deberemos conocer con precisi´on los patrones de movilidad y las estructuras sociales de la poblaci´on, tal y como se ha podido comprobar en el estudio de la persistencia en metapoblaciones. Bibliograf´ıa 34 Bibliograf´ıa [1] Mark Newman. Networks: an introduction. Oxford University Press, 2010. [2] S G´omez, A Arenas, J Borge-Holthoefer, S Meloni, and Y Moreno. Discrete-time Markov chain approach to contact-based disease spreading in complex networks. EPL (Europhysics Letters), 89(3):38009, 2010. [3] B´ela Bollob´as. Modern graph theory, volume 184. Springer, 1998. [4] Supriyo Datta. Electronic transport in mesoscopic systems. Cambridge University press, 1997. [5] Erd¨os P. and R´enyi A. On random graphs i. Publ. Math. Debrecen, 6:290–297, 1959. [6] Statista. Average number of facebook friends of u.s. users in 2014, by age group. http://www. statista.com/statistics/232499/. [7] Albert-L´aszl´o Barab´asi and Jennifer Frangos. Linked: The New Science Of Networks Science Of Networks. Basic Books, 2002. [8] R´eka Albert, Hawoong Jeong, and Albert-L´aszl´o Barab´asi. Internet: Diameter of the world-wide web. Nature, 401(6749):130–131, 1999. [9] Albert-L´aszl´o Barab´asi and R´eka Albert. Emergence of scaling in random networks. Science, 286(5439):509–512, 1999. [10] John MD Coey. Magnetism and magnetic materials. Cambridge University Press, 2010. [11] Size of the world wide web. http://www.worldwidewebsize.com/. [12] Albert-L´aszl´o Barab´asi and Eric Bonabeau. Scale-free networks. Scientific American, 2003. [13] Guido Caldarelli. Scale-free networks: complex webs in nature and technology. OUP Catalogue, 2007. [14] Mapping the internet. http://nicolasrapp.com/?p=1180. [15] Manfred S Green, Tiberio Swartz, Elana Mayshar, Boaz Lev, Alex Leventhal, Paul E Slater, and Joshua Shemer. When is an epidemic an epidemic? The Israel Medical Association Journal: IMAJ, 4(1):3–6, 2002. [16] David Easley and Jon Kleinberg. Networks, crowds, and markets. Cambridge Univ Press, 6(1):6–1, 2010. [17] Alain Barrat, Marc Barthelemy, and Alessandro Vespignani. Dynamical processes on complex networks, volume 574. Cambridge University Press Cambridge, 2008. [18] The Pennsylvania State University. Epidemic - the dynamics of infectious diseases. Coursera, 2013. [19] Albert-L´aszl´o Barab´asi. Class 17: Epidemic modeling. University Lecture, 2012. [20] Stefano Boccaletti, Vito Latora, Yamir Moreno, Mario Chavez, and D-U Hwang. Complex networks: Structure and dynamics. Physics reports, 424(4):175–308, 2006. [21] James D. Murray. Mathematical Biology: I. An Introduction (Interdisciplinary Applied Mathematics) (Pt. 1). Springer, 2007. [22] James Holland Jones. Notes on R0. Standford University, May 2007. [23] Beniamino Guerra and Jes´us G´omez-Garde˜nes. Annealed and mean-field formulations of disease dynamics on static and adaptive networks. Physical Review E, 82(3):035101, 2010. Bibliograf´ıa 35 [24] Yamir Moreno, Romualdo Pastor-Satorras, and Alessandro Vespignani. Epidemic outbreaks in complex heterogeneous networks. The European Physical Journal B-Condensed Matter and Complex Systems, 26(4):521–529, 2002. [25] Romualdo Pastor-Satorras and Alessandro Vespignani. Epidemic spreading in scale-free networks. Physical review letters, 86(14):3200, 2001. [26] Neil M Ferguson, Matt J Keeling, W John Edmunds, Raymond Gani, Bryan T Grenfell, Roy M Anderson, and Steve Leach. Planning for smallpox outbreaks. Nature, 425(6959):681–685, 2003. [27] Andrea Apolloni, Chiara Poletto, Jos´e J Ramasco, Pablo Jensen, Vittoria Colizza, et al. Metapopulation epidemic models with heterogeneous mixing and travel behaviour. Theoretical Biology and Medical Modelling, 11(1):3, 2014. [28] Eckehard Sch¨oll. Nonlinear spatio-temporal dynamics and chaos in semiconductors, volume 10. Cambridge University Press, 2001. [29] Vittoria Colizza, Romualdo Pastor-Satorras, and Alessandro Vespignani. Reaction–diffusion processes and metapopulation models in heterogeneous networks. Nature Physics, 3(4):276–282, 2007. [30] Vittoria Colizza and Alessandro Vespignani. Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: Theory and simulations. Journal of theoretical biology, 251(3):450– 467, 2008. [31] Sandro Meloni, Nicola Perra, Alex Arenas, Sergio G´omez, Yamir Moreno, and Alessandro Vespignani. Modeling human mobility responses to the large-scale spreading of infectious diseases. Scientific reports, 1, 2011. [32] T D´eirdre Hollingsworth, Neil M Ferguson, and Roy M Anderson. Will travel restrictions control the international spread of pandemic influenza? Nature medicine, 12(5):497–499, 2006. [33] Ben S Cooper, Richard J Pitman, W John Edmunds, and Nigel J Gay. Delaying the international spread of pandemic influenza. PLoS Medicine, 3(6):e212, 2006. [34] Michele Tizzoni, Paolo Bajardi, Chiara Poletto, Jos´e J Ramasco, Duygu Balcan, Bruno Gon¸calves, Nicola Perra, Vittoria Colizza, and Alessandro Vespignani. Real-time numerical forecast of global epidemic spreading: case study of 2009 a/h1n1pdm. BMC medicine, 10(1):165, 2012. [35] Global epidemic and mobility (GLEAM) computational model. http://www.mobs-lab.org/ global-epidemic-and-mobility-model.html. [36] SocioEconomic Data and Application Center at Columbia University. Gridded population of the world (GPW). http://sedac.ciesin.columbia.edu/data/collection/gpw-v3. [37] Michele Catanzaro, Mari´an Bogu˜n´a, and Romualdo Pastor-Satorras. Generation of uncorrelated random scale-free networks. Physical Review E, 71(2):027103, 2005. [38] Alberto Aleta, Andreia Hisi, Chiara Poletto, Sandro Meloni, Vittoria Colizza, and Yamir Moreno. in preparation.