Abstract: In this study, a system of SIR models depended on effective contact rate of the disease and potential size of the community is developed for simulating patterns of pandemic by influenza A during the Haj. Indeed, we consider three different SIR models such that each of them is calibrated to simulate a separate part of the ceremony and then we link these models together to obtain system simulating whole the Haj ceremony. Considering different numerical values of effective contact rate and potential size of the community, we estimate the portion of susceptible hosts in the community and also the portion of total infection during the ceremony and portion of infective individuals at any time of the ceremony. Our goal is to suggest doable optimal situations keeping the pandemic under the control and the number of infected people as low as possible.
INTRODUCTION
For who are interested in study of infectious diseases, the influenza is of great importance. Since, influenza virus continues always evolving, so there exists a menace of pandemic by mutant influenza virus and it is still an important problem for the public health systems. In 1918, the great pandemic of influenza on modern history occurred. About 40000000 humans died of this pandemic all over the world. After 1918, a pandemic of influenza occurred twice in 1957 and 1968 and today again, the World Health Organization (WHO) has warned of a substantial risk of pandemic of influenza A throughout the globe. Apart from the mortality associated with complications due to influenza, the costs of workdays lost by the infected people, vaccines, antiviral drugs and so on must be taken in account among the consequences of probably pandemic by influenza. Meanwhile, Islamic countries with near to 2 billions populations, mostly dwelled in Asia, are threaten by pandemic of influenza during the Haj ceremony when almost 2 millions of Moslems gather in Mecca every year. Islamic republic of Iran, as country with more than 70 millions populations joints this ceremony every year with tens of thousands of Iranian people. In 2009, because of the threats existing for pandemic by influenza A, the government put some limitation as they did not let the people with more that 63 years old and children to joint the ceremony. Also they stated a restriction of number of people jointing the ceremony. Similar policies took by a few other Islamic countries in this regards. These facts convinced the author that any effort on developing models for illustrating, estimation and simulation patterns of infection within the Haj ceremony may be helpful.
In this regards, mathematical models are considered as easily accessible tools for simulating patterns of infection by influenza. Especially, when these models are supported by exact statistical data, the consequent results are applicably reliable. Many mathematical models have been proposed in the study to describe the inter-pandemic ecology of influenza A in humans. A typical approach is modeling of the interactions between individuals who are (or have been) infected by one or different viral strains (Castillo-Chavez et al., 1989; Brauer et al., 2008; Brauer and Castillo-Chavez, 2001). An advanced modeling is done by Gupta et al. (1998), where the researchers showed that the simultaneous circulation of several antigenic variants of the same pathogen can give rise to complex dynamics, including sustained oscillations and chaos in disease. This may occur when a genetic mutation for influenza A causes simultaneous existence of two or more various strain. Then the results were illustrated more specifically by Lin et al. (1999), where an epidemiological model consisting of a linear chain of three co-circulating influenza A strains that provide hosts exposed to a given strain with partial immune cross-protection against other strains is analyzed. They could find a sub-model that exhibits sustained oscillations. Since, genetic mutation for influenza A is occurred usually ones within per several years, so the models which simulate mid time patterns, consider pandemic by a single strain. Such simulations have been done by Hethcote (1976), Casagrandi et al. (2006) and Thieme and Yang (2002), where a SIR or SIRC model develops a simple ordinary differential equation for studying patterns of infection in a certain mid-time interval. For an especial case by Andreasen et al. (1997), a SIRC model is considered to simulate the epidemiological consequences of the drift mechanism for influenza A viruses. Their model could predict a rich variety of complex temporal patterns that are realistic for influenza A. One year later, a mathematical model is proposed to interpret the spread of avian influenza from the bird world to the human world (Iwami et al., 2007). They gave their results on the spread of avian influenza and mutant avian influenza and measures to control the spread of avian influenza. Since, the complete immunity is usually lost within a more that several weeks, so the short time patterns are usually simulated by simple SIR models (Brauer and Castillo-Chavez, 2001; Hethcote, 2000; Anderson and May, 1991).
In this study, we consider a set of SIR models which simulate patterns of infection by influenza A during the stages of the Haj ceremony. These models depend on several parameters, including the effective contact rate and the potential size of the community. The parameters are tried to be well calibrated for the Iranian Haj ceremony; however; these systems can be implemented simply for other communities or diseases by selecting the values of the parameters appropriately.
MATERIALS AND METHODS
Introducing the SIR System and Parameters
In a short period of time, infectious diseases follow epidemiological patterns
of infection which is called in brief SIR models. That is : individuals in susceptible
subgroup are infected by the disease and expend a period of time as member of
infective subgroup and finally a portion of them are recovered from the disease
and the rest of them died. In a short time, recovered individuals are assumed
to have complete immunity against the disease. Figure 1 shows
the corresponding model.
Since, the Haj ceremony is usually accomplished within a period of 30 days, so SIR models are appropriate for modeling of pandemic by influenza A. Therefore, we divide the Haj community at time t into three compartments: susceptible individuals S(t) of who have no specific immunity against influenza A; the fraction I(t) of those individuals that are infected and infective; the fraction R(t) of those individuals who are recovered from the disease and gain complete immunity (at least for a short period of time which we consider of). Also we assume that the community of people in Haj has no natural mortality or birth. On the other hand, the Haj community is constituted within three steps:
| Step 1: | With a period of 13 days in which the people receive with a constant rate to the place of the ceremony |
| Step 2: | With a period of 4 days, in which the major traditions of the ceremony is accomplished |
| Step 3: | With a period of 13 days, in which the people return with a constant rate and fixed schedule from the ceremony to the country |
| Fig. 1: | A simple figuration of subgroups interacting in a SIR system |
If N0 is the number of the people contributing in the ceremony then we can say that N0 is the potential size of the Haj community. Thus, we can define the corresponding SIR models for step 1-3, respectively by the following set of three ordinary differential equations.
| (1) |
| (2) |
| (3) |
Equation 1 explains that the susceptible hosts arrive to the ceremony place by constant rate N0/13 and they are infected at rate αS1I1; then, the infective hosts are recovered at rate εI1, die at rate μI1nd are generated at rate αS1I1; also the recovered hosts are generated at rate εI1. The Eq. 2 and 3 use the same structure depending on the stage of the ceremony that they simulate. These Eq. 1-3 are respectively subjected to the following initial conditions.
| (4) |
The condition I1(1) = 1 is considered because of the fact that any pandemic begins by existence at least one infective individual and we assume that the first infection appears in the first day of the ceremony, also in the first day of the ceremony, the total size of the community is N0/13. The parameter μ stands for the rate of mortality due to the infection and the parameter ε stands for the rate of recovering from the disease. Therefore, ε+μ shows the rate of individuals who leave the infective subgroup. Therefore, if the average number of days which a specific infective individual expends to be recovered or die is l, then μ+ε = 1/l day-1 (Brauer and Castillo-Chavez, 2001; Hethcote, 2000; Andreasen et al., 1997). Thus, ε = (1/l)–μ. In this study, according to the statistical data (Castillo-Chavez et al., 1989; Lin et al., 1999; Andreasen et al., 1997), we assume that the average time expended by the individuals in the infective subgroup is 4 days i.e., l = 4; also the mortality rate is considered 1% i.e., μ = 0.01/l = 0.0025 and therefore ε = 0.99/l = 0.2475.
Computation of the parameter α is more difficult. This parameter is called effective contact rate and it must be measured as effective contacts per unit time. So, it defers depending on the size and social specifications of the given population. This may be expressed as the total contact rate (the total number of contacts, effective or not, per unit time, denoted γ), multiplied by the risk of infection, given contact between an infectious and a susceptible individual. This risk is called the transmission risk and is denoted by p; (Brauer and Castillo-Chavez, 2001). Thus: α = γp. On the other hand p = x/N, where N denotes the total size of the community and x denotes the number of secondary infected cases in the infective subgroup. Thus: α = γx/N, which means that α is of scale of (1/N) day-1 and so it is usually considered a small parameter. The total size of the community at any time t is obtained from N = S+I+R. Thus, as time goes farther, the total size of the community for systems 1-3 holds, respectively for Eq. 5 below.
| (5) |
In order to bring the systems Eq. 1-3 into dimensionless standard form, first we divide variables N, S, I, R to N0. This makes the variables dimensionless. Then we rescale the time with τ = 2(μ+ε)t/15 = t/30. Thus τ is dimensionless. Therefore the dimensionless SIR systems are obtained, respectively as below:
| (6) |
| (7) |
| (8) |
These Eq. 6-8 are subjected, respectively to the following new initial conditions.
| (9) |
These three sets of ordinary differential equations subjected to the corresponding initial conditions are considered as the standard dimensionless SIR models for simulating the patterns of infection by influenza A. Since, the total size of the community N(t) is always less than N0, so in systems 6-8, we have S(t)+I(t)+R(t)≤1. Furthermore, according to the way we rescaled parameter t, when t varies from 0 to 30 days then τ varies from 0 to 1. This means that for systems Eq. 6-8, we have respectively 0≤τ≤13/30, 13/30≤τ≤17/30, 17/30≤τ≤1. Now we will study the patterns of infection predicted by these models and try to find critical values for effective contact rate α.
RESULTS
Estimating Patterns of Infection
In this study of a pandemic by the infectious disease, there are three matters
of importance: First of all, the patterns by which the infection develops in
the community; the second one is the total number (or total portion) of all
people in the community who get infection at a certain time. This is especially
important because it helps us to know how much of vaccines, antiviral drugs
or other necessary medical cares we should be able to supply during this certain
time. Finally, the third one is the number (or portion) of the people in the
community who are infected at time t. This helps us to know how much of medical
cares is needed at any time t. In this section we try to indicate these three
matters of importance. We begin with the first one.
It is obvious that here the goal is to keep the number of people who get infection during the Haj ceremony as low as possible. Therefore, the patters which describe the number (or portion) of the individuals in each of susceptible, infective or recovered subgroup are interesting to us; especially, the portion of susceptible hosts in compare with the portion of infective and recovered hosts is an appropriate variable which describes the pattern by which infection develops during the Haj ceremony. Thus, we define two new variables A(τ) = S(τ)/N(τ) and B(τ) = (I(τ)+R(τ)/N(τ)). By this definition, the variable A(τ) shows the portion of non-infected people in the community and the variable B(τ) shows the portion of the people who belong to infective or recovered subgroups in the community. These variables are nonnegative and hold for the equation.
Thus, the graphs of these two new variables can describe development of infection during the Haj ceremony. In the Fig. 2, the numerical solutions of systems 6-7, computed for different values of the parameters α and N0, are applied for plotting the graph of A(τ) and B(τ). In the Fig. 2a-r, the green and black curves are, respectively stand for A(τ) and B(τ), where the horizontal axis stands for τ and the vertical axis stands for the values of A(τ) and B(τ); also the vertical dash lines separate three distinct parts of the Haj ceremony. The values of the parameters α and N0 for each of the Fig. 2 are also indicated in its below. In the Fig. 2a-f, the parameter α is constantly 5x10-6 while the parameter N0 varies from 5x104 to 105. This means that the average number of secondary infection by each infective individual is 5 person for per 1 million susceptibles. As it can be seen, the portion of non-infected individuals decays by growth of size of the community. However, as the Fig. 2 show the portion of infected people is very low. This means that, as long as the effective contact rate α is lower that 5x10-6, the pandemic will be under control.
In the Fig. 2g-l, the effective contact rate α is constantly 1x10-5 and the parameter N0 varies from 5x104 to 105. As it can be seen, the portion of infected individuals increases fastly by growth of size of the community. Especially, for N0 = 105 the green curve intersects decreasingly the black curve which means that the portion of infected people will be more that the portion of the non-infected people after the intersection point. In this case, the value of τ which corresponds to the intersection point can be considered as the time of breaking out of infection.
| Fig. 2: | (a-r) The graph of variables A(τ) and B(τ) for various range of the parameters α, N0 according to the text. Here the green and the black curves are stand by respectively for A(τ) and B(τ) and the dashed lines indicate three separate parts of the Hajj ceremony |
Finally, in the Fig. 2m-r, the effective contact rate is constantly 1.5x10-6 and the parameter N0 varies from 5x104 to 105. As it can be seen, the pandemic occurs of N0≥7x104. Figure 2a-r show that the time of breaking out decays by growth of size of the community. Moreover, for N0>9x104, the breaking out time may decay to the last days of the second part of the ceremony (Fig. 2q , r) which is the most important part. This means that increasing the value of effective contact rate to 1.5x10-5 and N0 to more that 9x104 disturbs the ceremony.
Thus, the model estimates that, in order to have the pandemic under the control, the effective contact rate α must not exceed than 5x10-6 (Fig. 2a-f) and if it exceeds up to 1x10-5 then the total size of the community N0 is better to be less than 9x104 (Fig. 2g-k). Furthermore, if α exceeds to 1.5x10-5 then, in order to control the pandemic, N0 must be less that 6x104 (Fig. 2m, n).
The second important matter is the portion of total infected individuals during the ceremony. This is the portion of all individuals in the community who are infected or recovered at day 30 of the ceremony plus to the number of individuals who died due to the infection. To obtain this portion, firstly we have to compute the number of all infected people up to t-th day of the ceremony, then total number of all infected individuals during the Haj ceremony is obtained for t = 30 and the portion of total infection is obtained by rescaling the variables and putting τ = 1.
From the Eq. 5 and the initial conditions Eq. 4, we have,
In each of the three equation above,
shows the number of individuals in the infective and recovered subgroups who have left the ceremony (and returned to the country) up the t-th day of the third part of the ceremony. By the same way,
shows the number of individuals in the susceptible subgroup who have left the ceremony (and returned to the country) up the t-th day of the third part of the ceremony. This yields that the total number of infected individuals at each part of the ceremony, which we denote it by Γ, is obtained by the following equations:
Thus, the dimension less variable
Thus, the portion of total infection which we denote it by Ω is equal
to
| Fig. 3: | (a-f) The total portion of infected individuals for various range of N0 according to the text. Here, the vertical axis stands for Ω and the horizontal axis stands for α |
As it can be seen through Fig. 3a-d, for the effective contact rate α≤10-6and N0≤8x104 the portion of all infected individuals is less that 5% of the community however, this portion may suddenly rise up to 20% of the community for N0≥8x104 (Fig. 3e) and even to 40% of the community for N0≥9x104(Fig. 3f). This means that the critical value for the potential community size N0 is approximately 8x104.
| Fig. 4: | (a-r) The graphs of portion of infective individuals according to the text. Here, the vertical axis stands for total portion of infected individuals and the horizontal axis stands for τ |
This is because of the fact that any exceeding from this value may cause a sudden breaking out of the infection. One more important fact which can be observed is the existence of a relatively sudden rise up in the portion of total infection by exceeding the values of α from 10-5. This means that the critical value for effective contact rate is 10-5. Therefore, in order to keep the portion of total infection lower than 25%, the secondary infections must be keep lower that one person in per 100000 susceptible individuals.
At the end of the third important matter i.e., the number of the people how are in infective subgroup at t-day of the ceremony. This is directly obtained from the values of the variable I(t) in Eq. 1-3. Similarly, the dimensionless variable portion of the infective individuals can be obtained through the values of I(τ) in Eq. 6-8. The Fig. 4a-r show this portion for various values of parameters α and N0, where the horizontal axis stands for τ and the vertical axis stands for the portion of individuals in infective subgroup.
As it can be seen, probably except for the Fig. 4a, we have an increase in portion of infective individuals within the second part of the ceremony and this increase may continue even up to the last day of the ceremony (Fig. 4i-k and m, n).
The peak point in each of the Fig. 4 indicates the time in which the portion of the infective hosts is in its maximal portion. Since, accomplishing the traditions of the ceremony is hard for the individuals how is being infected by the disease, so it is better to keep the peak of the Fig. 4 far from the second part of the ceremony, especially for those figures which are according to a great portion of total infection in Fig. 3. For example, Fig. 3d-f with α = 1.5x10-4 and N0≥8x104 and Fig. 4p-r. Also, according to the Fig. 3, 4, we must expect a flash rise up in portion of infective individuals right after the second part of the ceremony. For α≤10-4, this rise up is less that 2 percents of the community (Fig. 4a-l), but if α increases it may even reach to 50% of the community (Fig. 4q, r).
DISCUSSION
We developed a system of SIR models for simulating the probably pandemic by influenza A during the Haj ceremony. We developed SIR models to three sets of ordinary differential Eq. 1-3 and in order to avoid complexity due to scales, we rescaled these equations to dimensionless differential Eq. 6-8, which are subject to appropriate initial conditions (9). Then we considered three matter of importance: the portion of susceptible hosts in compare with the portion of infective and recovered hosts in the community; the total portion of infected individuals and the portion of infective hosts at any time.
First in the Fig. 2, we saw that in order to have the portion of infective and recovered hosts low in compare with the susceptible hosts, the effective contact rate α must be less that 5x10-6 (Fig. 2a-f). We also saw that of the effective contact rate exceeds up to 10-5 then it is necessary to keep the number of the people how join the ceremony N0 less than 8x104 (Fig. 2g-j). Also, we saw that for α≈1.5x10-5 there exists a breaking out of the disease which may begin even from the last days of the second part of the ceremony (Fig. 2q, r).
In the Fig. 3, we plotted the total portion of infection and discussed about different cases which occurs. Especially we indicated the critical values for α and N0; for example, we saw that if the value of α exceeds from 10-5, then we expect a relatively sudden rise up in the total portion of infection.
Finally, in the Fig. 4, we plotted and discussed about the different cases of portion of infective hosts at time t of the ceremony. We distinguished an increase in portion of infective individuals within the second part of the ceremony and we saw that for α≥10-5 this increase may continue even up to the last day of the ceremony (Fig. 4i-k and m, n). Also, we found the cases for which the peak of infection occurs in earlier days of the third part of the ceremony (Fig. 4q, r).
From all discussion we conclude that the these SIR models predicts that the optimal situation for the Haj ceremony which is obtainable in real world can be occurred for N0 at most up to 8H104 and medical cares which keep the effective contact rate α less than 10-5. If it happens then, according to Fig. 2-4, we expect that the portion of infective or recovered individuals and the total portion of infection hold less that 10%. Also in this case the peak of infective individuals is less that 1.5% and it occurs in last days of the ceremony which is not much important because at that time much part of the community has left the ceremony and returned to the country.