Avian Influenza Type A-H5N1 Epidemiological Model: Puerto Rico as a Case Study

Keywords

SIR; Avian Influenza; Epidemiology; Metapopulation; Vaccination; Basic Reproductive Number; Simulations Spread of Diseases

Abstract

Our research focused on Avian Influenza Type A-H5N1, specifically on an epidemiological model centered in Puerto Rico. Our main goal is to address the following: first, to determine the potential outbreaks of this disease in Puerto Rico using as a base the location of the poultry industry as a hub, we are interested in the repercussions of the infection among the human-to-human potential interaction. The second goal centers on the possibility of vaccination to mitigate an epidemic among humans. In order to address these goals and future ones, we will construct a mathematical model and use parameters according to two cases; the first is a single population model and the second one is a metapopulation model involving 5 cities in Puerto Rico. Our research will specifically target the spread of this particular disease, to investigate possible alternatives to mitigate the spread using measures of immunization. Our results show that a 30% vaccination regime will eradicate the disease in cities that are immunized.

Citation

Collazo-Rivera M and Cruz-Aponte M. Avian Influenza Type A-H5N1 Epidemiological Model: Puerto Rico as a Case Study. SM Vaccine Vaccin. 2015;1(1):1005.

Introduction

There are three types of influenza viruses: A, B and C. They are divided into subtypes on the basis of two proteins on the surface of the virus: Hemagglutinin (HA) that has 17 subtypes known and Neuraminidase (NA) that has 10 known subtypes [1]. Epidemiologists need to make an educated guess for the inclusion of such variants in annual vaccines in order to assure proper immunization of the population and mitigate a possible epidemic among the population. However, due to the lack of immunity in humans against new mutations of the virus, epidemics or even pandemics can emerge resulting in high morbidity and mortality [2]. In this article we focus on Avian Influenza Type A-H5N1 that is a zoonotic disease (i.e. spread from animal to humans) it is an acute and recurring respiratory disease, occurring in particular during winter months and straining the public health system worldwide. It can only spread from infected poultry to humans that had been in contact with poultry or infected soil. Avian influenza has high mortality rate, as high as 60% [2-6] hence it is a mayor concern if it mutates and can spread from human to human. As any type of influenza it is fatal for immune compromised individuals; children, elderly, patients with chronic illnesses and pregnant women. The initial symptoms are: high temperature 38oC or 100.4oF, upper respiratory tract symptoms, diarrhea, vomiting, abdominal pain, inflammation of the lungs (pleuritic pain) and nose bleeding [3-7].

We implemented modification of a simpler SIR model (Susceptible-Infected-Recovered) developed by Kermack and McKendrick [8], then implemented an SEIR (adding Exposed individuals) for a single population to study the effect of different levels of vaccination in the population that can potentially contain or eradicate the disease in a hypothetical scenario of an epidemic in Puerto Rico. We will base our models in past epidemics for different countries that have suffered from epidemics or outbreaks of Avian Influenza Type A-H5N1 [2]. In 1959 chickens spread the disease in Scotland, in 1992 the disease was spread in England by turkeys, in 1997 chickens also spread the disease in Hong Kong [3-6]. In different countries in Asia and the Middle East, between 2005 and 2006 the avian species that spread the disease predominantly were domestic and wild poultry [3-6].

This article will specifically target the spread of infectious diseases, which will be held in perspective with a mathematical-epidemiological model to determine possible alternatives to mitigate the potential spread of Avian Influenza Type A-H5N1 using measures such as immunization. Particularly, we address the spread of influenza type A-H5N1 in the country of Puerto Rico in the town of Cayey, using as a possible site of infection the city of Aibonito. The hub of infection will be set in Cayey where there is a poultry factory located near. We would also work on constructing a metapopulation model involving 5 cities of Puerto Rico connected by the main highways. A meta population model consists of a group of interacting spatially separated populations of the same species, it is defined as a set of differential equation coupled together. After our first simple model we implemented a SIR-type mobility model for five cities in Puerto Rico to investigate the potential outbreak of this disease in the linked areas of: San Juan, Caguas, Cayey, Ponce and Humacao as seen in Figure 1 in the methodology section.

Our focus is on the behavior of the epidemic of human interaction and the effect on the entire population. However, we implemented the mathematical model to address these specific aims: first, to determine the outbreaks of this disease in Puerto Rico using as a base the location of the poultry industry as a hub. The second goal centers on the possibility of vaccination to mitigate an epidemic among humans. Third, for the meta population networking approach we want to determine the potential outbreaks of this disease in Puerto Rico and how different vaccination implementations mitigate or reshape the epidemic outbreak as time goes by.

The research questions we want to address are: What mechanisms are effective to contain a potential Avian A-H5N1 epidemic? Which dynamics of disease spreading affect an outbreak among cities? What mechanisms are effective to contain a potential avian flu epidemic among cities? In order to address this question and further ones we will construct a mathematical model and use a different set of parameters, for now we are focusing first on the mitigation of the epidemic using vaccination strategies. Important to mention, we are interested in the repercussion of the infection among human-to-human interactions for the networking approach.

Methodology

We implemented the modification of a simpler SIR epidemiological model evolving to a SEIRV and use parameter values of outbreaks that occurred in different countries in the past to modify them to the case of Puerto Rico [9,10]. We want to explore what would happen in Puerto Rico if the disease of Avian Influenza Type A-H5N1 emerged, taking into consideration what happened in different countries such as: Europe, Latin America and Middle East [11]. In order to present our simple SEIRV mathematical model, let’s define our epidemiological classes as shown in Table 1 and parameters as shown in Table 2.

Table 1: Epidemiological status of humans and poultry

Table 2: Parameter values of the mathematical model.

We show in Figure 1 a schematic of the transition of individuals into the different epidemiological classes. Then we present our system of ordinary differential equations and our simulations. Afterwards we implemented a meta population model modification of a simpler mobility model for five cities in Puerto Rico on a SIR-type model, for future work we will consider an SIRV network model [12-14]. In this meta population approach we want to address the potential outbreak of this disease in the areas of: San Juan, Caguas, Cayey, Ponce and Humacao that are linked as shown in Figure 1.

Figure 1: Google Map of Puerto Rico showing the traffic flow of the 5 connected cities studied in our model. (https://www.google.com.pr/ maps/@18.1986249,-66.5863458,9z?hl=en)

SEIRV single city epidemiological model

Let’s describe our modeling approach, starting with the description of the model shown in Figure 2 and the equations that follow. The model classifies the subjects into the different compartments according to species and disease stages. The model classifies the population into a unique epidemiological class or compartment in general as described in Table 1: susceptible subjects (Si ), exposed individuals (Ei ), infected (Ii ), recovered (Ri ) and vaccinated (Vi ) individuals. Ni is the total population; where i∈{p,h} the subscripts p and h indicate poultry or human subjects respectively. Here we assume that the total population of poultry is 60,000 and the population of the city of Cayey (that we are using for a one city model) is 48,119 according to the 2010 census data.

Figure 2: Schematic of the mathematical model for Avian Influenza A-H5N1.

The parameters used in the mathematical model for Avian Influenza Type A-H5N1 are shown in Table 2 in general (recall that subscripts of p are h are imposed to identify poultry or humans respectively) they are defined as:µi that represent the demographics which is defined as the rate of birth and death for chickens and b is the “birth rate” when the chickens are massively produced and distributed to factories, the βi parameter is the ratio of effective contact between the interaction of a susceptible and an infected subject, the recovery rate is αi where (1/αi ) is the number of days the person was ill, 1/σI is the incubation period and πh is the vaccination rate for humans.

We present the schematic flow of the model on Figure 2 and below we present the ordinary differential equations that define them. Each box represents a compartment where subjects are categorized into the epidemiological classes described in Table 1 at the rates described and shown in Table 2. The differential equation system is as follows:

dSp /dt = – βp Sp Ip /Np + bNp + µNp + µSp

dIp /dt = βp Sp Ip /Np + αp Ip– δp Ip - µIp

dRp /dt = αp Ip – µRp

dSh /dt = – (γh βh Sh Eh + βh Sh Ih + βhp Sh Ip ) / Nh – πh Sh

dEh /dt = (γh βh Sh Eh + βh Sh Ih + βhp Sh Ip ) / Nh – σh Eh– πh Eh

dIh /dt = σh Eh – αh Ih– δh Ih

dRh /dt = αh Ih – πh Rh

dVh /dt = πh (Sh + Eh + Rh )

An arrow entering a compartment indicates inflow of individuals (or poultry), to address this in the equations this quality has a plus sign (+). Similarly an arrow leaving a compartment indicates that the sign is negative (-).

SIR epidemiological mathematical model with respect to the poultry: Entering Sp are susceptible chickens, the parameter b or µ, indicates (with a + sign) the demographics (birth rate µNp and death rate µSp ) of chickens on a free range environment or import of massive amounts of chickens to a factory indicated with the factor bNp . The factors -µSp , -µIp , -µRp are the natural death of the poultry. The susceptible chickens that get into the Ip class are indicated with the factor -βp Sp Ip /Np with a - sign is going out of Sp into the Ip class it represents the contact between a susceptible and an infected subject, hence it has a + sign getting into the infected poultry Ip . The parameter αp Ip indicates the recovery ratio of the poultry with a - sign getting out of the infected class and a + sign getting into the recovery class. Chickens that died from the disease are indicated with the factor -δ_p I_p .

SEIRV epidemiological mathematical model with respect to the human population: In this model since we are running the simulations and considering the time of the epidemic for less than a year we are not including demographics for the human population. The most important factor in this model is the force of infection:

(γ_h β_h S_h E_h +β_h S_h I_h +β_hp S_h I_p ) / N_h

that represents how a susceptible human gets infected with A-H5N1 by means of interacting with exposed individuals (γ_h β_h S_h E_h ) /Nh that can transmit the disease at a lower rate than the interaction with infected individuals βh Sh Ih /Nh or by means of interaction with infected poultry β_hp S_h I_p /Nh in a factory. This factor is negative leaving the susceptible individuals and positive enter into the infected class. The factors πh S_h , π_h E_h and πh R_h are the individuals that get vaccinated. Notice that we vaccinate not only the susceptible class, but also the recovered and exposed individuals, since they are not aware of their epidemiological status and hospitals or medical personal will vaccine individuals that are not symptomatic, hence infected individuals will not get vaccinated. The importance of these distinctions is that we can get a rough estimate of the wasted vaccines (i.e. vaccines given to individuals that had already immunity in the recovered class) or will get sick independently since they where incubating the virus already. The factor σh E_h are the individuals that where incubating and now leave the exposed class because they are symptomatic. The factor αh Ih are the individuals that recover from the disease, leaving the class Ih getting into the recovered class with a positive sign. Humans that died from the disease are indicating with the factor -δ_h I_h .

To study the epidemic further, we compute the basic reproductive number R0 , that is the number of people a sick individual infects when inserted into a fully susceptible population. In this case, foran epidemic to occur, we need an R0 > 1. In order to compute R₀ , we will be implementing the second-generation operator developed by Van Driessche and Watmough, in 2002 [15]. The R0 for this model is defined as: R_0p = [β_p / (α_p +δ_p )], R_0h = [β_h / σ_h ] + [β_h / (α_h +δ_h )] that represent respectively the R0 for the poultry epidemic by itself and the R0 for the human epidemic by itself (see Appendix 1 for the computation process).

Results: Simulations of the disease Avian Influenza Type A-H5N1

The parameter values used in the simulations are as described in Table 2, where we use an incubation period for humans of 9 days and change the percentage of vaccinated individuals for the simulations shown and the starting day of the vaccination campaign for Figures 3 and 4.

Figure 3: Vaccinating 10% of the population we vary the initial day of vaccination from 0 (blue), 10 (green), 30 (red) to 50 (turquoise) days. As we delay the vaccination, the number of deaths in the human population changes where on day 30 increase dramatically. The infected population also increases as the vaccination start later changing the duration and the peak of the epidemic.

Figure 4: Vaccinating 30% of the population, we vary the initial day of vaccination from 0 (blue), 10 (green), 30 (red) to 50 (turquoise) days. As we delay the vaccination the number of deaths in the human population changes where on day 30 increase dramatically. The infected population also increases as the vaccination start later changing the duration and the peak of the epidemic. We notice a difference in the epidemic size for starting days 0 and 10 in comparison with Figure 3.

We focus on the behavior of the epidemic entirely by human interaction and the effect of either a single individual or a population. On the results of the simulations when Ro >1, exists the possibility of an epidemic, we focused on the total population of Cayey that is 48,119 people in 2010 according to the census. The mortality rate in infected individuals is greater than 60% according to the World Health Organization (WHO) [16].

The parameter values used in the simulations from Figures 3 and 5 are as described in the Table 2 where we used an incubation period of 9 days for humans and changed the percentage of individuals vaccinated for the simulations shown and the starting day of the vaccination campaign for Figures 3 and 4. In Figure 3, we delay the vaccination from 0 to 30 days after the start of the epidemic and observed that the number of deaths in the human population increases dramatically if vaccination started 30 days after the epidemic. The infected population also increases as vaccines are delayed. In Figure 4, as we delay the vaccination the number of deaths in the human population on day 30 changes where it increased dramatically. We notice a difference in the epidemic size for starting days 0 and 10 in comparison with Figure 3. In Figure 5, as we increased the percentage of people vaccinated the epidemic morbidity and mortality changes as well, it has been shown computationally that when we vaccinate at least 30% of the population we can eradicate the epidemic.

Figure 5: We start the vaccination at the same time the epidemic starts. We change the vaccinated percentage from 10% to 70%. As we increase the vaccination percentage, the epidemic decreases as well as the number of deaths but when we vaccinate at least 30% of the population we can eradicate the epidemic.

SIR five city meta population model

We implemented the modification of a simpler mobility model for five cities in Puerto Rico on a SIR-type model, to investigate the potential outbreak of this disease in the areas of: San Juan, Caguas, Cayey, Ponce and Humacao. The model classifies the population into a unique epidemiological class or compartment per city as: susceptible subjects (S_i ), infected subjects (I_i ), recovered subjects (R_i ) and vaccinated individuals (V_i ). Ni is the total population; where i∈{1,2,3,4,5} indicate the city where the human subjects are located.In order to present our mathematical model let’s define the schematic of the model as shown in Figure 6, the network or connectivity of the cities in Figure 7 and the parameters as shown in Table 3. For the networking approach we assume that the disease mutates such that it is transmitted from human to human and we only focus on the human interaction and the population traveling from one city to another.

Figure 6: Schematic of the mathematical model for avian flu Type A-H5N1.

Figure 7: Schematic of the mobility network.

Table 3: Parameters values of the meta population SIR model for the five cities.

Following the flow of the schematic and descriptions above, we are ready to write down our SIRV system of ordinary differential equations (where the summations Σj pij run for j = 1…5) as follows:

dS_i /dt = µN_i – µSi – βS_i I_i / N_i – πS_i +∑_j p_ij S_j – ∑_j p_ij S_i

dI_i /dt = βS_i I_i / N_i– (µ + α + )I_i +∑_j p_ij I_j– ∑_j p_ij I_i

dR_i /dt = – (µ + π)Rh + αI_i +∑_j p_ij R_j– ∑_j p_ij R_i

dV_i /dt = π(S_i +R_i ) +∑j p_ij V_j– ∑_j p_ij V_i

Similar to the SEIRV model we presented earlier in our meta population model schematic (Figure 6) each box represents a compartment where subjects are categorized into the epidemiological classes using the parameter rates described in Table 3. As before, an arrow entering a compartment indicates inflow of individuals represented in the ordinary differential equations coupled system with a plus sign. Similarly an arrow leaving a compartment indicates that the sign is negative. Important to mention is that the i subscript represents each city.

SIR epidemiological mathematical model networking approach with respect to human population: Entering the susceptible class S i we have the population birth rate µNi for city i and leaving the death rate of human population – µS_i , the infected individuals βS_i I_i / N_i that enter into the infected class I_i . The factors µIi and µRi are the natural death of the human population for the infected and recovered individuals. The recovery ratio of the human population αIh has a negative sign getting out of the infected class and a plus sign getting into the recovery class. Humans that die from the disease are indicated with the factor δIh and vaccinated individuals are indicated by πi getting out of the class S_i and R_i into the V_i class. We vaccinate all individuals that are not symptomatic, since the individual is not aware of his/her epidemiological class. The CDC recommended in the 2008 influenza pandemic that since people are unaware of what strain of the disease they might have suffered prior they should get vaccinated against the disease [17]. j p_ij ’s are the proportion of the population that travels from city i to city j. Hence, the summation of the pij in the differential equations ∑ p_ij , indicates the individuals traveling from other cities to city i. Note that the epidemiological state of the individuals does not change with the mobility because it is an instantaneous change of location.

The mobility component of the model: In this model we introduced traveling to five cities in Puerto Rico as seen in Figure 1 and the schematic of Figure 7 with the intention to extend our model in the future to the 72 municipalities of the island. It is important to mention that the individuals don’t change the epidemiological state while they travel. Note that the parameter values of instant transportation are: p12 that represents the proportion of the population that travels from San Juan to Caguas, p23 is the proportion of the population that travels from Caguas to Cayey, p34 is the proportion of the population that travels from Cayey to Ponce, and finally p25 is the proportion of the population the travels from Caguas to Humacao. The big cities are San Juan with a population of 395,324 individuals, Caguas with 142,893 individuals and Ponce with 166,327 individuals. The smaller cities are Humacao with 58,466 individuals and Cayey with 48,119 individuals [18].

Results: Simulations for the Metapopulation model of Avian flu Type A-H5N1

The parameter values used in the simulations are as described in Table 3 with variations in the initial conditions or vaccination rate.

Multiple peaks of the same epidemic outbreak can be seen in the total population as shown in Figure 8 and Figure 12, depending on the city of the initial outbreak and the connectivity with the other cities in the network we can observe the different shapes of the epidemic. The important result in our simulations is that with the parameters used from Table 3 avian influenza will remain endemic in the population if there is a mutation that will permit the disease to be transmitted from human to human.

Figure 8: The epidemic starts in the city of San Juan and spreads to Caguas, Humacao, Cayey and finally to Ponce. Notice that the shape of the total epidemic changes depending on the initial conditions, giving the notion of having multiple waves of the epidemic as a whole. As time progresses there is a secondary wave with a pick around day 200 and the epidemic remains endemic.

Figure 9: The epidemic starts in the city of Caguas and spreads to San Juan, Humacao, Cayey and finally to Ponce as time progresses there is a secondary wave with a pick around day 200 and we can see that the epidemic remains endemic.

Figure 10 (A): The epidemic starts in the city of Cayey and spreads to Ponce, Caguas, San Juan, and finally to Humacao as time progresses there is a secondary wave with a pick around day 200 and we can see that the epidemic remains endemic.

Figure 10 (B): We start the vaccination at the same time the epidemic starts. We change the vaccinated percentage from 10% to 70%. As we increase the percentage the epidemic decreases as well as the number of deaths, but computationally it can be shown that when we vaccinate at least 30% of the population we can eradicate the epidemic.

Changing the percentage of vaccination

The parameter values used in the simulations are as described in the Table 3 where we use an infection period for humans of 7 days and changed the percentage of vaccinated individuals for the simulations shown and the starting day of the vaccination. We used a non-democratic plan of vaccination [19] by vaccinating individuals in the two main cities in Puerto Rico Caguas and San Juan. In Figure 10, we started the vaccination at the same time the epidemics start. We change the vaccinated percentage between 0%, 10%, 30% and 70%. As we increase the percentage the epidemic decreases as well as the number of deaths (not shown). We compared the vaccination percentage from 10% to 30% only in San Juan and Caguas, but when we vaccinate at least 30% of the population we reduced considerably the morbidity of the epidemic. In fact, we eradicated the disease as shown in Figure 11 in the cities vaccinated (San Juan and Caguas). We must take into consideration the people who were not vaccinated in the cities of Cayey, Ponce and Humacao so there is an epidemic in those cities.

Figure 11 (A): The epidemic starts in the city of Ponce and spreads to Cayey, Caguas, San Juan, and finally to Humacao as time progresses there is a secondary wave with a pick around day 230 and we can see that the epidemic remains endemic.

Figure 11 (B): We start the vaccination at the same time the epidemics starts only in San Juan and Caguas. We vaccinate 30% of the population in these two cities. The epidemic is eradicated in the cities that are vaccinated, but there is an outbreak in the other three cities that are the ones that contribute to the total epidemic waves shown in Figure 10.

Figure 12: The epidemic starts in the city of Ponce and spreads to Cayey, Caguas, San Juan, and finally to Humacao as time progresses there is a secondary wave with a pick around day 230 and we can see that the epidemic remains endemic.

Conclusion

To address our research questions: What mechanisms are effective to contain a potential Avian A-H5N1 epidemic? We discover that at least a 30% vaccination coverage reduces the morbidity of the epidemic significantly. Which dynamics of disease spreading affect an outbreak among cities? The connectivity of the cities, especially where the epidemic starts shape the overall morbidity of a country’s epidemic. What mechanisms are effective to contain a potential avian flu epidemic among cities? If vaccination is not administered democratically and epidemic cannot be contained in the whole country having a lower morbidity over all but an epidemic on cities that have been neglected.

Future Work

Currently, in our research, we focus on working on avian influenza type A-H5N1 that has evolved over the years, producing more and more serious outbreaks that affect a larger number of birds. The increase in those outbreaks is due to the development of the poultry industry in recent years. The research specifically targets the spread of infectious diseases, which will be held in perspective with mathematical epidemiological models to investigate possible alternatives to mitigate the spread using measures such as treatment, immunization and as for the importance to educate the public or the community about the conduct of the epidemic in humans and how to prevent the spread of the disease. Based on our objectives and specific approaches, we want to study and focus on the behavior of the epidemic entirely by human interaction and the effect on the population. We also want to extend our meta population model on the 72 municipalities of Puerto Rico. Our future plans are to continue this research beyond the preliminary results shown in this article to find scenarios where the disease might be eradicated.

Acknowledgment

Funding for this publication was made possible by the Institute of Interdisciplinary Research at the University of Puerto Rico at Cayey and the Building Research Infrastructure and Capacity (BRIC) program from the National Institute on Minority Health and Health Disparities (P20 MD006144) and also by the Research Initiative for Scientific Enhancement (RISE) program (5R25GM059429-17). The views expressed do not necessarily reflect the official policies of the Department of Health and Human Services; nor does mention by trade names, commercial practices, or organizations imply endorsement by the US Government. We would like to thank students and professors at UPR Cayey for their insights on this work.

References

1. Ramis A, van Amerongen G, van de Bildt M, Leijten L, Vanderstichel R, Osterhaus A, et al. Experimental infection of highly pathogenic avian influenza virus H5N1 in black-headed gulls (Chroicocephalus ridibundus). Veterinary research. 2014; 45: 84.

2. García-García J, Ramos C. Influenza, an existing public health problem. Salud Pública de México. 2006; 48: 244-267.

3. Claas EC, Osterhaus AD, van Beek R, De Jong JC, Rimmelzwaan GF, Senne DA, et al. Human influenza A H5N1 virus related to a highly pathogenic avian influenza virus. The Lancet. 1998; 351: 472-477.

4. Montalvo-Corral M, Reséndiz M, Santos-López G, Vallejo-Ruiz V, Reyes-Leyva J, Hernández J. Standardization of a molecular detection method of highly pathogenic avian influenza virus (H5N1). Acta Bioquímica Clínica Latinoamericana. 2009; 43: 49-52.

5. Ghosh A, Nandy A, Nandy P. Computational analysis and determination of a highly conserved surface exposed segment in H5N1 avian flu and H1N1 swine flu neuraminidase. BMC Structural Biology. 2010; 10: 6.

6. Kane MJ, Price N, Scotch M, Rabinowitz P. Comparison of ARIMA and Random Forest time series models for prediction of avian influenza H5N1 outbreaks. BMC Bioinformatics. 2014; 15: 276.

7. Beigel JH, Farrar J, Han AM, Hayden FG, Hyer R, De Jong MD, et al. Avian influenza A (H5N1) infection in humans. New England Journal of Medicine. 2005; 353: 1374-1385.

8. Kermack W, McKendrick A. Contributions to the mathematical theory of epidemics. Proceedings of the Royal Society of London. 1927; 115: 700-721.

9. Montesinos-López OA, Hernández-Suárez CM. x Modelos matemáticos para enfermedades infecciosas. Salud Pública de México. 2014; 49: 218-226.

10. Boyev BV, Reyes TG, Gómez AG. Modelos Analíticos de Epidemias con Fines de Pronóstico II. Foro-Red-Mat: Revista electrónica de contenido matemático. 2005; 16: 1.

11. Nyuk Sian Chong, Jean Michel Tchuenche, Robert J. Smith? A mathematical model of avian influenza with half-saturated incidence. Theory in Biosciences. 2014; 133: 23-38.

12. Wang L, Li X. Spatial epidemiology of networked metapopulation: An overview. Chinese Science Bulletin. 2014; 59: 3511-3522.

13. Arino J, Van den Driessche P. A multi-city epidemic model. Mathematical Population Studies. 2003; 10: 175-193.

14. Arino J. Diseases in metapopulations. Modeling and dynamics of infectious diseases. 2009; 11: 65-123.

15. Van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences. 2002; 180: 29-48.

16. World Health Organization.

17. Centers For Disease Control and Prevention: Vaccine against 2009 H1N1 influenza virus.

18. Herrera-Valdez MA, Cruz-Aponte M, Castillo-Chavez C. Multiple outbreaks for the same pandemic: Local transportation and social distancing explain the different “waves” of A- H1N1-pdm cases observed in México during 2009. MBE. 2011; 8: 21-48.

19. Cruz-Aponte M, McKiernan EC, Herrera-Valdez MA. Mitigating effects of vaccination on influenza outbreaks given constraints in stockpile size and daily administration capacity. BMC infectious diseases. 2011; 11: 207.