SATURATED TREATMENTS AND MEASLES RESURGENCE EPISODES IN SOUTH AFRICA: A POSSIBLE LINKAGE

We consider the case of measles in South Africa to show that an high vaccination coverage may be not enough alone to ensure measles eradication. The occurrence of certain epidemic episodes may in fact be encouraged by delays in the treatments or by not adequately fast clinical case management, which may be related to the backward bifurcation phenomenon as well as to an intriguing spiking dynamics which appears in the system for specific ranges of parameter values.

1. Introduction and research motivations.Epidemic models are mostly expressed as nonlinear dynamical systems.One of the most striking features of nonlinear dynamical systems is the large number of different dynamical behaviors that may be obtained -even from relatively simple models -by varying significant parameters in the system.Bifurcation theory is a fundamental tool to get insight such variety of behaviors since, detecting the specific type of bifurcation which occurs in a certain system, provides the ability to recognize the dynamical behavior in the neighboring of the bifurcation point as well as the threshold values separating different dynamical regimes.Kermack and Mc Kendric [29] -just in 1927 -developed ideas as epidemic threshold or endemicity persistence.Today the concept of epidemic threshold has become one of the cornerstones of theoretical epidemiology and is intimately related to the basic reproduction number R 0 which commonly indicates the expected number of new infections produced by a single infective individual introduced into a diseasefree population [44].R 0 is hence a clear indicator of the strength of the disease and provides a direct measure of the control effort required to eradicate the disease.As a matter of fact, it may be considered the single-most useful quantity to estimate when modeling the population dynamics of an infectious disease.
For an endemic infection, the value R 0 = 1 defines a threshold, separating the stability and instability regimes of the disease-free equilibrium.In many cases, condition R 0 > 1 corresponds to disease endemicity whereas condition R 0 < 1 ensures disease eradication.However, since such threshold on the reproduction number is based on conditions at the disease-free equilibrium, it does not necessarily indicate the ability of the disease to persist at endemic levels [44].

DEBORAH LACITIGNOLA
Commonly two are the possible bifurcation scenarios at R 0 = 1: (i) forward bifurcation (ii) backward bifurcation.The former is surely the simplest from the dynamical point of view and implies disease eradication below the threshold R 0 = 1.
The backward bifurcation scenario involves instead multiplicity of endemic equilibria and subcritical persistence of the disease.This phenomenological framework includes a saddle-node bifurcation at R 0 = R sn 0 < 1 and a subcritical transcritical bifurcation at R 0 = 1.In this case, the classical method to reduce R 0 below 1 is not sufficient to eradicate the disease and a further effort should be done until R 0 is lowered below the critical value R sn 0 .In the perspectives of disease control, checking the occurrence of backward bifurcations has hence a primary importance.
A way to get insight on disease persistence is to perform a bifurcation analysis in order to establish the nature of the bifurcation at R 0 = 1.This goal may be achieved through specific bifurcation approaches [8,14,30,31] which make use of two powerful tools of dynamical systems theory: center manifold and normal forms [20].Such techniques often require preliminary transformations to be applied and tend to be technically enough complicated.More recently [3] has proposed a simple approach -based on algebraic elementary techniques -to get insight into the bifurcation behavior near R 0 = 1, without making use of the above specific bifurcation methods.While the center manifold and normal forms theory -even through complicated technicalities -are however the only possible choice to study bifurcation behavior of models with many compartments or strong nonlinearities, the simpler approach proposed by [3], turns out to be effective in many commonly used epidemic models.
Although many papers have been focused on detecting backward bifurcations in a variety of epidemiological models, few papers in literature have effectively investigated the role of the backward scenario in concrete disease-control settings, e.g.[6,19,22].In this paper we choose to consider the case of measles in South Africa to show how situations of saturated treatments may be responsible for measles resurgence episodes through the following two mechanisms: (i) subcritical persistence of the disease, which is related to the backward bifurcation phenomenon (ii) a spiking dynamics emerging in the system for specific ranges of parameter values when R 0 > 1.
Two questions might arise at this point: why measles?And why South Africa?Different motivations have supported this choice.Measles is a disease -mathematically well described by SEIR models [1] -for which endemic equilibria have been observed in many places, frequently with sustained oscillations about the equilibrium.Although its global incidence has been substantially reduced with vaccination, it remains an important public health problem, especially in developing countries.However -even in industrialized countries -several resurgences of measles have been observed in the last decades, despite the high vaccination coverage.A clinical study of measles in Poland reports an epidemic outbreak between 1997 and 1998 (with 2255 cases), despite a 95% vaccination coverage since the 1980s, [27].Persistence of measles in spite of the high vaccination coverage, raises questions about the existence of further mechanisms in order to explain why reducing R 0 below the transcritical bifurcation threshold may fail to control the spread of the disease.Saturated treatment mechanisms may offer interesting answers in this direction.
On this line, we have focused on South Africa because (i) it exhibits an high overall national rate of measles vaccination coverage (ii) despite of this, measles resurgences have been recently observed, e.g. the outbreaks of 2003 and 2009 (iii) investigations during recent measles outbreaks have revealed the considerable need to improve therapeutic treatments because vaccination of contacts alone was shown to be ineffective in stopping transmission of measles in resource-limited settings [36].We also think that using mathematical models with parameters and data proper of the South-African context is of a certain importance since there are still few measles studies performed with parameters proper of the developing countries, principally due to a lack of data [4].In South Africa there is a good availability of data thanks to the NICD (National Institute for Communicable Diseases), the national organ for public health surveillance of communicable diseases born with the aim to collect, analyze and interpret communicable diseases data on an ongoing and systematic basis.Data used in our case-study are mainly coming from NICD sources, as the Communicable Diseases Surveillance Bulletin.
The rest of the paper is structured as follows.In Section 2 we introduce a SEIR model potentially suitable for the case study, which we show to be structurally equivalent to the model considered in [53].In Section 3 we recall its basic properties and revisit the problem of subcritical endemicity showing that backward bifurcation conditions may be very easily obtained through the simple approach proposed in [3].In Section 4, we explicitly consider the case of measles in South Africa.The obtained theoretical results are validated by considering epidemiological parameters proper of the South African context in the case of slow treatments, interestingly showing that backward scenario may be detected within such field-parameter setting.The considered model seems however to be dynamically richer than this.In fact, we numerically show the existence -for a certain range of the parameter values -of a spiking attractor involving the disease free equilibrium when R 0 > 1, which seems to be related to an homoclinic bifurcation.In Sections 5-6, the role of backward scenario and spiking mechanism as direct fruits of saturated treatments is discussed in order to provide operative solutions to the problem of measles resurgence episodes in South Africa.
2. The SEIR model.We consider the following SEIR model with treatment, describing the dynamics of susceptible (S), latent or exposed (E), infectives (I) and removed (R) individuals: where S(t) + E(t) + I(t) + R(t) = N (t).The involved parameters, all positive, have the following meaning: A denotes the influx or recruitment of susceptibles, p is the proportion of the recruited individuals who are vaccinated, i.e. 0 ≤ p < 1, d is the natural death whereas α is the disease-caused death; β is the contact rate; and γ denote the transfer rates between the corresponding compartments.Heuristically, 1 stands for the mean latent period and 1 γ for the mean infectious period.The requirement of a constant recruitment A of susceptibles has a particular sense in the perspective of using model (1) to analyze endemic diseases in the context of less developed countries.In fact to model a disease which may be endemic, both the shorter epidemic time scale and the longer demographic time scale have to be considered.This means that births and deaths have to be included in the model.A simple way to allow this, is to assume that the number of births per unit of time equals the number of deaths, so that the total population size is constant.However, such an assumption is plausible if there are no deaths due to the disease (α = 0) but it is not appropriate in presence of significant disease-caused deaths (α = 0).In less developed countries, where there is often a very high mortality rate for certain kind of infectious diseases, the total population must vary in time: this may be simply obtained by requiring a constant recruitment rate.
In model (1), F (I) represents the treatment term.In classical epidemic models, the treatment rate of infectives is assumed to be directly proportional to the number of infectives and hence modeled through a linear function.In practice, this feature is enough unsatisfactory because the resources required for treatment should be quite large.On the contrary, each city or country has its maximal capacity for the treatment of a disease.At this regard, in [46], a staged treatment function is introduced for which the treatment rate is proportional to the number of infectives when the capacity treatment is not reached and takes the maximal saturated level otherwise.This accounts for example situations where -to treat patients -the number of hospital beds is limited or medicines are not sufficient.
A similar idea was proposed in [11] and [51], where saturation recovery of infectives was modeled through a smooth Verhulst-type function [38].We adopt this idea and consider the following treatment function: where r > 0 is the treatment rate and b 0 > 0 measures the extent of the effect of the infectives being delayed for treatment.To show how the saturation coefficient b 0 may be related to delays in the treatment, we derive the essential features of the treatment term (2).At this aim, we assume the total time T spent for treatment of infective individuals as divided in two main kinds of activities: (i) collecting infectives (ii) managing infectives.Hence T = T collect + T manage , where T collect and T manage are respectively the times spent for collecting and managing infectives.Now, denoting with I c the number of infectives collected to be cured during the time T , we can also assume I c as proportional to the time spent to collect them: As a consequence, T collect = I c r I .On the other hand, it is also plausible to assume that the management time T manage is proportional to the number of infectives I c collected to be cured: T manage = I c T h , where T h is the time spent on managing one infective individual.It thus follows Since the number of infective individuals treated per unit of time is I c = F (I) T , one has: As a consequence (2) is obtained with the saturation coefficient b 0 given by b 0 = rT h and hence related to the time spent in managing the infectives.Higher values of b 0 may hence be related to higher collecting rates as well as higher times spent in managing individuals.Treatment term (2) is thus a plausible choice for modeling situations when the number of infectives is getting larger and the medical condition is limited; it also explicitly includes the effect of delayed treatments.This is an important point to consider since the efficiency for treatment may be seriously affected if the infective individuals are delayed for treatment.
With regard to disease transmission, model (1) implicitly assumes that new cases are generated through homogeneous mixing (between susceptible and infective individuals) so that the mass-action incidence term βIS is considered.The hypothesis of homogeneous mixing may be inaccurate, particularly under specific circumstances, i.e. situations including saturation or thresholds effects.In this context we have used a mass-action term to better elucidate the contribute of the treatment term for the occurrence of the backward bifurcation phenomenon.We are however aware that a more general contact rate would have provided a more reliable model.3) hence exhibits the same mathematical structure of the general SEIR model with constant recruitment introduced in [32], enriched with the saturated treatment function modeled in [51].
A number of studies have shown that multigroup structure with asymmetry between groups, multiple interaction processes (e.g.demographic and epidemiological ones), pair formation, macroparasite infection, age structure, are all mechanisms that can cause backward bifurcations, [2,14,23,25,26,34].
In this line, a certain attention has been recently devoted to the role of the treatment mechanism: saturated-type treatment functions have been indicated as responsible for the occurrence of backward bifurcations in SIR [46,51,52] as well as in SIS [11,47] and SEIR models [53].
At this regard, model (3) may be considered a clear example supporting the general circumstance that saturated-type treatments may be one of the causes of backward bifurcation in epidemic models.In fact investigations performed in [32] about the SEIR model (without treatment) indicate the occurrence of a very classical forward scenario at R 0 = 1: global stability for the disease-free equilibrium in the range R 0 < 1 and global stability of the endemic equilibrium in the range R 0 > 1. Hence the occurrence of backward bifurcation in the same model enriched with treatment, means considering the backward scenario as a direct fruit of the saturated treatment mechanism.Moreover, the SEIR model (3) with treatment mechanism ( 2) is structurally equivalent to the SEIR model in [53] where local and global equilibria stability properties as well as the occurrence of backward bifurcation have been investigated.We here briefly revisit some of their results in view of the model application in Section 4.
3.1.Analytical results.For the sake of simplicity, in the rest of the paper we use model ( 3) with dropping the prime notation.The first three equations in system (3) are independent of R and one considers the reduced system: whose dynamics can be studied in the feasible region Γ = {(S, E, I) ∈ R 3 + : S + E + I ≤ A d }, which is positively invariant with respect to (4).A first effect that the treatment mechanism (2) produces on the dynamics of model ( 4), is the chance to have multiple endemic states.In fact introducing the basic reproductive number, we observe that system (4) admits (i) the disease-free equilibrium P 0 = A d , 0, 0 which is locally asymptotically stable for R 0 < 1 and unstable for R 0 > 1 (ii) the endemic equilibrium P = (S * , E * , I * ) where and I * represent the real positive solutions of the algebraic equation Here, ) Looking at (7) and observing that it follows that, when R 0 > 1 then C 1 < 0 and hence ∆ = B 2 1 − 4A 1 C 1 is a positive quantity.According to the Descartes's rule of signs, the algebraic equation ( 6) admits a unique real positive root, independently of the sign of B 1 .When R 0 < 1, the matter is more complicated and is worthy of separate investigations as detailed in Appendix A. Let A * , β * 2 and β * be the quantities respectively defined in ( 13),( 11), (8).The following result holds: 4) admits the disease-free equilibrium and no endemic equilibrium.If β * 2 < β v < β * , then model ( 4) admits the disease-free equilibrium together with two endemic equilibria.If β v > β * then model ( 4) admits the diseasefree equilibrium and a unique endemic equilibrium.
The above framework clearly indicates a backward bifurcation scenario, suggesting the occurrence of a saddle-node bifurcation at 2 and of a transcritical bifurcation at 3.2.Backward vs forward scenario revisited.Differently from [53] where the occurrence of backward bifurcation has been established through the bifurcation method introduced in [8],we make here use of the simple approach proposed in [3] which is based on an accurate glance at the equilibrium curve in the neighboring of the bifurcation point R 0 = 1.
Recalling (5), it is easy to observe that all coefficients in equation ( 6) may be regarded as functions of the parameter β v , i.e.
. By (7), it follows that when R 0 = 1 then C 1 = 0, i.e.C 1 (β * ) = 0. Hence equation ( 6) becomes: which admits as roots I = 0 and I = − B1(β * ) A1(β * ) .The root I = 0 is related to the disease free equilibrium whereas I = − B1(β * ) A1(β * ) corresponds to a positive (endemic) equilibrium only if B 1 (β * ) and A 1 (β * ) have opposite sign and hence B 1 (β * ) < 0 must hold.Further, implicit differentiation with respect to β v of the equilibrium curve (6) leads to: Now, looking at the equilibrium I = 0 at β v = β * , one obtains: since by (7) the quantity dC 1 dβ v is negative.One may hence conclude that the slope of the bifurcation curve at I = 0 has the same sign as the coefficient B 1 (β * ).As a consequence, if B 1 (β * ) < 0 then a backward bifurcation occurs at In our case condition B 1 (β * ) < 0, which is necessary and sufficient for the occurrence of the backward bifurcation at β v = β * , is given by which means: The above results are summarized in the following theorem: This feature strongly enlightens the role of the delayed treatment for the occurrence of backward scenario.In fact, when the effect of the infectives being delayed for treatment is weak, i.e. when b < b * , then R 0 has to be classically lowered below 1 to eradicate the disease.On the contrary, when this delayed effect is strong, i.e. when b > b * , driving R 0 below 1 is not enough to eradicate the disease.4. Measles in South Africa: A case of study.The above results may be fruitfully used to get insights on recent measles outbreaks in South Africa.Measles is a highly contagious, vaccine preventable, viral disease caused by paramyxovirus of the genus Morbillovirus and transmitted through aereosol particles.The incubation period lasts about 8 days after which individuals are infectious for nearly 5 days.Typical signs and symptoms include fever, a non blistering rash, coryza and conjunctivitis, [24].Although an effective vaccine is available, measles remains among the leading causes of vaccine preventable deaths in children under five years of age, especially in developing countries, [48].
In South Africa, routine immunization campaigns were first introduced in 1975 and measles has been a notifiable disease since 1979.In the 1980's, from 15000 to 20000 measles cases were counted each year and the reported measles-caused deaths ranged annually from 250 to 500.During the 1990's, measles remained endemic and epidemic continued to occur periodically even if at the beginning of the 1990's, the case fatality ratio (CFR) sharply declined and stabilized at low levels, [43].Since 1995 seven southern South Africa nations, including South Africa, launched measles elimination initiatives -in accordance with the recommendations of the World Health Organization (WHO) Regional Office for Africa [9] -which resulted in the virtual elimination of measles in Southern Africa, [36] The main goals of such recommended programs were (i) to achieve and sustain routine immunization coverage of 95% (ii) to implement a one-time national 'catchup' measles vaccination campaign (iii) to implement periodic national 'follow-up' campaigns (iv) to establish case-based measles surveillance with laboratory confirmation [9].The NICD was accredited by WHO to perform measles IgM testing for national case-based surveillance, as part of the measles elimination strategy.In South Africa, a first and a second dose of measles vaccine (MCV1, MCV2) are provided to children at 9 and 18 months respectively.Moreover, since 1996, South Africa has been at the forefront of implementing regular national or subnational supplementary immunization activities (SIAs), following an approachimplemented by the Pan American Health Organization -which is thought to have contributed to measles eradication in the Americas [13].Over the period 1996-2010, SIA coverage in South Africa was executed in each of the country's nine provinces and has remained at high levels, which indicates that a considerable amount of resources have been employed for the SIA policy, [45].In 1996-1997, the overall measles vaccination campaign coverage was estimated at 85% [43] whereas in 2000, a nationwide campaign reported 92% administrative coverage.Measles incidence reflected the success of the above campaigns: in the period 1999-2002, less than 60 measles cases were reported annually, [36].Nevertheless in 2003, a measles virus -introduced to Mpumalanga and Gauting provinces from Mozambique -lead an epidemic lasting more than 2 years.To face such an emergence, a national measles vaccination campaign was conducted in 2004, with a consequent decreasing in the measles virus transmission, [36].For the period 2006-2008, numbers of measles IgM positive cases remained relatively low.However, a widespread outbreak of measles occurred in 2009.
This seems to indicate that even in a setting of low reported cases, communities still remain vulnerable to measles importation so that it is always necessary to maintain immunity among the population even in the absence of a circulating virus.Efforts to maintain the vaccine coverage at 95% -as recommended by WHO -have hence been emphasised [10].However such a measure, alone, might be not enough.
We will use model ( 4) with epidemiological parameter values proper of the South African context, to show that even in the presence of an high routine vaccination coverage, the occurrence of some epidemic episodes may be encouraged by features -as the not adequately fast clinical case management -that might be related to backward bifurcation or to spiking phenomena.4.1.Parameters setting.As a first step, we aim to choose parameter values so that model (4) might adequately describe the South African case.

Demographic indicators
For the choice of demographic parameter values, we explicitly refer to [42].In 2009, the crude birth rate (i.e. the annual number of births per 1000 population) for South Africa was set at 22. We then assume the recruitment A of susceptible individuals to be A = 22000 [individual][year] −1 discounted by a 95% vaccination rate (p = 0.95) The crude death rate (i.e. the annual number of deaths per 1000 population) was instead established at 15.We then choose the death rate d to be d = 0.015 [year] −1 .

Measles indicators
The mean latent period for measles was recognized to be 8 days.Since we choose year as time unit, we can assume the parameter -representing the rate at which exposed individuals become infectious -as = 365/8 [year] −1 , that is = 45.6250[year] −1 .The mean infectious period is 5 days.Hence we can assume the parameter γ -representing the rate at which infective individuals recover -as γ = 365/5 [year] −1 , that is γ = 73 [year] −1 , [24].Finally, measles mortality rate -mortality rate per 1000 live births -is set at 168 so that we assume the disease-caused death rate α = 168/1000 = 0.168 [year] −1 , [49].

Treatment indicators
South Africa, has a well-developed, predominantly hospital-centered health care system, in which approximately 80% of the population receives their health care, mainly in the public sector, [43].Accordingly, we assume the treatment rate to be r = 0.8 [year] −1 .
The WHO and the United Nations Children's Fund (UNICEF) Integrated Management of Childhood Illness guidelines, recommend for measles treatment, antibiotic therapy and two doses of Vitamin A given 24 hour apart.Moreover, a prompt isolation of measles cases strongly limits nosocomial measles transmission.Investigations during measles outbreaks in South Africa [36] have however revealed that less than half of the children presenting to health facilities reported any treatment with vitamin A and that isolation measures in health facilities were inadequate.Motivated by the above field observations, in the sequel we consider the contact rate β and the saturation coefficient b as bifurcation parameters, and exploit the analytical results of Section 3 to study the dynamics of system (4) with the parameter values described above and summarized in Table1.1. Parameters values used in simulations and the related sources.

4.2.
A qualitative model testing.The chosen parameter values allow to estimate the basic reproductive ratio during the 2009's outbreak in South Africa, by using the theoretically derived relationship R 0 = g/ā where g is the reciprocal of the per capita birth rate and ā is the mean age at first infection [1,4].For values as in Table1, g = 1/0.022= 45.4545whereas the parameter ā may be derived by the reported ages of the infectives which provide us with a field estimate of the mean age at infection of nearly 11.2927 years, Fig1.This provides a field basic reproductive number R 0 = 4.0251.To qualitatively validate model (4), we illustrate its dynamical behavior when the basic reproductive number is R 0 = 4.0251 and p = 0.85.This choice for p is since the reported vaccination coverage for South Africa in 2008, was 85%, which is less than the target of 90% that forms part of the measles elimination strategy.According to the reported cases [36], we also consider enough slow treatments for the infectives and choose b = 15.This point surely deserves a comment since it requires the choice of a specific value of the parameter b, which is not available in literature.Such numerical choice for b has been performed according to the considerations summarized in Fig. 2. For parameter values as in Table1, the threshold b * separating the backward and forward scenario is b * ≈ 0.3.Recalling that b ∝ rT h , we observe that higher values of b may be related, for a fixed treatment rate r, to higher times T h spent in managing individuals.We qualitatively consider as very fast treatments the ones for which b < b * = 0.For the purposes of the backward bifurcation, the effect of b-values higher than 30 are similar enough to those obtained for b = 30.Hence we will consider 0 < b ≤ 30.For the above field-value of R 0 , bifurcation analysis prescribes endemicity of the disease.Analysis in the time dependent regimes for initial conditions near the endemic equilibrium are shown in Fig. 3 and reveal for the damped oscillations, a period of oscillation T ≈ 5 which may be in a certain accordance with field data estimation.Time dependent investigations for the same parameter values and initial conditions in the neighboring the disease-free equilibrium are depicted in Fig. 4 which shows how trajectories tend asymptotically to the endemic equilibrium but also reveals a peculiar structure of the transient oscillations.We will come back on this point at the end of this section.

4.3.
Disease subcritical persistence.Basing on the WHO recommendations inside the measles elimination programs, in the sequel we increase the vaccination coverage p at the recommended value p = 0.95 to show that -despite of the high vaccination coverage -subcritical persistence of the disease via backward bifurcation, may be related to sufficiently high delays in the treatment of the infectives.
At first, we provide some explicit remarks on the different thresholds of model (4) considered with the above chosen parameter values.The threshold b * , separating the backward and the forward scenarios, is b * = 0.311095 whereas the basic reproductive number is given by R 0 = 999.8929β.R 0 does not depend on the value we choose for the parameter b as well as the transcritical threshold value β * ,    1. where . In this range, two endemic equilibria P * 1 and P * 2 exist which coalesce at β = β * 2v via saddle-node bifurcation, Fig. 6.With regard to stability properties, the endemic equilibrium P * 1 -characterized by a small number of infectives -is always unstable whereas the disease-free equilibrium P 0 is locally asymptotically stable for β < β * v .The endemic equilibrium P * 2 -characterized by a large number of infectives -may instead change its stability properties by varying the bifurcation parameter.4.4.Slow treatment settings and disease control.Suppose now to consider a situation involving a fixed, sufficiently slow treatment b and investigate the effects produced on the dynamics of system (4) by varying the contact rate β.  1.
Table 2. Feasible equilibria and the related stability properties for b = 15 and β as bifurcation parameter.The numerical values for the other parameters are as in Table 1.
By further increasing the parameter β, the endemic equilibrium P * 2 becomes an unstable focus, so that initial conditions near P * 2 make the system approach the disease free equilibrium P 0 through oscillations of increasing amplitude, Fig. 7(a).At β = β H ≈ 0.0010029, the endemic equilibrium P * 2 changes its stability because of a subcritical Hopf bifurcation.Hence, for β > β H , P * 2 is a stable focus: initial conditions near P * 2 lead system's trajectories approach P * 2 with damped oscillations.Fig. 7(b) well underlines how trajectories starting near P * 2 , initially seem to follow the unstable periodic orbit in the surrounding of such equilibrium before approaching P * 2 , so that the transient nearly resembles a periodic motion.Existence and stability properties of the equilibria are specifically summarized in Table2.These results put into evidence that -although multiplicity of endemic equilibria are possible in the range β * 2v < β < β * v -when treatments are sufficiently slow, subcritical persistence is really possible only in the range β H < β < β * v .This should suggest that -to control disease -it might be enough to lower β below β H .However, in this case the disease-free equilibrium -which is the only attractor for the system -might be approached via oscillations of large amplitude.Hence a more effective way to eradicate disease turns out to lower β below the saddle-node bifurcation threshold β * 2v .Such findings well elucidate the key role of accurate timedependent investigations for gaining correct informations on disease-thresholds in way to efficaciously perform the most suitable control measure.This aspect is more strongly confirmed by looking at the range β ≥ β * v .For β = β * v the transcritical bifurcation takes place so that the disease free equilibrium P 0 changes its stability becoming unstable.In this range the endemic equilibrium P * 2 is locally asymptotically stable and initial conditions near P * 2 make system trajectories to damply oscillate toward P * 2 .However a different set of initial conditions may produce a different output.Numerical simulations show in fact that initial conditions in the neighboring of the disease-free equilibrium P 0 lead system trajectories toward a spiking attractor, represented in Fig. 8 and Fig. 9 where temporal diagrams and the related phase ports are plotted for values of β increasing from β * v .It clearly appears how the spiking effect increases with increasing the values of β whereas for β → β * v the spiking attractor is likely to culminate into an homoclinic orbit with beginning and end at P 0 .The spiking attractor obtained for initial conditions near P 0 when β = 0.00106 is specifically reported in Fig. 10 for different time intervals, in order to better elucidate the spiking mechanisms.Such an attractor seems to loose stability for β = 0.002 so that initial conditions near P 0 lead system trajectories toward the stable endemic equilibrium P * 2 , Fig. 11.Numerical investigations in the case β > β * v have hence revealed an intriguing bistability situation in the range β * v < β < 0.002 where -according to the initial conditions -the endemic equilibrium P * 2 or the spiking attractor involving P 0 may be achieved.The former dynamical case may be epidemiologically interpreted as persistence of the disease whereas the latter corresponds to a cyclic apparent eradication of the disease with sudden epidemic bursts since it occurs superthreshold when initial infective population size is close enough to zero.Moreover the frequency of the bursts increases with β increasing from β * v .We also stress that such bistability situation appears to be the fruit of including saturated mechanisms in the model since, for the case b = 0, the endemic equilibrium P * 2 is a globally stable attractor for the system, as shown in [32].Theoretical investigations to accurately get insight on the arising and the nature of the above spiking attractor will be developed in future studies.

5.
Discussions.The importance of detecting backward bifurcation is a matter of fact in terms of disease management and control.Differently from the forward scenario -for which the disease control appears relatively smooth -the backward phenomenon requires a precise strategy for disease eradication.The transcritical threshold is not a strict threshold for the disease control and the bifurcation parameter has to be lowered below the saddle-node bifurcation value in order to avoid subcritical persistence of the disease.As a consequence, the skill to connect the backward phenomenon with the occurrence of specific biological mechanisms -as saturated-type treatment of the infectives -may be a great practical advantage to set the most suitable disease elimination strategy.
Necessary and sufficient conditions for backward bifurcation are here obtained for model (4) by exploiting considerations on the equilibria curve at R 0 = 1 as  recently suggested in [3].The full agreement of such findings with the general results obtained in [53] clearly reveals that such simple approach may be considered as an efficacious and more direct alternative to the classical bifurcation techniques, for models with not many compartments and without strong nonlinearities.
We have applied the above theoretical results to the case of measles in South Africa, for which certain saturated treatment mechanisms have been identified.Field investigations during measles outbreaks in South Africa [36] have in fact revealed that less than half of the children presenting to health facilities reported any treatment with vitamin A and that isolation measures in health facilities were nor prompt nor adequate.As a consequence, the strong need for more extensive therapeutic use of vitamin A and the greater use of respiratory isolation have been highly stressed since vaccination of contacts alone was shown to be ineffective in stopping transmission of measles in resource-limited settings.
The study of model ( 4) with parameter values proper of the South African context, has revealed that -although in presence of an high vaccination coverage, sufficiently high delays in the treatment of infectives may encourage the occurrence of a backward bifurcation phenomenon and lead to possible subcritical persistence of the disease.Time dependent investigations clearly elucidate the role of b on the dynamics of system (4) and strictly confirm the importance to strongly reduce delays in the disease management.
Our results show that a very effective method to prevent the backward scenario is to reduce the parameter b which means -in practice -to give patients very timely treatments with improving medical technologies or investing much more in medicines, beds, etc.In this case the control of the disease may straightly follow the classic road of reducing the contacts β below the transcritical threshold β * v .On the contrary if not enough fast treatment may be provided, then the backward scenario cannot be avoided and has to be properly managed.In this case a reduction of the contact rate β under the saddle-node threshold β * 2v is required to obtain the eradication of the disease.At this aim, a more adequate use of isolation measures in health facilities as well as specific public education programs may be employed since properly reducing contacts they can reduce disease transmission.
Time dependent investigations have also revealed that saturated treatment mechanisms (b = 0) may be responsible for an intriguing bistability situation above the threshold β = β * v where in a specific range of the parameter β, system trajectories tend either toward the stable endemic equilibrium or to a spiking attractor.Hence, according to the initial population sizes, disease may persist at an endemic level or exhibit recurrent epidemic bursts.
Looking at the peculiar structure of the spiking attractor, it appears reliable to conjecture that such latter scenario might be the effect of an homoclinic bifurcation involving the non hyperbolic disease-free equilibrium at β = β * v .The homoclinic bifurcation is a global bifurcation, not easy to be detected since it is not revealed by local stability analysis.In an homoclinic bifurcation of an (hyperbolic/non hyperbolic) equilibrium, a periodic attractor emerges from an homoclinic orbit, i.e. a trajectory that tends to the same equilibrium both in direct and in reverse time [17].For their own structure, homoclinic bifurcations can epidemiologically be related to sudden epidemic bursts and -as well as backward bifurcations -they join the group of the so called dangerous bifurcations.A dangerous bifurcation implies discontinuous behavior namely sudden appearance or disappearance of the attractor and a fast dynamic jump toward a different, unrelated attractor.
Such a feature clearly reveals how homoclinic phenomena may have a remarkable potential for applications in epidemic context and justifies the increasing theoretical interest for homoclinic bifurcation in epidemics [12,15,40,50].Nevertheless, to the best of our knowledge, there is no example of concrete linkage between homoclinic phenomena and epidemic bursts in a field-case of study.At this regard the case of measles in South Africa is worth to be deeply investigated since it might interestingly offer such kind of example.
6. Concluding remarks.In conclusion, the application of model ( 4) in the South African context strongly suggests that the occurrence of certain epidemic episodes might be discouraged through an overall improvement of the health system as a whole.On this line the backward bifurcation phenomenon along with the spiking mechanism may help in the understanding of that complex picture which sees measles resurgence episodes in South Africa -like the outbreaks of 2003 and 2009 -as open challenges to the basic measles control in a country characterized by an high overall country rate of vaccination coverage.In fact, giving a satisfactory solution to the problem of measles eradication in South Africa is far from being a straightforward issue since the interplaying of many factors has to be properly taken into account.As well indicated in a recent study on measles control measures in South Africa [45], the high national rates of measles vaccination coverage masks a substantial heterogeneity across the 9 provinces and 52 districts which ultimately facilitates the resurgences of outbreaks.
Verguet and coauthors also stress that, although supplementary immunization initiatives can reduce part of this heterogeneity with raising the overall coverage in all the districts, SIAs can also have the collateral effect of affecting performance of the health system as a whole by interfering with delivery of other health services which may be interrupted during SIAs because of staff shortage and inadequate preparation [45].On the other hand in a resource-limited setting, possible malfunctioning in the health services are often associated to saturation mechanisms and may favor delays in treatments as well as not adequately fast case managements, and hence subcritical persistence of the disease via backward bifurcation or spiking dynamics via likely homoclinic phenomena.
Our findings hence outline one of the reasons why South Africa remains vulnerable to large measles outbreaks and suggest possible prevention measures which are in accordance with field results: focusing on rapid case notifications and on adequate and fast clinical case management including provision of vitamin A and antibiotics where indicated; strongly improving the use of isolation measures in health facilities and -last but not least -keeping the high recommended vaccination coverage by carefully maintaining at their best health system services, especially during SIAs campaigns.
To deeper investigate on conditions (i)-(ii), we preliminary introduce the following quantities: We also observe that equation ∆ = B 2 1 − 4A 1 C 1 = 0 may be regarded as where: B 2 = −2dΓ 2 [d (d + γ + α + + r) + (α + γ + bA + r)] ; Moreover, ∆ 2 = B 2 2 − 4A 2 C 2 is such that ∆ 2 = 16d 2 b Ar( + d)Γ 2 2 .∆ 2 is a positive quantity and coefficients (10) always have definite sign.As a consequence, according to the Descartes's rule of signs, equation ( 9) admits two real positive roots, β * 1 and β * 2 : It follows that ∆ = B 2 1 − 4A 1 C 1 > 0 if and only if β v < β * 1 or β v > β * 2 .We now explicitly search for conditions on the system parameters, ensuring that B 1 < 0 in the range R 0 < 1.This is certainly true if where A = (d + γ + α + r)( + d) b , β = dΓ 2 Γ 1 and β * -already defined in (8)rearranged as β * = Ab − Γ 1 A .Note that (12).1 ensures Γ 1 to be a positive quantity whereas (12).2 prescribes B 1 to be negative in the range R 0 < 1.Moreover, in order (12).2 to hold, we must require that It is also easy to note that A * > A so that A > A * ⇒ A > A. By combining the above observations with (12), we can conclude that in the range R 0 < 1, the coefficient B 1 is a negative quantity if The above results may be summarized as follows: The last step is hence to make order among the above β v -thresholds. .But β > β * is equivalent to R 0 > 1, so that no real root of equation ( 6) should exist in this range of R 0 .By the Descartes' rule of sign, we know it is absurd.Hence β * 2 < β * .

Remark 2. Let
Now, by considering in model (1) the change of variables S = S (1 − p), E = E (1 − p), I = I (1 − p) and R = R (1 − p) + (A/d)p, the dynamical equations for the primed variables are given by:

Figure 2 .
Figure 2. Qualitative explanation for the numerical choice of the parameter b, as related to fast or slow treatments.
3.Moreover, choosing b as b = 3 means that the related treatments are 10 times slower with respect to b * and we classify as not so fast the ones whose b * < b ≤ 3. Treatments whose b are such that 3 < b ≤ 15 are instead classified as slow since for b = 15 they are 50 times slower with respect to b * .Following the same reasoning, treatments for which 15 < b ≤ 30 are classified as very slow.

Figure 3 .
Figure 3.The case R0 = 4.0251, p = 0.85, b = 15.Initial conditions are chosen in the neighboring of the endemic equilibrium.The numerical values for the other parameters are as in Table 1 (a) Time dependent behavior (b) Phase diagram.

Figure 4 .Figure 5 .
Figure 4.The case R0 = 4.0251, p = 0.85, b = 15.Initial conditions are chosen in the neighboring of the disease-free equilibrium.The numerical values for the other parameters are as in Table 1 (a) Details of the temporal dynamics of the exposed population (b) Details of the temporal dynamics of the infective population.

Figure 6 .
Figure 6.Bifurcation diagram in the plane (b, β, I * ) showing the coalescence of the two endemic equilibria I * .The numerical values for the other parameters are as in Table1.

Figure 10 .
Figure 10.The attractor related to the case β = 0.00106 for different time intervals (a) Projection of the attractor in the (S, I) phase space; t f inal = 5000 (b) Projection of the attractor in the (S, I) phase space; t f inal = 10000 (c) (S, E, I) phase port; t f inal = 20000 (d ) Details of the temporal dynamics of the infective population.

Figure 11 .
Figure 11.The case β = 0.002.Initial conditions are chosen in the neighboring of the disease-free equilibrium P0.The numerical values for the other parameters are as in Table 1 (a) Time dependent behavior (b) Phase diagram.