Computational and Mathematical Methods in Medicine

Volume 2013 (2013), Article ID 648291, 19 pages

http://dx.doi.org/10.1155/2013/648291

Research Article

Modeling the Spread of Tuberculosis in Semiclosed Communities

Mauricio Herrera,¹ Paul Bosch,² Manuel Nájera,³ and Ximena Aguilera³

¹Facultad de Ingeniería, Universidad del Desarrollo, Santiago 7620001, Chile

²Facultad de Ingeniería, Universidad Diego Portales, Ejército 441, Santiago 8370179, Chile

³CEPS, Facultad de Medicina, Clínica Alemana, Universidad del Desarrollo, Santiago 7710162, Chile

Received 23 December 2012; Revised 11 March 2013; Accepted 11 March 2013

Academic Editor: Nestor V. Torres

Copyright © 2013 Mauricio Herrera et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

We address the problem of long-term dynamics of tuberculosis (TB) and latent tuberculosis (LTB) in semiclosed communities. These communities are congregate settings with the potential for sustained daily contact for weeks, months, and even years between their members. Basic examples of these communities are prisons, but certain urban/rural communities, some schools, among others could possibly fit well into this definition. These communities present a sort of ideal conditions for TB spread. In order to describe key relevant dynamics of the disease in these communities, we consider a five compartments SEIR model with five possible routes toward TB infection: primary infection after a contact with infected and infectious individuals (fast TB), endogenous reactivation after a period of latency (slow TB), relapse by natural causes after a cure, exogenous reinfection of latently infected, and exogenous reinfection of recovered individuals. We discuss the possible existence of multiple endemic equilibrium states and the role that the two types of exogenous reinfections in the long-term dynamics of the disease could play.

1. Introduction

Dynamics of tuberculosis (TB) spread has been the subject of a considerable body of theoretical and mathematical work. For review, see, for example, [1, 2] and references therein. The choice of a particular model is strongly connected to the questions we want to answer, and in the present work we will address the problem of long-term dynamics of tuberculosis and latent tuberculosis (LTB) in semiclosed communities.

For semiclosed communities we mean not strictly closed communities with certain mobility of their members out or into the community, with a recruitment of new members and departure of others. But, essentially these communities are congregate settings with the potential for sustained daily contact for weeks, months, and even years between community members. Basic examples of these communities are prisons, but certain urban/rural communities, schools, among others could possibly fit well into this general definition. These communities present a sort of ideal conditions for frequent TB outbreaks, enhanced TB transmission, and accelerated spread of the disease.

The basic characteristics of such settings including the possibility of high concentrations of infectious individuals and immunodeficient hosts, improper precautions taken for protection, delay in diagnosis, sustained contact with the index case, and inadequate ventilation and/or overcrowding make them well suited for TB transmission, creating this way genuine high transmission pockets of TB inserted in the general population [3, 4].

In fact, prisons are especially high burden communities, in which incidence and prevalence of TB are very high, and consequently the frequency of infections and reinfections considerably increases in comparison with population at large; see the works by Chiang and Riley [5] and by Baussano et al. [6].

Studying the dynamics of the TB spread in semiclosed communities is an interesting and significant topic by itself; however, there is an important phenomenon due to which the study of these types of communities is essential in the context of TB spread. This phenomenon has been called the Reservoir Effect [6, 7]. Indeed, semiclosed communities such as prisons represent a reservoir for disease transmission to the population at large and should be a source of public concern. For example, TB infection may spread into the general population through prison staff, visitors, and close contacts of released prisoners. The transmission dynamics between prisoners and the general population [6], together with immigration from developing countries with high prevalence of TB [8, 9], has been hypothesized to play a key role in driving overall population-level TB incidence, prevalence, and mortality rates.

In a recent work [4] the authors have even gone further in relation to this effect and have named these communities Institutional Amplifiers of TB Propagation. Some examples of communities given by these authors are poor hospitals in which dozens of patients share poorly ventilated communal rooms, crowded prison cell blocks, and mining barracks among others.

The transmission and progression of TB infection has been relatively well understood on a population scale. Generally, it is assumed that once an individual is infected with TB, he or she is immune from further infection events. Moreover, it was proposed what came to be known as the unitary concept of pathogenesis [10], which states that TB always begins with primary infection, and subsequent episodes of active TB are due to reactivation of dormant bacilli from this primary infection. However, a persistent evidence has recently been shown (see [5] for a review) that the paths to TB infection are not as linear as was suggested by the unitary concept of pathogenesis. The availability of individual, strain-specific infection histories (see, e.g., [11–13]) has made it clear that exogenous reinfection in people with previously documented TB infection does occur. The important question is whether reinfection occurs commonly enough to have an effect on the overall infection dynamics of the population [14].The relative importance of these pathways to the development of active disease has significant implications for treatment and control strategies, most notably in deciding whether latently infected and treated individuals are at risk of reinfection [15].

Several authors [15–20] have declared that exogenous reinfection plays an important role in the disease progression and that the inhalation of tubercle bacilli by persons who have had a primary TB infection previously for more than five years represents an increasing risk to develop active TB soon after reinfection. A study from South Africa [21] has demonstrated that the rate of reinfection by TB after successful treatment could be higher than the rate of new TB infections. In this study the reinfection rate after successful treatment was estimated at 2.2 per 100 person-years, which was approximately seven times the crude incidence rate (313 per 100 000 population per year) and approximately four times the age-adjusted incidence rate of new TB (515 per 100 000 population per year). So, ignoring exogenous reinfection when modeling TB spread in high-incidence and high-prevalence community setting such as semiclosed communities has been seen to be inappropriate. (Henao-Tamayo et al. in [22] recently published a mouse model of TB reinfection that could help to explain immunological aspects of reinfection risk in high-incidence areas.)

We will use an SEIR standard compartmental model; see for example the works by Blower et al. [23] and more recently by Liao et al. [24] with some modifications explained bellow that turn out to be quite useful in the study of the particularities of TB spread at this type of communities. This model assumes that the population in the community is homogeneous that it does not consider the heterogeneities in the social structure between community members, and it is based on the so-called mass action or fully mixing approximation. This means that individuals with whom a susceptible individual has contact are chosen at random from the whole community. It is also assumed that all individuals have approximately the same number of contacts in the same time and that all contacts transmit the disease with the same probability.

The model we use in this work takes into account the following relevant facts in the context of semiclosed communities.

(1)The overcrowding in the community can increase (compared to what occurs in the population at large) the likelihood of exogenous reinfection due to repeated contacts with active infected individuals. That is, besides primary infection the model considers the possible reinfection of individuals with LTB (individuals who are assumed to be asymptomatic and noninfectious but capable of progressing to active TB) and recovered individuals (individuals who have been treated for TB in the past and been declared cured). If latently infected or recovered individuals remain in the community, they could be infected again.(2)At present, it is not completely clear whether in all cases previous infections with Mycobacterium TB with or without subsequent recovery offer some protection that could be translated into a reduced susceptibility to reinfection [5, 21, 22, 25]. So, we will be open at exploring different situations with regard to this fact in the model.(3)Poor nutrition, immunodepression, and other diseases increase the likelihood of accelerated progression to active TB.

We will see that considering exogenous reinfection to describe TB spread produces a richer and more complex dynamics than the one observed in previous models (see e.g., [23, 25, 26]). In particular, unlike the model published by Feng et al. in [26], which uses a single parameter for exogenous reinfection, our model uses two parameters related to two possible reinfections (reinfection of latently infected and reinfection of recovered individuals).

2. Basic Epidemiology of TB Sources and Probability of Infection in Semiclosed Communities

The risk of infection with Mycobacterium tuberculosis, the bacterium causing TB, depends mainly on two factors: first, significant exposure to a source of infection and second, the probability of getting infection if there is exposure.

TB is mostly transmitted through the air; tubercle bacilli, that depends on host and agent factors, is distributed in tiny liquid droplets that are produced when someone with clinical or active TB coughs, sneezes, spits, or speaks, allowing infected individual to infect others. In closed places the bacteria are expelled into a finite volume of air unless there is ventilation, see [27]. In these conditions they may remain viable and suspended in the air for a prolonged period of time. But, the number of bacilli excreted by most persons with active pulmonary TB is relatively small [16], so the probability of TB transmission per contact, per unit of time is in general quite low. The risk of infection is very small during a single encounter with an infectious individual [28]. However, the probability of TB transmission can be enhanced by systematic and long exposure of susceptible individuals to particular infectious individuals.

The risk of TB transmission is particularly high in settings with poorly ventilated areas (places with reduced air volume per occupant, with ventilation systems which recirculate the air, or with poorly filtered air exchanges) and/or closed areas in which people are in close and frequent contact. Closed regime prisons are examples of these high-risk areas. In effect, the occurrence of TB in prisons for example is usually reported to be much higher than the average levels reported for the corresponding general population [6].

Although most exposed individuals develop an effective immune response to the initial infection [17], there is another factor that raises the chances of TB contagion, the fact that TB is an opportunistic disease. Indeed, infected individuals with weakened immune systems are at significant risk of developing clinical TB disease (active TB). High TB prevalence is therefore observed in individuals with HIV infection, poor nutritional status, alcoholism, drug abuse, concurrence of other pathology, and psychological stress decrease immune response levels. These conditions occur frequently in imprisoned peoples.

TB is usually described as a slow disease because of its long and variable period of latency and because of its short and relatively narrow infectious period distribution. Long periods of latency (inactive TB or latent TB or LTB) imply that new cases of infection are not clinically noticeable and therefore remain unobserved for a period of time. Immune response of susceptible individuals can restrict proliferation of the bacilli leading to what seems to be long-lasting partial immunity against reinfection or a response capable of stopping the progression from LTB to active TB.

Exposed individuals may remain in the latent stage for long and variable periods of time. In fact, it often happens that the host dies without ever developing active TB. The progression from latent to active TB is uncommon in the population at large. It is estimated that only about 5 to 10 percent of LTB individuals develop clinical or active TB [16], but due to the above described extreme conditions at semiclosed communities such as prisons, persons lived in these communities may be at risk of rapid progression from LTB to active TB following recent infection or reactivation of latent infection, or reinfection, see [6].

Some additional known epidemiological facts to be considered for TB disease are the following.

(1)Most of the secondary infections generated by an infected individual do take place within the first months following TB activation [29].(2)In the work by Styblo [16] it was noted that nearly 60 percent of the new cases arose during the first year following infection, while the cumulative number of cases generated over the first five years after infection accounted for nearly 95 percent of the total observed cases. People ill with TB can infect up to 10–15 other people through close contact over the course of a year [30].(3)Case fatality among untreated pulmonary TB cases is around 66.6 percent [30].(4)Recovered individuals, naturally or from treatment, may develop active TB again, a phenomenon known as TB relapse. (Recurrent cases (formerly relapse cases) have been treated for TB in the past and been declared successfully treated (cured/treatment completed) at the end of their treatment regimen. Recurrent cases include relapses due to the same Mycobacterium tuberculosis strain as for the previous episode as well as new episodes of TB due to reinfection.) (5)Individuals with LTB may progress to active TB due to reexposure and reinfection. The extent to which latent tuberculosis infection could reduce the risk of progressive disease following reinfection is not known [31].

3. A Compartmental Model for the TB Spread

In order to describe key relevant dynamics in the study of the TB spread in semiclosed communities, we consider five compartments SEIR model represented in Figure 1.

Figure 1: Flow chart of TB compartmental model.

The compartments are uninfected individuals (susceptible), the class; the latent class , that is, individuals who are assumed to be asymptomatic and noninfectious but capable of progressing to the clinical disease or active TB; the infectious class is subdivided into two subclasses: (a) infected and infectious individuals and (b) infected and noninfectious individuals ; and the class of recovered by treatment, self cure, or quarantine.

Every individual in the , , and classes is considered infected. There are five possible routes toward TB infection according to this model: primary infection after a contact with infected and infectious individuals (fast TB), endogenous reactivation after a period of latency (slow TB), relapse by natural causes after a cure, exogenous reinfection of latently infected, and exogenous reinfection of recovered individuals.

The and are probability of developing infectious TB if one develops fast and slow TB, respectively, is the relapse rate to active TB. Uninfected individuals are recruited at the rate , and is the natural mortality rate. Individuals with TB experience a death rate due to TB infection.

After infection, a fraction of individuals progresses to disease relatively soon (i.e., within the first two years) after infection; the remaining fraction of infected individuals become latently infected. Of newly infected individuals who thus progress quickly to disease, represents the fraction that develops infectious disease, and represents the fraction that develops noninfectious disease.

The class, latently infected individuals, does not shed bacilli and is not infective to others. In some latently infected individuals, the infection remains latent and it may persist for life. But, in a small minority of latently infected individuals, reactivation of the latent infection leads to the disease. The coefficients and denote the treatment rates for infected and infectious individuals class and infected and noninfectious individuals class, respectively. The model does not consider unsuccessful treatments.

The parameter is the primary TB transmission rate; this parameter summarizes socioeconomic and environmental factors that affect primary TB transmission. We assume that transmission rates are determined by broad demographic and social contexts, as well as by characteristics of both the transmitter and recipient (i.e., the number, viability, and virulence of the organisms within sputum droplet nuclei, immune status of the recipient, etc.)

TB transmission rate in case of reinfection might be different than the transmission rate of primary infection. The quantities that take into account these differences in case of reinfection of latently infected individuals and reinfection of recovered individuals are given by the dimensionless parameters and , respectively. The parameter is the proportion in TB transmission due to exogenous reinfection of latently infected individuals, and is the proportion in TB transmission due to exogenous reinfection of recovered individuals. Thus, is the exogenous reinfection rate of latently infected, and is the exogenous reinfection rate of recovered individuals.

A conservative point of view will consider that biologically plausible values for the reinfection parameters and are given within the intervals , . In this case, the parameters and can be interpreted as factors reducing the risk of reinfection of an individual who has previously been infected and has acquired some degree of protective immunity. However, studies on genetic predisposition [22] or in communities with cases as those reported in [21] have gathered some evidence that in certain situations there may be some increased susceptibility to reinfection. Therefore, we are willing to explore in the next sections other mathematical possibilities where the reinfection parameters can take even less usual values and .

However, recurrent TB due to endogenous reactivation (relapse) and exogenous reinfection could be clinically indistinguishable [32]; they are independent events. For this reason, beside primary infection we will include in the model the possibility of endogenous reactivation and exogenous reinfection as different way toward infection. So, we have the following.

(1)TB due to the endogenous reactivation of primary infection (exacerbation of an old infection) is considered in the model by the terms and .(2)TB due to reactivation of primary infection induced by exogenous reinfection is considered by the terms and .(3)Recurrent TB due to exogenous reinfection after a cure or treatment is described by the term .

The parameters of the model, its descriptions, and its units are given in Table 1.

Table 1: Parameters of the model, its descriptions, and its units.

Parameter
Description
Unit

β
Transmission rate
1/year

Recruitment rate
1/year

Natural cure rate
1/year

Progression rate from latent TB to active TB
1/year

Natural mortality rate
1/year

Mortality rate or fatality rate due to TB
1/year

Relapse rate
1/year

Probability to develop TB (slow case)
—

Probability to develop TB (fast case)
—

Proportion of new infections that produce active TB
—

Exogenous reinfection rate of latent
1/year

Exogenous reinfection rate of recovered
1/year

Treatment rates for
1/year

Treatment rates for
1/year

All these considerations give us the following system of equations:

Adding all the equations in (1) together, we have

where represents the total number of the population, and the region

is positively invariant of system (1).

It is a common practice in epidemic research to introduce the basic reproduction number , defined as the average number of secondary infections produced by an infected individual in a completely susceptible population, as the measure for the epidemic thresholds, if an epidemic will arise.

We have calculated for this model using the Next Generation Method [35] and it is given by

where

3.1. Steady-State Solutions

In order to find steady-state solutions for (1) we have to solve the following system of equations:

Solving system (6) with respect to we have the following equation:

The coefficients of (7) are all expressed as functions of the parameters listed in Table 1. However, these expressions are too long to be written here. See Appendix A for explicit forms of the coefficients.

3.1.1. Disease-Free Equilibrium

For we get the disease-free steady-state solution:

Lemma 1. The disease-free steady-state solution is locally asymptotic stable for .

Proof. The Jacobian matrix of system (6) evaluated in is

The characteristic equation for have the form

Given the polynomial

in the special case when , roots of the polynomial have negative real part and if

(i), the 4th root, or largest eigenvalue, is zero,(ii), all eigenvalues are negative,(iii), the largest eigenvalue has positive real part.

Thus, the stability of disease-free steady-state solution is determined solely by the sign of the constant term of the polynomial [36].

The coefficients are all decreasing functions of and they are linear functions with respect to the parameter , so they all take the general form with and . We can define . It is easy to see that for we have .

For example,

By straightforward calculations and reminding that the coefficients are decreasing functions of , we found that

This way, and if we take , all the coefficients are positive. But, for , we can see that the condition is fulfilled. Indeed, the constant term of the polynomial can be written as

Using this form for the coefficient we can see that if , then so .

Remark 2. For we have , and the disease-free steady-state solution is unstable. Indeed, if , and are the roots of polynomial , we have , and it is impossible that the four roots have negative real part.

3.1.2. Endemic Equilibria

From (7) we get that nontrivial solutions are possible if

The explicit expressions for the coefficients and are

where parameters , , , and are defined as in (5).

The coefficients and can be written in the following general form:

where

The coefficients and are all positive and depend on the parameters given in Table 1. See Appendix A for the explicit form of these coefficients.

Changes in the signs of the coefficient and as function of transmission rate can be explained using the above defined functions and , respectively. The functions and both are linear and decreasing functions of .

Consider the polynomial function

From (17) we can see that for the coefficient is negative, so we have . On the other hand, because the coefficient is always positive, there must be a value such that, for , it holds that . Since function is continuous, this implies the existence of solution for the equation .

To determine how many possible endemic states arise, we consider the derivative , and then we analyse the following cases.

(1)If , for all , then is monotonically increasing function and we have a unique solution, that is, a unique endemic equilibrium.(2)If , we have solutions of the equation given by

and for all and . So, we need to consider the positions of the roots and in the real line. We have the following possible cases.

(i)If , then for both cases and , we have , and for all . Given that , this implies the existence of a unique endemic equilibrium.(ii)If and , then both roots and are negative and for all .(iii)If and , then both roots and are positive and we have the possibility of multiple endemic equilibria. This is a necessary condition, but not sufficient. It must be fulfilled also that .

Let be the value of such that and the value of such that . Moreover, let be the value for which the basic reproduction number is equal to one (the value of such that coefficient becomes zero).

Lemma 3. If the condition is met, then system (1) has a unique endemic equilibrium for all (Table 3).

Proof. Using similar arguments to those used in the proof of Lemma 1,  we have, given the condition , that for all values of such that , all polynomial coefficients are positive; therefore, all solutions of the polynomial are negative and there is no endemic equilibrium (positive epidemiologically meaningful solution).

For the coefficients and are both positive, while the coefficient is negative; therefore, appears only one positive solution of the polynomial (the greatest one), so we have a unique endemic equilibrium.

For , the coefficient is negative and is positive. According to the cases studied above we have in this situation a unique endemic equilibrium.

Finally, for the coefficients and are both negative, and according to the study of cases given above we also have a unique positive solution or endemic equilibrium.

Let us first consider biologically plausible values for the reinfection parameters and , that is, values within the intervals , . This means that the likelihood of both variants of reinfections is no greater than the likelihood of primary TB. So, we are considering here partial immunity after a primary TB infection.

Lemma 4. For biologically plausible values system (1) fulfils the condition .

Proof. Using straightforward but cumbersome calculations (we use a symbolic software for this task), we were able to prove that if we consider all parameters positive (as it is the case) and taking into account biologically plausible values , then and and it is easy to see that these inequalities are equivalent to .

We have proven that the condition implies that the system can only realize two epidemiologically meaningful (nonnegative) equilibrium states. Indeed, if we consider the disease transmission rate as a bifurcation parameter for (1), then we can see that the system experiences a transcritical bifurcation at , that is, when (see Figure 2). If the condition is met, the system has a single steady-state solution, corresponding to zero prevalence and elimination of the TB epidemic for , that is, , and two equilibrium states corresponding to endemic TB and zero prevalence when , that is, . Moreover, according to Lemma 4 this condition is fulfilled in the biologically plausible domain for exogenous reinfection parameters . This case is summarized in Table 2.

Table 2: Qualitative behaviour for system (1) as a function of the disease transmission rate , when the condition is fulfilled.

Interval
Coefficients
Type of equilibrium

, , ,
Disease-free equilibrium

, , ,
Unique endemic equilibrium

, , ,
Unique endemic equilibrium

, , ,
Unique endemic equilibrium

Table 3: Qualitative behaviour for system (1) as function of the disease transmission rate , when the condition is fulfilled. Here, is the discriminant of the cubic polynomial (20).

Interval
Coefficients
Type of equilibrium

, , ,
Disease-free equilibrium

, , ,
Two equilibria () or none ()

, , ,
Two equilibria () or none ()

, , ,
Unique endemic equilibrium

Figure 2: Bifurcation diagram (solution of polynomial (20) versus ) for the condition . is the bifurcation value. The blue branch in the graph is a stable endemic equilibrium which appears for .

From Table 2 we can see that although the signs of the polynomial coefficients may change, other new biologically meaningful solutions (nonnegative solutions) do not arise in this case. The system can only display the presence of two equilibrium states: disease-free or a unique endemic equilibrium.

The basic reproduction number in this case explains well the appearance of the transcritical bifurcation, that is, when a unique endemic state arises and the disease-free equilibrium becomes unstable (see blue line in Figure 2).

However, the change in signs of the polynomial coefficients modifies the qualitative type of the equilibria. This fact is shown in Figures 5 and 7 illustrating the existence of focus or node type steady-sate solutions. These different types of equilibria as we will see in the next section cannot be explained using solely the reproduction number .

In the next section we will explore numerically the parametric space of system (1), looking for different qualitative dynamics of TB epidemics. We will discuss in more detail how dynamics depends on the parameters given in Table 1, especially on the transmission rate , which will be used as bifurcation parameter for the model.

Let us consider here briefly two examples of parametric regimes for the model in order to illustrate the possibility to encounter a more complex dynamics, which cannot be solely explained by changes in the value of the basic reproduction number .

Example I. Suppose , this implies that and ; therefore, we have the equation:

It is easy to see that besides zero solution, if , and , (22) has two positive solutions and . So, we have in this case three nonnegative equilibria for the system.

The condition for means , and this in turn implies that . On the other hand, the condition implies and therefore .

Gathering both inequalities we can conclude that if , then the system has the possibility of multiple equilibria.

Since the coefficients and are both continuous functions of , we can always find a neighbourhood of , such that the signs of these coefficients are preserved. Although in this case we do not have the solution , we eventually could still have two positive solutions and consequently, multiple equilibrium states; see the green line in Figure 3.

Figure 3: Polynomial for different values of with the condition . The graphs were obtained for values of and . The dashed black line indicates the case . The figure shows the existence of multiple equilibria.

Example II. Suppose we take numerical values for the parameters in Table 1 such that the condition is fulfilled.

If , then all coefficients of the polynomial (20) are positive and there is not nonnegative solutions. In this case, the system has only a disease-free equilibrium.

For and the signs of the coefficients of the polynomial are , , , and , , , , respectively.

In both cases the polynomial has two possibilities:

(a)three real solutions: one negative and two positive solutions for ,(b)one negative and two complex conjugate solutions for .

Here is the discriminant for the polynomial (20).

In the (a) case we have the possibility of multiple endemic states for system (1). This case is illustrated in numerical simulations in the next section by Figures 8 and 9.

We should note that the value is not a bifurcation value for the parameter .

If , then , , , and . In this case we have

The discriminant is a continuous function of , for this reason this sign will be preserved in a neighbourhood of .

We should be able to find a bifurcation value solving numerically the equation

where can be bounded by the interval (see Figure 4).

Figure 4: Bifurcation diagram for the condition . is the bifurcation value. The blue branch in the graph is a stable endemic equilibrium which appears even for < 1.

Figure 5: Numerical simulation for , , , and . The system goes toward a focus type stable stationary equilibrium.

4. Numerical Simulations

In this section we will show some numerical simulations with the compartmental model (1). This model has fourteen parameters that have been gathered in Table 1. In order to make the numerical exploration of the model more manageable, we will adopt the following strategy.

(i)First, instead of fourteen parameters we will reduce the parametric space using four independent parameters , , , and . The parameters , , and are the transmission rate of primary infection, exogenous reinfection rate of latently infected, and exogenous reinfection rate of recovered individuals, respectively. is the value of such that basic reproduction number is equal to one (or the value of such that coefficient in the polynomial (20) becomes zero). On the other hand, depends on parameters given in the list . This means that if we keep all the parameters fixed in the list , then is also fixed. In simulations we will use instead of using basic reproduction number .(ii)Second, we will fix parameters in the list according to the values reported in the literature. In Table 4 are shown numerical values that will be used in some of the simulations, besides the corresponding references from where these values were taken. Mostly, these numerical values are related to data obtained from the population at large, and in the next simulations we will change some of them for considering the conditions of extremely high incidence/prevalence of TB in semiclosed communities. In any case, these changes will be clearly indicated in the text.(iii)Third, for any pairs of values and we can compute and , that is, the values of such that and , respectively, in the polynomial (20). So, we have that the exploration of parametric space is reduced at this point to the study of the parameters , , , and . According to the chosen values for , , and , we have six possible orderings for the parameters , , and (see Appendix B).The dynamic behavior of system (1) will depend of these orderings. In particular, from Table 5, it is easy to see that if then the system has a unique equilibrium point, which represents a disease-free state, and if , then the system has a unique endemic equilibrium, besides an unstable disease-free equilibrium.(iv)Fourth and finally, we will change the value of , which is considered a bifurcation parameter for system (1), taking into account the previous mentioned ordering to find different qualitative dynamics.

Table 4: Numerical values for the parameters in the list . Some of the given numerical values for the model parameters are mainly related to the spread of TB in the population at large and are basically taken as reference. Other values assuming for the parameters, different than those given in this table will be clearly indicated in the text.

Parameter 
Description
Value

Recruitment rate
200 (assumed)

Natural cure rate
0.058 [23, 33, 34]

Progression rate from latent TB to active TB
0.0256 [33, 34]

Natural mortality rate
0.0222 [2]

Mortality rate due to TB
0.139 [2, 33]

Relapse rate
0.005 [2, 33, 34]

Probability to develop TB (slow case)
0.85 [2, 33]

Probability to develop TB (fast case)
0.70 [2, 33]

Proportion of new infections that produce active TB
0.05 [2, 33, 34]

Treatment rates for
0.50 (assumed)

Treatment rates for
0.20 (assumed)

Table 5: Different possible orderings for , , and . In every case , is the cubic discriminant of the equation , is the discriminant of the quadratic equation , where is the polynomial (20).

Interval
Coefficients
Equilibria

, ,
Disease-free equilibrium

, ,
Unique endemic equilibrium

, ,
Unique endemic equilibrium

, ,
Unique endemic equilibrium

, ,
Disease-free equilibrium

, ,
Two equilibria if or ; none if or

, ,
One equilibrium for or

, ,
Unique endemic equilibrium

, ,
Disease-free equilibrium

, ,
Two equilibria if or ; none if or

, ,
Two equilibria () or none ()

, ,
Unique endemic equilibrium

, ,
Disease-free equilibrium

, ,
Unique endemic equilibrium

, ,
One equilibrium (), three equilibria ()

, ,
Unique endemic equilibrium

, ,
Disease-free equilibrium

, ,
Two equilibria () or none ()

, ,
One equilibrium (), three equilibria ()

, ,
Unique endemic equilibrium

, ,
Disease-free equilibrium

, ,
Two equilibria () or none ()

, ,
Two equilibria () or none ()

, ,
Unique endemic equilibrium

It is especially interesting to explore the consequences of modifications in the values of the reinfection parameters without changing the values in the list , because in this case the threshold remains unchanged. Thus, we can study in a better way the influence of the reinfection in the dynamics of the TB spread.

The values given for the reinfection parameters and in the next simulations could be extreme, trying to capture this way the special conditions of high burden semiclosed communities.

Example I (Case ,  , ). Let us consider here the case when the condition is met. We know from the previous section that this condition is met under biologically plausible values .

According to Lemmas 3 and 4, in this case the behaviour of the system is characterized by the evolution towards disease-free equilibrium if and the existence of a unique endemic equilibrium for . Changes in the parameters of the list alter the numerical value of the threshold but do not change this behaviour.

First, we consider the following numerical values for these parameters: , , and . We also fix the list of parameters according to the numerical values given in Table 4.

The basic reproduction number for these numerical values gives . The initial conditions considered were

We also have the following values:

These values clearly meet the condition , and according to Lemma 3 the system must have in this case a unique endemic equilibrium for all .

Figure 5 shows that under the above described situation, the system will converge to an endemic equilibrium given by the focus type stationary steady solution:

By straightforward calculations we can show that this focus is stable, and no matter what initial conditions are taken for the system, the solutions always evolve to this endemic state.

Figure 6 shows the trajectories of the system for multiple initial conditions in a three-dimensional phase space in which the horizontal axes are susceptible and recovered individuals, while the vertical axis is the prevalence .

Figure 6: Phase space representation of the evolution of the system toward a stable focus type equilibrium. In this representation were used multiple initial conditions and the following values: , , , and .

Figure 7: Numerical simulation for , , , and . In this case the system converges to a stable node type equilibrium.

Figure 8: Numerical simulation for , , and . The system can evolve to two different equilibria (red lines) or (dark green lines) according to different initial conditions.

Figure 9: Numerical simulation for , , and . Phase space representation of the system with multiple equilibrium points.

Example II (Case ,  , ).  For our next numerical simulation we consider the following values for the used parameters: , , , and as before the list of parameters is fixed according to Table 4.

The basic reproduction number for these parameters as before gives the same value . The used initial conditions were

We also have the following values:

These values meet the condition , and as in the previous simulation the system evolves toward a unique endemic equilibrium, but this time the dynamical properties of the equilibrium have changed.

In fact, Figure 7 shows the evolution of the system toward a stable node type endemic equilibrium:

In our model, considering biologically plausible domain for exogenous reinfection parameters , the condition is fulfilled. Under this condition we have a unique endemic equilibrium for . The emergence by a transcritical bifurcation of this endemic state is properly explained by the basic reproduction number . However, changes in the reinfection parameters , can modify the qualitative nature of the dynamics of the disease, in addition to changing the numbers of individuals in the different compartments of the model in the endemic equilibrium state, without having any change in the value of the basic reproduction number , which in this case fails to describe these variations in the dynamics of the disease.

Example III (Case , , ).  There is now evidence that rates of secondary tuberculosis in high endemic communities (for example semiclosed communities), in patients with LTB or/and already treated for primary disease, are actually higher than in people presenting with primary infection [21, 22]. Taking this into consideration we consider now the following numerical values for the parameters: , , . In this case the basic reproduction number takes the value . Additionally we have

For these parameter we have that the condition is fulfilled and the system has the possibility of multiple equilibria. In fact, we have in this case the following stationary points :

is a stable disease-free equilibrium point (stable node), is a stable endemic equilibrium (stable focus), and is an unstable equilibrium point (saddle point).

Figure 8 shows the convergence to or to according to with different initial conditions.

In Figure 9 is shown another representation (phase space) of the evolution of the system toward or to according to different initial conditions. The representation is a three-dimensional phase space in which the horizontal axes are susceptible and recovered individuals, while the vertical axis is the prevalence .

For the previously numerical values, the system experiences a backward bifurcation [37] at the value with . For , the system possesses two stable equilibrium points and one unstable (see Figure 4).

Example IV (Case , , ). Consider now a more extreme situation with , , and (the other parameters kept the same values given in Table 4). In this case the condition is fulfilled.

This example is shown in order to illustrate more complex and rich dynamics that might admit system (1), which is mathematically possible and could in principle be a model case for an extreme hypothetical situation in a semiclosed high burden community. For these parameters we have

which clearly satisfy the condition . Therefore, as was explained in the previous section, the system has the possibility of multiple equilibria.

In fact, for the bifurcation value of the disease transmission rate, which satisfies the condition , the system acquires two positive equilibria, apart from the disease-free equilibrium.

When appear three positive equilibrium points and the disease-free equillibrium becomes unstable. For with the system admits a unique and stable endemic equilibrium (see Figure 10).

Figure 10: Bifurcation diagram (solution of polynomial (20) versus ) for the condition . The system experiences multiple bifurcations at , , and .

We take now the value , which satisfies the condition .

With these numerical values the basic reproduction number is , and therefore, the disease-free equilibrium is stable.

We have in this case the following stationary points :

is the stable disease-free equillibrium point (stable node), is an unstable equilibrium point (saddle point), and is a stable endemic equilibrium (stable focus). Figure 11 shows the convergence to or to according to the initial condition.

Figure 11: Numerical simulation for , , and . The system can evolve to two different equilibria or according to the initial condition.

In Figure 12 is shown another representation (phase space) of the evolution of the system toward or to according to the initial conditions.

Figure 12: Numerical simulation for , , and . Phase space representation of the system with multiple equilibrium points.

Let us take now the value , which satisfies the condition . In this case, the basic reproduction number has the value . We still have that the condition is fulfilled and the system in this case has four equilibrium points :

is the unstable disease-free equillibrium point (saddle point ), is a stable endemic equilibrium point (node), is an unstable equilibrium (saddle point), and is a stable endemic equilibrium point (focus).

Figure 13 shows the phase space representation of this case.

Figure 13: Numerical simulation for , , and . The system can evolve to two different equilibria (stable node) or (stable focus) according to the initial condition. and are unstable equilibria.

For further numerical analysis, we set all the parameters in the list according to the numerical values given in Table 4, leaving free the parameters , , and related to the primary transmission rate and reinfection rates of the disease.

We will explore the parametric space of system (1) and relate it to the signs of the coefficients of the polynomial (20).

In Figure 14, we consider values of such that . We can observe from this figure that as the primary transmission rate of the disease increases, and with it the basic reproduction number , the system under biological plausible condition, represented in the figure by the square , evolves such that initially (for lower values of ) coefficients and are both positive, then remains positive and becomes negative and finally both coefficients become negative.

(a)

(b)

(c)

(d)

Figure 14: Signs of coefficients and as functions of exogenous reinfection rate of latent and exogenous reinfection rate of recovered for . The parameter has the values: (a) , (b) , (c) , and (d) .

This change in the coefficients signs as the transmission rate increases agrees with the results summarized in Table 2 when the condition is fulfilled.

Next, in order to explore another mathematical possibilities we will modify some numerical values for the parameters in the list in a more extreme manner, taking a hypothetical regime with , , , , , , , , , .

In Figure 15 besides signs of and we consider also the signs of the discriminant of the quadratic equation , where is the polynomial (20).

(a)

(b)

(c)

(d)

Figure 15: Signs of coefficients , , and discriminant as functions of exogenous reinfection rate of latent and exogenous reinfection rate of recovered for . The parameter has the values: (a) , (b) , (c) , and (d) .

From Figure 15, in particular we can see that the domain with , , , and represented in red allows the possibility of multiple endemic equilibria for system (1). However, despite the extreme numerical values for the parameters taking in , this domain is still far from the domain of biologically plausible values for and , represented in the figures by the square .

As the transmission rate of the disease increases, and with it the basic reproduction number , this red domain of the parametric space becomes increasingly smaller until finally it disappears. In fact, the red domain is only significant when the basic reproduction number is near one.

In Figure 16, we show some numerical simulation where basic reproduction number is less than one. The red domain indicates the possibility of multiple endemic equilibria for the system even for . We can see that this domain is far from the domain of biologically plausible values for and represented in the figure by the square . As the transmission rate of the disease decreases, and with it the number , this parameter domain moves away from the square. In all represented cases and inside the square .

(a)

(b)

(c)

(d)

Figure 16: Signs of coefficients , , and discriminant as functions of exogenous reinfection rate of latent and exogenous reinfection rate of recovered for . The parameter has the values: (a) , (b) , (c) , and (d) .

5. Discussion and Conclusions

In order to consider high incidence and prevalence of TB and LTB in semiclosed communities, we have used in this work a compartmental SEIR model with five possible pathways to TB disease. The extra nonlinear terms considered in the model lead to a more complex mathematical treatment in comparison with previously used models (see e.g., [23, 25–29, 32, 35, 36]). But the special form of some coefficients obtained from the analysis of standard SEIR models with constant transmission rate allowed us to move forward with some analytical results that were confirmed later by numerical simulations.

In this paper we follow a new methodology for numerical exploration of this kind of models, which allowed us to handle the high dimensionality of its parameter spaces and thus study its different dynamic behaviors. We found that the transmission rate can be used as a bifurcation parameter and that the system undergoes qualitative changes in the dynamics when this parameter varies. In this context the analysis of the parametric space is reduced to the study of variations of the parameters , , , and . We divided the parametric space into six possible arrangements for these parameters, which in turn determine all possible different qualitative behaviours of the system dynamics.

From model (1) we can see that reinfection requires latently infected, recovered, and actively infectious individuals. The basic reproduction number has to do solely with infections produced when an infectious individual is introduced into an uninfected population, but reinfection does not alter the dynamics of an uninfected population. So, the number does not completely describe the dynamics of the model when the reinfection is incorporated, as was noted before by Feng et al. in [26]. Unlike the model published by these authors, which uses a single parameter for exogenous reinfection, in our model we use two parameters related to two possible pathways of reinfection (reinfection of latently infected and reinfection of recovered individuals).This is a reason why our model shows a more complex and richer dynamics.

We have showed through theoretical analysis in Section 3 and numerical simulations in Section 4 that if we accept as valid the plausible assumption that exposure to mycobacterium induces an immune response, which is partially protective against reinfection, then the system for semiclosed communities (1) reproduces well, common observed trends in TB epidemiology that are similar to what happens in population at large, which is basically that, for , there is only one disease-free status, while for , there exists a unique endemic state with nonzero prevalence. For occurs a transcritical bifurcation from which emerges an endemic stable state.

Moreover, according to Lemmas 3 and 4, any values of reinfection parameters in this parametric regime: would lead to the same qualitative dynamics and will not affect this already classical behavior in SEIR models. In this case only one of the aforementioned arrangements emerges as valid under this biologically plausible condition.

Since the two parameters related to exogenous reinfection of latently infected and recovered individuals do not affect the value of the number , even under the plausible assumption of partial immunity, variation of reinfection parameters can make that for the same value of the number , the quality of dynamics and the number of affected by disease individuals (incidence and prevalence) drastically change. For example, Figures 5 and 7 show two types of dynamics, that is, convergences to different stationary points, a focus and a node for the same basic reproduction number . Some evidence of this variability in tuberculosis epidemiology due to dynamic balance between primary infection and reinfection has been presented in several works (see e.g., [26, 38]).

Taking less plausible assumption, but already evidenced in several works [5, 21, 22, 26], of an increased susceptibility to reinfection over primary infection in some cases leads us to a further study of model (1). For and system (1) experiences a rich and complex dynamics with successive and different kind of bifurcations as the transmission rate changes. These cases incorporate possible multiple endemic states, regardless of whether the values for the basic reproduction number were less than or greater than 1. So, these behaviors cannot be explained using only this number. It is in this context that the use of the disease transmission rate as bifurcation parameter instead of acquires real usefulness.

Some important implications of the simulations with and lie in the fact that many of the measures taken to stop and control an epidemics are designed to reduce the value of the basic reproduction number such that disease-free status for    is achieved. However, in this parametric regime, reinfection might cause the system to fall into a state unable to eliminate endemic disease, although it fulfills that . Thus, semiclosed communities with this kind of regime will become in genuine high transmission pockets of TB inserted in the general population [4]. Indeed, semiclosed communities such as prisons might become in a reservoir for disease transmission to the population at large and should be a source of public concern [4, 6, 7].

The theoretical approach and numerical simulations presented in this paper for the study of the impact of reinfection on TB dynamics in semiclosed communities could have important implications at multiple levels, including vaccine design, control program design, epidemiology of tuberculosis in regions where the risk of reexposure is high, and for systems-based computer models which to date assume that primary infection will confer at least some degree of (stable) memory immunity to a secondary infection, but that in fact also have to consider less plausible assumptions about an increased susceptibility to reinfection.

Appendices

A. Explicit Form of Coefficients

Introducing the notations:

We have

The coefficients and can be written in the following general form:

where

B. Summary for Possible Orderings for the Different Transmission Rate Parameters

Table 5 shows different orderings for the parameters , , and and their implications concerning the possible equilibrium or stationary points of system (1). This table is complemented with the analysis given in Section 3.1.2.

Acknowledgments

This work was supported by Concurso Nacional de Proyectos de Investigación en Salud, FONIS project/Code SA11I2073. The authors wish to thank Claudia González, Macarena Hirmas, Patricia González, Juan Carlos Hormazabal, Iris Delgado, and the students Matias Meyer and Ilani Kauffman for their valuable collaboration to the research.

References

C. Castillo-Chavez and B. Song, “Dynamical models of tuberculosis and applications,” Mathematical Biosciences and Engineering, vol. 1, pp. 361–404, 2004.

C. Ozcaglar, A. Shabbeer, S. L. Vandenberg, B. Yener, and K. P. Bennett, “Epidemiological models of Mycobacterium tuberculosis complex infections,” Mathematical Biosciences, vol. 236, pp. 77–96, 2012.

J. Raffalli, K. A. Sepkowitz, and D. Armstrong, “Community-based outbreaks of tuberculosis,” Archives of Internal Medicine, vol. 156, no. 10, pp. 1053–1060, 1996.

S. Basu, D. Stuckle, and M. McKee, “Addressing institutional amplifiers in the dynamics and control of tuberculosis epidemics,” The American Journal of Tropical Medicine and Hygiene, vol. 84, no. 1, pp. 30–37, 2011.

C. Y. Chiang and L. W. Riley, “Exogenous reinfection in tuberculosis,” The Lancet Infectious Diseases, vol. 5, no. 10, pp. 629–636, 2005.

I. Baussano, B. G. Williams, P. Nunn, M. Beggiato, U. Fedeli, and F. Scano, “Tuberculosis incidence in prisons: a systematic review,” PLoS Medicine, vol. 7, no. 12, Article ID e1000381, 2010.

FONIS Project, Code SA11I2073, Ongoing research in: determinants of TB transmission in the prison population and its impact as a reservoir for the general population of Chile. CEPS, Facultad de Medicina, Clinica Alemana, Universidad del Desarrollo

Z. W. Jia, G. Y. Tang, Z. Jin et al., “Modeling the impact of immigration on the epidemiology of tuberculosis,” Theoretical Population Biology, vol. 73, no. 3, pp. 437–448, 2008.

Y. Zhou, K. Khan, Z. Feng, and J. Wu, “Projection of tuberculosis incidence with increasing immigration trends,” Journal of Theoretical Biology, vol. 254, no. 2, pp. 215–228, 2008.

W. W. Stead, “Pathogenesis of a first episode of chronic pulmonary tuberculosis in man: recrudescence of residuals of the primary infection or exogenous reinfection?” American Review of Respiratory Disease, vol. 95, no. 5, pp. 729–745, 1967.

R. Sahadevan, S. Narayanan, C. N. Paramasivan, R. Prabhakar, and P. R. Narayanan, “Restriction fragment length polymorphism typing of clinical isolates of Mycobacterium tuberculosis from patients with pulmonary tuberculosis in Madras, India, by use of direct-repeat probe,” Journal of Clinical Microbiology, vol. 33, no. 11, pp. 3037–3039, 1995.

L. K. Fitzpatrick, A. Okwera, R. Mugerwa, R. Ridzon, J. Ellner, and I. Onorato, “An investigation of suspected exogenous reinfection in tuberculosis patients in Kampala, Uganda,” International Journal of Tuberculosis and Lung Disease, vol. 6, no. 6, pp. 550–552, 2002.

A. Krüüner, L. Pehme, S. Ghebremichael, T. Koivula, S. E. Hoffner, and M. Mikelsaar, “Use of molecular techniques to distinguish between treatment failure and exogenous reinfection with Mycobacterium tuberculosis,” Clinical Infectious Diseases, vol. 35, no. 2, pp. 146–155, 2002.

P. D. van Helden, “Molecular epidemiology of TB: challenging dogmas and asking new questions,” IUBMB Life, vol. 53, no. 4-5, pp. 219–223, 2002.

P. E. M. Fine and P. M. Small, “Exogenous reinfection in tuberculosis,” The New England Journal of Medicine, vol. 341, no. 16, pp. 1226–1227, 1999.

K. Styblo, “Epidemiology of tuberculosis,” Selected Papers 24, Royal Netherlands Tuberculosis Association, The Hague, The Netherlands, 1991.

P. G. Smith and A. R. Moss, “Epidemiology of tuberculosis,” in Tuberculosis, Pathogenesis, Protection, and Control, B. R. Bloom, Ed., ASM Press, Washington, DC, USA, 1994.

J. A. Romeyn, “Exogenous reinfection in tuberculosis,” American Review of Respiratory Disease, vol. 101, no. 6, pp. 923–927, 1970.

A. van Rie, R. Warren, M. Richardson et al., “Exogenous reinfection as a cause of recurrent tuberculosis after curative treatment,” The New England Journal of Medicine, vol. 341, no. 16, pp. 1174–1179, 1999.

A. de Boer, R. van Soolingen, and M. Borgdorff, “Recurrent tuberculosis due to exogenous reinfection,” The New England Journal of Medicine, vol. 342, no. 14, pp. 1050–1051, 2000.

S. Verver, R. M. Warren, N. Beyers et al., “Rate of reinfection tuberculosis after successful treatment is higher than rate of new tuberculosis,” American Journal of Respiratory and Critical Care Medicine, vol. 171, no. 12, pp. 1430–1435, 2005.

M. Henao-Tamayo, A. Obregn-Henao, D. J. Ordway, S. Shang, C. G. Duncan, and I. M. Orme, “A mouse model of tuberculosis reinfection,” Tuberculosis, vol. 92, no. 3, pp. 211–217, 2012.

S. M. Blower, A. R. McLean, T. C. Porco et al., “The intrinsic transmission dynamics of tuberculosis epidemics,” Nature Medicine, vol. 1, no. 8, pp. 815–821, 1995.

C. M. Liao, Y. H. Cheng, Y. J. Lin et al., “A probabilistic transmission and population dynamic model to assess tuberculosis infection risk,” Risk Analysis, vol. 32, no. 8, pp. 1420–1432, 2012.

J. P. Aparicio and C. Castillo-Chavez, “Mathematical modelling of tuberculosis epidemics,” Mathematical Biosciences and Engineering, vol. 6, no. 2, pp. 209–237, 2009.

Z. Feng, C. Castillo-Chavez, and A. F. Capurro, “A model for tuberculosis with exogenous reinfection,” Theoretical Population Biology, vol. 57, no. 3, pp. 235–247, 2000.

E. A. Nardell, J. Keegan, S. A. Cheney, and S. C. Etkind, “Airborne Infection: theoretical limits of protection achievable by building ventilation,” American Review of Respiratory Disease, vol. 144, no. 2, pp. 302–306, 1991.

T. Baxter, “Low infectivity of tuberculosis,” The Lancet, vol. 342, no. 8867, article 371, 1993.

D. S. Barnes, The Making of a Social Disease: Tuberculosis in the Nineteenth-Century France, University of California Press, Berkeley, Calif, USA, 1995.

WHO TB fact sheets, 2012, http://www.who.int/mediacentre/factsheets/fs104/en/

J. R. Andrews, F. Noubary, R. P. Walensky, R. Cerda, E. Losina, and C. R. Horsburgh, “Risk of progression to active tuberculosis following reinfection with Mycobacterium tuberculosis,” Clinical Infectious Diseases, vol. 54, no. 6, pp. 784–791, 2012.

M. L. Lambert, E. Hasker, A. van Deun, D. Roberfroid, M. Boelaert, and P. van der Stuyft, “Recurrence in tuberculosis: relapse or reinfection?” The Lancet Infectious Diseases, vol. 3, no. 5, pp. 282–287, 2003.

M. A. Sanchez and S. M. Blower, “Uncertainty and sensitivity analysis of the basic reproductive rate: tuberculosis as an example,” American Journal of Epidemiology, vol. 145, no. 12, pp. 1127–1137, 1997.

J. Legrand, A. Sanchez, F. Le Pont, L. Camacho, and B. Larouze, “Modeling the impact of tuberculosis control strategies in highly endemic overcrowded prisons,” PLoS ONE, vol. 3, no. 5, Article ID e2100, 2008.

O. Diekmann, J. A. Heesterbeek, and J. A. Metz, “On the definition and the computation of the basic reproduction ratio R ₀ in models for infectious diseases in heterogeneous populations,” Journal of Mathematical Biology, vol. 35, pp. 503–522, 1990.

J. M. Heffernan, R. J. Smith, and L. M. Wahl, “Perspectives on the basic reproductive ratio,” Journal of the Royal Society Interface, vol. 2, pp. 281–293, 2005.

J. Dushoff, W. Huang, and C. Castillo-Chavez, “Backwards bifurcations and catastrophe in simple models of fatal diseases,” Journal of Mathematical Biology, vol. 36, no. 3, pp. 227–248, 1998.

M. G. M. Gomes, A. O. Franco, M. C. Gomes, and G. F. Medley, “The reinfection threshold promotes variability in tuberculosis epidemiology and vaccine efficacy,” Proceedings of the Royal Society B, vol. 271, no. 1539, pp. 617–623, 2004.