Nonlinear Pulse Vaccination in an SIR Epidemic Model with Resource Limitation

and Applied Analysis 3 Table 1: Parameter definitions, values and source for model (4). Parameter Definition Value Source b Rate of natural birth or death rate of the population (year) 0.02 [33] β Probability of transmission per contact (year) 1800 [33] γ Natural recovery rate of the infective individuals (year) 100 [33] pmax Themaximal vaccination proportion (year ) [0, 1) [8, 9] θ The half saturation constant [0, 1] Assumption T Period of pulse vaccination (year) [0, 20] [8, 9] The formula for p(t) plays a key role in investigating the effect of the limited capacity for treatment on the spread of infectious disease. To do this we assume, without loss of generality, that the PVS is implemented in a developing country but one in which vaccines are relatively abundant. Note that there are often major shortages of medical personnel in rural areas and in certain specialities such as surgery, paediatrics, and obstetrics [14]. Let vaccination pulses occur every T time units. When a pulse day occurs, the medical personnel have to go out into the rural population and vaccinate as many susceptibles as they can find on that particular day. If there are very few susceptibles, then it will be hard to find them, so the vaccination coverage will be low. The more susceptibles there are, the easier it is to find them, so the coverage per pulse will be greater. But there is a limit to how fast the medical team can find and handle each susceptible, and this causes the vaccination rate p(t) to saturate, that is, a saturation phenomenon of limited medical resources. By employing the Holling type II functional response for a predator-prey model, we can define p(t) as follows: p (t) = pmaxS (t) S (t) + θ , 0 ≤ pmax < 1, (3) where pmax denotes the maximal vaccination proportion and θ is the half saturation constant. Note that the vaccination coverage is less than 100% and in practice a vaccine efficacy is also less than 100%, which lead to the parameter pmax < 1. Furthermore generalized formulae such as the Hill function [29] can be employed to characterize the saturation phenomenon of the limited medical resources. Taking into account the resource limitation and saturation effects, for convenience and simplification, then model (2) can be rewritten as dS (t) dt = b − βS (t) I (t) − bS (t) , dI (t) dt = βS (t) I (t) − γI (t) − bI (t) , t ̸ = nT, n ∈N, S (t + ) = (1 − pmax) S 2 (t) + θS (t)


Introduction
Epidemiology is the study of the spread of diseases with the objective of tracing factors that are responsible for or contribute to their occurrence and serves as the foundation and logic of interventions made in the interest of public health and preventive medicine. Mathematical models describing the population dynamics of infectious diseases have played an important role in better understanding epidemiological patterns and disease control for a long time. Various epidemic models have been proposed and explored extensively and considerable progress has been achieved in the studies of disease control and prevention (see [1][2][3] and the references therein).
Outbreaks of infectious diseases have not only caused the loss of billions of lives but have often also rapidly damaged social economic systems, bringing about much human misery. Consequently, the focus of our research has been on how to prevent and cure infectious diseases effectively. It is well known that one of the most important concerns in the analysis of epidemic logical models is the efficacy of vaccination programmes. This subject gained prominence as a result of highly successful application of vaccinations for the worldwide eradication of smallpox and the restriction of diseases such as poliomyelitis, measles, tetanus, diphtheria, pertussis, and tuberculosis. Vaccination is one of the most cost effective of health investments as it can prevent or ameliorate morbidity due to infections. The effectiveness of vaccination has been widely studied and verified for the influenza vaccine [4], the HPV vaccine [5], the chicken pox vaccine [6], and others.
In particular, Agur et al. [7] first proposed a pulse vaccination strategy (PVS), which consists of periodical repetitions of impulsive vaccinations of all the age cohorts in a population, which has been confirmed as an important and effective strategy for the elimination of infectious diseases. At each vaccination time a constant fraction of susceptible people is vaccinated. This kind of vaccination is called impulsive since all the vaccine doses are applied within a very short time, with respect to the dynamics of the target disease. PVS allows the eradication of a disease with some practical advantages, as discussed in [7][8][9]. Recently, epidemic models with pulse vaccination have been the subject of intense theoretical analysis. For example, the dynamical behaviors of various epidemic models with PVS are studied in [8,[10][11][12][13]. The theoretical results show that PVS can be distinguished from the conventional strategies in leading to disease eradication at relatively low values of vaccination (see [7]).
Traditional epidemic models with PVS of population dynamics have assumed that the pulse vaccination proportion is constant, which implies that medical resources such as drugs, vaccines, hospital beds, and isolation units are sufficient for the infectious disease in question. However, in reality, every community or country has an appropriate or limited capacity for treatment, especially for emerging infectious diseases, and so understanding resource limitation is critical to effective management.
During the last decade, the outbreaks of SARS in 2003 and avian influenza among humans (H5N1 in 2003, H1N1 in 2009, and H7N9 in 2013) emphasized the need to enhance the capacity to fight emerging infectious diseases, which remains a challenge for public health. Several different vaccines have been developed specifically targeting avian influenza among humans, but limited production has exposed weaknesses [14][15][16]. A lack of sufficient vaccine at the outset of a new pandemic is to be expected, given the crucial lag between the timing of the emergence of a new strain and when a vaccine has been developed and is ready for distribution. Although most vaccination programmes aim to vaccinate as many susceptible individuals as possible, this is not always the case for emerging infectious diseases due to limited vaccine availability, particularly in developing countries.
However, resource limitation is difficult to quantify partly because it is a dynamic process (see [17]). Recently, continuous SIR models concerning limited medical resources have been developed and investigated (see [18,19]). Chow et al. [18] explored the dynamics of an SIR epidemic model to understand how limited medical resources and their supply efficiency affect the transmission of infectious diseases. Zhou and Fan [19] studied a multigroup SIR epidemiological model to explore the effects of limited medical resources and group-targeted vaccination strategies on disease control and prevention. To the best of our knowledge, no work has been done for the effects of resource limitation on an SIR model with PVS. In order to investigate the effect of limited vaccine availability on the spread of infectious disease, a saturation phenomenon of limited medical resources is considered. That is, we will study the dynamic behavior of an SIR epidemic model with nonlinear pulse vaccination. Such mathematical models are suitable for simulating processes with short duration perturbations during their development (see [20]).
In this paper, the dynamical behavior of an SIR epidemic model with nonlinear pulse vaccination is proposed and analyzed. The main purpose is to address and understand how a limited vaccine resource affects the transmission and control of infectious diseases. The rest of this paper is organized as follows. In the next section, an SIR epidemic model with nonlinear PVS is introduced. In Section 3, we investigate the existence and stability of the disease-free periodic solution by using the method of differential inequality, qualitative analysis, a discrete dynamical system determined by the stroboscopic map, and comparison theorem. Meanwhile, Latin Hypercube Sampling (LHS)/Partial Rank Correlation Coefficients (PRCCs) uncertainty and sensitivity analysis techniques are employed to investigate the key control parameters which are most significantly related to threshold values (see [21][22][23]). Section 4 focuses on finding the sufficient condition under which the endemic periodic solution exists by using the bifurcation theorem. The paper ends with some interesting biological conclusions and numerical bifurcation analyses, which complement the theoretical findings.

Model Formulation
In the classical epidemiological model (see [24][25][26]), the total population is composed of three groups of individuals: susceptible ( ), infective ( ), and recovered ( ). Let ( ) denote the number of members of a population susceptible to the disease at time ; ( ) represents the number of infective members; and ( ) denotes the number of members who have been removed from the population. This leads to the following SIR model: (1) In model (1), we assume that the three classes of subpopulations have the same constant birth and death rates (i.e., the total population has a constant size), which can be normalized to unity; that is, ( ) = ( ) + ( ) + ( ) = 1, and the infectious individuals recover at a rate > 0, so that 1/ is the mean infectious period. Susceptibles become infected at a rate , where > 0 is the transmission rate. In practice, the equation for d ( )/d is redundant because ( ) can be obtained from ( ) = 1. A detailed description of model (1) and its dynamics may be found in [9,27,28]. Let > 0 be the time between two consecutive pulse vaccinations and let ( ) (0 ≤ ( ) < 1) be the fraction of susceptible subjects who are inoculated with vaccine, which depends on the number of susceptibles in the population at time . This yields the following model with PVS [8,9]:  [0, 20] [8,9] The formula for ( ) plays a key role in investigating the effect of the limited capacity for treatment on the spread of infectious disease. To do this we assume, without loss of generality, that the PVS is implemented in a developing country but one in which vaccines are relatively abundant. Note that there are often major shortages of medical personnel in rural areas and in certain specialities such as surgery, paediatrics, and obstetrics [14]. Let vaccination pulses occur every time units. When a pulse day occurs, the medical personnel have to go out into the rural population and vaccinate as many susceptibles as they can find on that particular day. If there are very few susceptibles, then it will be hard to find them, so the vaccination coverage will be low. The more susceptibles there are, the easier it is to find them, so the coverage per pulse will be greater. But there is a limit to how fast the medical team can find and handle each susceptible, and this causes the vaccination rate ( ) to saturate, that is, a saturation phenomenon of limited medical resources. By employing the Holling type II functional response for a predator-prey model, we can define ( ) as follows: where max denotes the maximal vaccination proportion and is the half saturation constant. Note that the vaccination coverage is less than 100% and in practice a vaccine efficacy is also less than 100%, which lead to the parameter max < 1. Furthermore generalized formulae such as the Hill function [29] can be employed to characterize the saturation phenomenon of the limited medical resources.
Taking into account the resource limitation and saturation effects, for convenience and simplification, then model (2) can be rewritten as where the parameters , , , and max , are defined in models (1) and (3), respectively. Note that model (4) is reduced to the classical SIR model with constant pulse vaccination when = 0, which was investigated by Agur et al. [7], Stone et al. [8], and Shulgin et al. [9]. The parameter definitions of model (4) are summarized in Table 1.

Existence of the Disease-Free Periodic Solution
We first demonstrate the existence of a disease-free periodic solution of model (4), in which infectious individuals are entirely absent from the population permanently; that is, ( ) = 0 for a sufficiently long time. Under some threshold conditions, we show below that the susceptible population oscillates with period in synchronization with the periodic pulse vaccination and the infectious class will die out eventually. Thus, in this section, we focus on determining the global attraction of this disease-free periodic solution. To address this, we first consider the following submodel over the time interval < ≤ ( + 1) : Integrating the first equation of model (5) between pulses yields for < ≤ ( + 1) . It follows from the second equation of model (5) that Substitution of (6) into (7) gives Abstract and Applied Analysis Denote ( + ) = then the previous equation can be rewritten as the following difference equation: which is the so-called stroboscopic map of model (5), and it describes the relations of the number of susceptibles in the population between any two successive pulse vaccinations. Consequently, the existence of the positive steady state of model (9) implies the existence of a positive periodic solution of model (5). Therefore, we first discuss the conditions which guarantee the existence of a positive steady state of (9). Taking the derivative of ( ) with respect to yields and it is easy to see that 0 < ( ) < 1 holds true. For convenience, we denote the positive fixed point of the stroboscopic map (9) bỹ, which satisfies the following equation: where Obviously, (11) has a unique positive root; that is, which is stable due to 0 < (̃) < 1.
By using Lemma 1, we have thus shown that the unique fixed point̃of (9) is globally stable.
In the following, we will present the sufficient condition for the global attraction of the disease-free periodic solution ( * ( ), 0) of model (4).
The proof of Theorem 2 is given in Appendix A. Note that although the threshold value 1 0 depends on all parameters of model (4), the most interesting parameters here are pulse vaccination period , maximum proportion of pulse vaccination max , and the parameter related to resource limitation, that is, . In order to address how those factors ( and ) affect the threshold value 1 0 , we first note that which indicates the importance of the frequency of vaccination in eradicating an infectious disease. Then letting and vary and fixing all other parameters as those in Figure 1, we see that the threshold value 1 0 is a monotonically increasing function of for fixed , and it is a monotonically increasing function of for fixed period . Moreover, there are some critical values of and such that 1 0 = 1, which indicate that there exists one maximum allowable period of pulse vaccination max such that 1 0 < 1 for all < max . Moreover, the larger the smaller the value that max has. Unfortunately, Abstract and Applied Analysis 5 we cannot get its analytical expression due to the complexity of 1 0 . So we turn to find some approximations for max in the following.
If the PVS is applied frequently enough such that d ( )/d ≤ 0 for all > 0 in model (4), then the number of infectious individuals is a decreasing function of time. It is possible to satisfy this condition if pulsing ensures that ( ) ≤ ≐ ( + )/ for all ; that is, pulse vaccination is applied every time once the ( ) approaches the threshold value (see [7,27] and the references therein).
It follows from the first equation of model (4) that Consider the following comparison equation: which gives Hence and solving the above inequality with respect to yields which indicates that Thus, if the nonlinear pulse vaccination strategy is applied periodically with a period less than max , then the number of infectious in the population will decrease forever and eventually die out. We also note that if = 0, the maximum allowable period max becomes 0 max = (1/ ) ln(1 + (( max )/(1 − ))), and this is the case studied in [7,9] for a constant proportion pulse vaccination strategy. Obviously, max is a monotonically decreasing function with respect to , which gives that max < 0 max . This theoretical result indicates that if we aim to eradicate an infectious disease, then it is necessary to carry out pulse vaccination programs more frequently under resource limitation than when necessary resources are available. When various nonlinear factors including max , are seriously considered, our results show that it is getting more and more difficult to eliminate the endemic diseases. It also confirms that getting vaccinated is by far the most effective action that susceptibles can take to protect themselves and their family from infectious disease. Other options available to health practitioners in response to limited medical resources include the release of prescriptive health education, training vaccinators, and promoting the use of vaccine-delivery patch with dissolving microneedles (see [14][15][16] for details).
A typical solution of the SIR model with a nonlinear pulse vaccination strategy is shown in Figure 2, where we observe how the variable ( ) oscillates in a stable cycle (as shown in Figure 2(a)). In contrast, the infection ( ) rapidly decreases to zero (Figure 2(b)) if we set = 0. That is, without resource limitation the infectious disease is eventually eradicated (Figure 2(b)). However, if there is resource limitation, that is, the parameter is larger than zero, then both susceptible and infected populations oscillate periodically, as shown in Figures 2(c) and 2(d) with = 0.2.
Therefore, if we have enough vaccine and vaccinators to combat the infectious disease, then it is reasonable that the susceptible population vaccinated each time is assumed to be proportional to the number of susceptible individuals (here max ). However, this approximation cannot reflect the real case if we aim to mitigate the emerging infectious disease. Usually, the number of susceptible individuals needing to be vaccinated may exceed the capacity of local medical conditions due to shortages of vaccine and doctors, especially in rural areas in many developing countries, where reaching all of the target population may be difficult. For instance in 2010 in Lesotho, diphtheria, tetanus toxoid and, pertussis (DTP3) immunization coverage among 1-year-olds was only 83% (but an improvement on 78% in 2001, see [31]) and both the measles and tuberculosis (BCG) vaccination programs were hampered by supplies of vaccine running out in some parts of the country and then allowing vaccinations only for children born in hospitals and not for those born elsewhere in the communities. Thus, it is important to design pulse vaccination campaigns carefully. For example, how frequently should the vaccination strategy be applied? And how can the size of the susceptible population be evaluated? Those are key issues for evaluating vaccination coverages and the distribution of vaccines. Therefore, with the help of mathematical models the vaccination strategy can be properly designed.
It is well known that in order to eradicate an infectious disease, the vaccination strategy must be designed such that the threshold value 1 0 is less than 1, which is equivalent to In the following we employ the formula 2 0 to investigate the important factors which affect the threshold values most significantly by using uncertainty and sensitivity analysis (see [21,22,32]). Sensitivity analysis of the most significant parameters including the birth or death rate ( ), the contact rate ( ), the maximum proportion vaccinated ( max ), and the pulse vaccination period ( ) was performed by evaluating the PRCCs for various input parameters against the threshold condition LHS with 3,000 samples, uncertainty and sensitivity analyses for all parameters in model (4) were determined. A uniform distribution function was used and tested for significant PRCCs for all parameters with wide ranges, such as ∼ (0.001, 0.1), ∼ (1000, 2500), ∼ (0. 1,20), and ∼ (0, 1), and the baseline values of all parameters are given in the figure legend of Figure 3. Note that for comparative purposes, we chose the same baseline parameter values as those used in literatures [8,9,33]. Unfortunately, we do not have the realistic ranges of parameter values and we chose those parameter ranges just for illustrative purposes only and then provided a scientific basis for disease control. Therefore, the parameters , , , and are responsible for increasing the values of 2 0 , so decreasing all those four parameters can reduce the threshold value and consequently are beneficial for infectious disease control. As the parameter maximum proportion of pulse vaccination max changes, it is negatively and significantly correlated with large PRCCs. Thus, increasing this parameter will result in significantly reducing the threshold value. The most significant control parameters are , , max , and (as shown in Figures 3(b), 3(c), 3(e), and 3(g)). It is interesting to note that the period, , and the maximum proportion vaccinated, max , are strongly correlated with the threshold value 2 0 , which indicates that carefully designing a pulse vaccination strategy is crucial for infectious disease control.

Existence of Endemic Periodic Solution and Complex Dynamics
An important issue is how the infectious disease develops if the threshold condition 1 0 is larger than one. In particular, what types of dynamic behavior may model (4) have once the period of pulse vaccination exceeds certain threshold levels. Therefore, in this section we first focus on the existence of an endemic periodic solution, that is, to investigate under what conditions will the infected population oscillate periodically with small amplitude, and then we choose the period of pulse vaccination as a bifurcation parameter to investigate the complex patterns of the infectious disease dynamics.
Next, we analyze the existence of an endemic periodic solution of model (4) near the disease-free periodic solution ( * ( ), 0) using bifurcation theory [34] and set the impulsive period as bifurcation parameter. For ease of narrative, we rewrite the model (4) According to the bifurcation theorem in [34], we can obtain the following result.  The existence of 0 has been indicated in Figure 1 and the proof of Theorem 3 is given in Appendix B. Theorem 3 shows that there exists an endemic periodic solution under some conditions, providing that the disease-free periodic solution becomes unstable. It follows from Theorem 3 that if > 0 and is close to 0 , the endemic periodic solution of model (4) is stable.
To investigate the complex dynamics that model (4) can have, we chose the pulse period as a bifurcation parameter and fixed all other parameters as those in Figure 4 for two different values. Figure 4(a) is a bifurcation diagram without resource limitation; that is, = 0, which was obtained by Shulgin et al. [9]. Figure 4(b) is a bifurcation diagram with resource limitation; that is, = 0.2. Comparing those two bifurcation diagrams we conclude that the nonlinear pulse vaccination can produce more complex dynamics than those for the model with linear pulse vaccination. Once the pulse period exceeds max defined in (22), both susceptible and infected populations can oscillate periodically with a large amplitude that corresponds to periodic outbreaks of epidemics. As the pulse period is further increased, a series of complex and interesting bifurcation phenomena are observed (see Figures 4(a) and 4(b)). Figures 4(b)-4(d) indicate that the dynamical behavior of model (4) is very complex, including period doubling bifurcation, chaotic attractors, multistability, periodic-adding, chaos crisis, and periodic windows. Meanwhile, bifurcation analyses (i.e., Figures 4(b)-4(d)) also indicate that the model (4) has several different attractors which can coexist for a wide range of parameters. For example, Figure 5 provides an example of two attractor coexistence when = 7.8. Multiple attractor coexistence indicates that the infectious outbreak patterns depend on the initial values, which may cause difficulties for infectious disease control.

Discussion and Biological Conclusions
In order to understand the effect of resource limitation, in particular lack of vaccine, on the transmission of infectious disease, we deliberately investigated the dynamical behavior of an SIR epidemic model which incorporates a nonlinear pulse vaccination strategy. To this end, we introduced a nonlinear form ( ) = ( max ( ))/( ( ) + ) as vaccination proportion. This resulted in interesting and dramatic changes in the dynamical behavior of solutions and they became more and more complicated, which means that it is a very difficult task to control infectious diseases, in particular emerging infectious diseases, under resource limitation. We discussed the control strategy based on the threshold value in more detail through theoretical analysis, numerical studies, and sensitivity analysis. As pointed out by Agur et al. [7], the costs, risks, and the efficacy of the pulse vaccinations are key factors in eliminating transmission of the measles virus. In order to depict those factors, the nonlinear resource limitations should be taken into account [9], and the main purpose of the present work is to formulate the disease model to address the effects of resource limitations on disease control.
The results indicate that vaccination as an approach to controlling epidemics must be committed to a long-term strategy, especially when resources arelimited. To understand the interactions amongst ( ) and ( ) which are important for successful control of the spread of infectious disease, we have proposed detailed modeling methods involving vaccination on ( ) and determining the most significant parameters for the basic reproductive number 2 0 by using LHS/PRCC uncertainty and sensitivity analysis techniques (see Figure 3). Our results indicate that we should pay more attention to the vaccination period and vaccination proportion which reduce the threshold value 2 0 or 1 0 to prevent the outbreak of disease. Our results clarify that it is more significant to improve immunization programs (i.e., vaccinations as mentioned in [15,16]) under conditions of resource limitation.
The results also indicate that the dynamic behavior of model (4) may be dramatically affected by small changes in the value of initial densities of susceptible and infected with resource limitation. Bifurcation diagrams shown in Figure 4 indicate that there are many forms of complexities in model (4), which are related to chaotic bands with periodic windows and attractor crises, and Figure 5 reveals that the occurrence of multiple attractors is a common property of the SIR model with resource limitation, which can help us to further understand the application of a nonlinear pulse vaccination strategy in an SIR model (for more information about bifurcation diagrams of impulsive control strategies, see [35,36]). Some complexities are related to the long-term behavior of population dynamics, characterized either by well behaving relatively simple or very complicated strange attractors. We also found that the routes to chaos are very complicated. For instance, with resource limitation there are several hidden factors that can adversely affect the control strategy. The increasing number of potential complexities predicted by the theory does not seem to make this task any easier. Nevertheless, identifying complicated, possibly chaotic, dynamics in population data has remained a major challenge in epidemiology studies. In particular, the suggestion that the dynamics of measles are chaotic has been much debated (e.g., see [37]) and, when seeking to establish whether measles dynamics in New York City were chaotic or not, Sugihara and May [38] omitted data from after 1963 when vaccinations were introduced, as the immunisations had altered the intrinsic dynamics. However, our results imply that even with immunisations, disease dynamics could be chaotic under some circumstances and so a reanalysis of the New York data including the vaccination period might be worthwhile.
In the present paper, we have focused on studying the effects of dynamical behaviors on the SIR model with nonlinear pulse vaccination. For our simple model, we have assumed that the maximal vaccination proportion max is a constant and the pulse period can be varied. However, due to limited vaccine availability, the maximal vaccination proportion max should depend on the vaccine obtained in each pulse period. Thus, to describe the impact of a limited vaccine stockpile, a separate differential equation for the amount of vaccine, denoted by ( ), could be included in the model, and then ( ) should be a function of ( ) and ( ). Meanwhile, it is essential to link the costs of developing and implementing controls to population dynamic modeling of disease epidemics in order to consider other resource limitations such as quantities of drugs, availability of isolation units, numbers of hospital beds, and equipment. These topics will be considered in further work in the future.

A. Proof of Theorem 2
Firstly, we prove the local stability of a periodic solution ( * ( ), 0), which may be determined by considering the behavior of small amplitude perturbations of the solution. Define here ( ), ( ) are small perturbations, which may be written as where the fundamental matrix Φ( ) satisfies and Φ(0) = is the identity matrix. So where is not necessarily computed in detail as it is not required in the following analysis.
The resetting impulsive conditions of model (4) become ≐ ( ) ( ( ) ( ) ) . (A.5) Hence, according to Floquet theory, if the moduli of both eigenvalues of the matrix are less than one, then the periodic solution ( * ( ), 0) is locally stable. In fact, two Floquet expressions multiplied are thus It is obvious that 0 < 1 < 1. Thus, the stability of ( * ( ), 0) is decided by whether 2 < 1 or not. So we conclude that the disease-free periodic solution which holds true due to 1 0 < 1. Therefore, we can draw a conclusion from the above analysis that if 1 0 < 1, the disease-free periodic solution ( * ( ), 0) of model (4) is locally asymptotically stable.

B. Proof of Theorem 3
In order to apply the bifurcation theory of [34] in the main text, we make the following calculations: which indicates that there exists a 0 such that 1 0 = 1 and the disease-free periodic solution = ( * ( ), 0) loses its stability.
Further, It follows that < 0 and we get the condition, that is, if the parameters satisfy condition < 0, then the model (4)  Obviously, < 0 due to (25). The proof is completed.