## Abstract

Transmissible spongiform encephalopathies (TSEs), such as kuru, are invariably fatal neurodegenerative conditions caused by a malformation of the prion protein. Heterozygosity of codon 129 of the prion protein gene has been associated with increased host resistance to TSEs, although the mechanism by which this resistance is achieved has not been determined. To evaluate the epidemiological mechanism of human resistance to kuru, we developed a model that combines the dynamics of kuru transmission and the population genetics of human resistance. We fitted our model to kuru data from the epidemic that occurred in Papua New Guinea over the last hundred years. To elucidate the epidemiological mechanism of human resistance, we estimated the incubation period and transmission rate of kuru for codon 129 heterozygotes and homozygotes using kuru incidence data and human genotype frequency data from 1957 to 2004. Our results indicate that human resistance arises from a combination of both a longer incubation period and reduced susceptibility to infection. This work provides evidence for balancing selection acting on a human population and the mechanistic basis for the heterozygote resistance to kuru.

## 1. Introduction

Prion diseases, or transmissible spongiform encephalopathies, are neurodegenerative disorders that include kuru and Creutzfeldt–Jakob disease (CJD) in humans, bovine spongiform encephalopathy in cattle, chronic wasting disease in deer and scrapie in sheep. Late-onset neurodegenerative conditions, such as Parkinson's and Alzheimer's diseases, may also have prion-like mechanisms [1]. Transmissible spongiform encephalopathies result from the accumulation [2–4] of post-translational conversion of a normal cellular protein (PrP^{C}) into an aggregated isomer (PrP^{Sc}) [5–7] and that leads to neuronal apoptosis [8–11].

Genetic susceptibility to human prion diseases is conferred by homozygosity of an amino acid located at codon 129 of the prion protein gene (*PRNP*). Codon 129 is polymorphic in human populations for methionine or valine [12–17]. Homozygosity of either allele has been associated with an earlier onset of kuru and a shorter clinical duration [12]. Furthermore, in the UK, where 51 per cent of the population has the heterozygote genotype frequency, codon 129 homozygotes constituted 21 of the 22 confirmed cases of sporadic CJD, 19 of the 23 suspected sporadic CJD cases [17], all 147 genotyped cases of variant CJD [18] and four of the seven iatrogenic CJD cases [13].

The kuru epidemic in Papua New Guinea that peaked in the 1950s most probably originated from a case of sporadic CJD, although this has not been confirmed [19]. The epidemic was propagated by ritualistic cannibalism of deceased tribe members [20]. In 1959, cannibalism was banned [19], and no individual born after 1959 has since been diagnosed with kuru [21]. However, individuals infected before 1959 continue to die from kuru after an extended incubation period [14,22,23]. The kuru epidemic in the Eastern Highlands of New Guinea has resulted in at least 1874 deaths, predominately in the Fore population, which had a size of 11 991 in 1958 [24,25].

Here, we develop a dynamic mathematical model of kuru transmission to first estimate key epidemiological parameters of kuru and then use those estimates to determine the mechanism of human resistance to kuru. This is the first use of transmission modelling to determine the incubation time of kuru. We also present the first estimate of the transmissibility of kuru, thereby allowing us to calculate the basic reproductive number of kuru. Furthermore, we evaluate the mechanism for genetic resistance to kuru, revealing a heterozygote advantage at the *PRNP* codon 129 that both increases the incubation period and reduces the susceptibility of infection to kuru.

## 2. Material and methods

### 2.1. Kuru model

A dynamic differential equation compartment model was used to analyse prion transmission. The population was stratified into three compartments: the number of susceptible individuals (*S*), infected (both symptomatic and asymptomatic) individuals (*E*), and dead individuals, infectious through cannibalism (*I*).
2.1
2.2
2.3where *N* = *S* + *E*.

The dynamics of the susceptible compartment are determined by parameters *α* *μ* and *β*. Parameter *α* is the *per capita* birth rate of susceptible individuals taking into account an annual population growth rate of 2 per cent [26]. The background *per capita* death rate, *μ*, of uninfected and alive infected individuals was derived from the approximately 50 year life expectancy of Papua New Guineans without kuru [27]. The transmission rate of kuru per susceptible individual per year is *β*. We assumed that no transmission occurred after 1959, as a result of the cannibalism ban, such that *β* = 0 after this date. Although it may be possible that cannibalism was practised after 1959, no deaths of people born after this date have occurred [21].

The dynamics of the exposed compartment are determined by parameters *β*, *μ* and *ν*. Parameter *ν* determines the virulence rate, and thus the rate at which individuals pass from the alive infected class (*E*) to the dead infectious class (*I*). This incubation period has been previously fitted to a Weibull distribution with mean 10 years [28]; however, a Weibull distribution does not fit the data significantly better than a more parsimonious exponential distribution (*p* = 0.05). Therefore, we assume an exponentially distributed incubation period with mean 1/*ν*.

The dynamics of the infectious compartment are determined by parameters *ν* and *δ*. Symptomatic and asymptomatic individuals are not distinguished in the *E* compartment because neurodegenerative symptoms do not have an impact on transmission. Individuals are assumed be infectious for a week, 1/*δ*. Unlike most transmission models, an infectious individual is already dead.

To estimate the transmission rate, *β*, and the virulence rate, *ν*, the model was fitted by the least-squares method to data on the annual mortality of kuru [29,30]. The distribution of the error between the best-fit model output and the incidence data yielded a standard deviation of 13 individuals (with a Lillie test accepting normality, *p* = 0.4) from a total of 1793 individuals. A credible interval for the model value was determined by estimating a joint posterior distribution for *β* and *ν* using Markov chain Monte Carlo; errors for the annual incidence of kuru were assumed to be normally distributed. Varying the standard deviation from 10 to 20 did not change the results. The Markov chains were run for a total of 5000 iterations, with a burn-in period of 5000 to estimate the posterior distribution of each parameter.

The exact date of the introduction of kuru is unknown [28]. However, cannibalism was introduced as a funeral rite of the southern Fore population around 1890 [31], providing a lower bound for the first transmission event. It was assumed that the epidemic occurred from a single index case before the first transmission event [19]. According to census data, the size of the Fore population was 11 991 in 1958 [24,25].

### 2.2. Epidemiological mechanism of host resistance

To integrate host heterogeneity into the model, the population was stratified by host genotype at codon 129: methionine/methionine (MM), methionine/valine (MV) or valine/valine (VV). Thus, in addition to equations (2.1)–(2.3) that model transmission dynamics, we incorporate equations (2.4)–(2.13) to include the population genetics of codon 129 such that
2.4
2.5
2.6
2.7
2.8
2.9
2.10where
2.11
2.12and
2.13It was assumed that both MM and VV homozygotes are epidemiologically identical (i.e. *β*_{VV} = *β*_{MM}), consistent with population genetic data [12,17]. We quantified the reduction in susceptibility to infection and duration of incubation period for a heterozygote relative to a homozygote. In the case that MV is at most as susceptible as VV or MM,
2.14

Conversely, in the case that MV is at least as susceptible as VV or MM,
2.15where 0 ≤ *f*, *g*, ≤ 1. Positive values of *f* and *g* in equation (2.14) or (2.15) correspond to relative resistance of a heterozygote or a homozygote, respectively. Under this suite of models (table 1), the host genotype affects either the transmission potential (setting *g* = 0, model T ‘transmission’) or the incubation period (setting *f* = 0, model I ‘incubation’) or both (model B ‘both’).

## 3. Results

### 3.1. Kuru model

To estimate the mortality rate and transmission rate of kuru, we fitted the dynamic model to kuru incidence data (figure 1). The mortality rate of kuru, *ν*, was estimated to be 0.081 (95% confidence interval (CI) 0.067–0.089) per year, corresponding to an exponentially distributed incubation period with a mean of 12.4 (95% CI 11.1–14.8) years. For a susceptible individual, the annual transmission rate, *β*, was estimated to be 0.015 (95% CI 0.014–0.018) per susceptible individual upon death of an infectious individual.

Because there is uncertainty regarding the size of the population at the start of the epidemic and the year of the index case, a sensitivity analysis was conducted using start dates between 1890 and 1930 and initial population sizes between 11 900 and 13 000. The best least-squares fit of the model to 1958 census data occurred when the death of the index case occurred in 1922 with a population size of 12 350. Therefore, we predict a population decline of 2.9 per cent over the period 1922–1958. To a precision of two significant figures, estimates for *β* and *ν* did not change during this sensitivity analysis of start date and population size.

In our model, the reproductive number for kuru is given by

Substituting the estimated values for *ν*, *β* and *N*(0) yields an *R*_{0} of 2.86 (95% CI 2.75–3.22) for the kuru epidemic.

If we considered the possibility that alive infected hosts who died of other causes could transmit, the model becomes
3.1
3.2
3.3In this case, *β* was estimated to be 0.012 per susceptible individual per year. The mortality rate *ν* remained at 0.081 per year, yielding *R*_{0} = *βN*(0)/*δ* = 2.85.

### 3.2. Epidemiological mechanism of host resistance

To infer the incubation periods and transmissibilities of the respective human genotypes, we fitted the epidemiological population genetic model for human resistance to the kuru incidence data (table 2). The *R*_{0} for kuru was calculated as the spectral radius of the next generation matrix of the model [32], assuming the initial population in 1922 was at Hardy–Weinberg equilibrium and the frequency of MV genotypes was 50 per cent [16]. All genotype models described the kuru mortality incidence data equally well, yielding similar corrected Akaike information criteria (AICc) values. Therefore, the genotype models could not be distinguished by the incidence data alone. However, these models could be distinguished in terms of their consistency with the population genetics data [16,33,34] (figure 2). The suite of eight models assumed resistance to kuru for heterozygotes or homozygotes manifested via either transmission rate (models Ti or Tii, respectively), or incubation period (Ii or Iii, respectively) or both mechanisms (Bi, Bii, Biii, Biv) with combinations of resistance/susceptibility for the heterozygotes and homozygotes.

Both the model structure (cf. figure 2*b*,*e*,*h*) and the susceptibility of heterozygotes compared with homozygotes (cf. figure 2*b*,*c*) affect the genotype frequencies qualitatively and quantitatively by resulting in different temporal trajectories and heterozygote percentages throughout the kuru epidemic. Furthermore, although the combinations of the free parameters in the model-fitting procedure significantly alter the expected gene frequencies, the free parameters can be selected such that the overall epidemic curve is maintained and the expected incidence of kuru is within the credible intervals we calculated (figure 2*a*,*d*,*g*). The relative frequencies of the genotypes affect the trajectory of the epidemic. For example, if heterozygotes are more resistant to infection than either homozygote (models Ti, Bi and Biii), less than 50 per cent of the kuru cases would be heterozygotes until transmission ceased in 1959. If genotype only affected susceptibility to infection, heterozygotes would constitute fewer than 50 per cent of kuru deaths before 1959. By contrast, if genotype also affected the incubation period, kuru incidence would increase among the heterozygotes over the latter half of the epidemic (1959–2010).

Accordingly, we compared the different model structures by evaluating whether the model results were consistent with the genotype frequency data. Specifically, we divided the results into three qualitative outcomes: whether (i) the fraction of MV genotypes is equal to about 0.5 (the genotype frequency in young modern Fore), or (ii) the fraction of MV genotypes is greater than 0.5, or (iii) the fraction of MV genotypes is less than 0.5, as corresponds to the snapshots of genotype frequencies in the early, mid and late epidemic (table 3). The model structures were then compared with these criteria at the times when they were sampled (table 4).

The genotype frequency snapshot in 1957–1958 from 142 individuals who died of kuru would have to show below 60 (42%) individuals or above 82 (58%) individuals who are heterozygous in order for the population to deviate significantly from 50 per cent heterozygosity (using a binomial model two-tailed test at 5%). Accordingly, we took the MV genotype frequency during these years to be approximately 50 per cent. Therefore, the genotype frequency snapshot in 1957–1958 was consistent with all of the models.

In contrast, the data from 1964 to 1988 included 13 individuals who died of kuru. There would have to be either below three (23%) heterozygotes or above 10 (77%) heterozygotes to conclude that the population deviates significantly from 50 per cent heterozygosity (using a binomial model two-tailed test at 5%). Accordingly, we conclude that the MV genotype frequency was below 50 per cent during this middle epidemic period.

Our results show that the genotype frequency snapshot in 1964–1988 was consistent with five models: (i) heterozygotes being resistant to infection (figure 2*b*, Ti), (ii) heterozygotes having a shorter incubation period than homozygotes (figure 2*f*, Iii), (iii) heterozygotes being resistant to infection and having a longer incubation period (figure 2*h*, Bi), (iv) heterozygotes being susceptible to infection but having a shorter incubation period (figure 2*k*, Biii), and (v) heterozygotes being susceptible to infection and having a shorter incubation period (figure 2*l*, Biv). Three of the models can be eliminated on the basis of this comparison (models Tii, Ii, Bii, respectively, figure 2*c*,*e*,*i*).

Given that in 1996–2001 10 individuals died from kuru, there would have to be either below two (20%) heterozygotes or above eight (80%) heterozygotes in order for the population to deviate significantly from 50 per cent heterozygosity. Therefore, these data suggest that more than 50 per cent of the individuals who died of kuru in those years were heterozygotes.

All together, these empirical genotype frequencies are narrowed down to four models: (i) heterozygotes being more susceptible to infection than homozygotes (figure 2*c*, Tii), (ii) heterozygotes having a longer incubation period than homozygotes (figure 2*e*, Ii), (iii) heterozygotes being resistant to infection and having a longer incubation period than homozygotes (figure 2*h*, Bi), and (iv) heterozygotes being resistant to infection but having a shorter incubation period than homozygotes (figure 2*i*, Bii). Three of these models (Tii, Ii, Bii) have already been eliminated as being inconsistent with the 1964–1988 data.

Overall, our findings collectively demonstrate that only one hypothesis (Bi) could account for all the data. This hypothesis predicts that the fraction of heterozygote kuru patients was slightly less than 50 per cent before 1988, rising to more than 50 per cent after 1996. Our calculations indicate that heterozygotes have a longer incubation period and lower infection rate than homozygotes, such that decreasing the resistance advantage of heterozygotes increases the percentage of heterozygotes throughout the epidemic period.

## 4. Discussion

To evaluate both the epidemiological parameters for kuru and the epidemiological role of human resistance to kuru, we combined a dynamic model of kuru transmission with a population genetics model of human resistance to kuru. We estimated the reproductive number for kuru within the Fore population of Papua New Guinea together with the most likely start date of the epidemic. Our findings show that heterozygotes were resistant to kuru with both a reduced susceptibility to infection and an extended incubation period compared with homozygotes.

We fitted our model to data on the trajectory of the kuru epidemic to estimate the incubation period, transmission rate and *R*_{0} of kuru in the Fore population over the last hundred years. We estimated that the first death due to kuru was in 1922, when the population size of the Fore people was around 12 350. Our estimate of the time of the index case is consistent with anecdotal reports of the first cases within the region [31]. Our calculations show that the average incubation period for kuru was 12 years and that the *R*_{0} was 2.86.

The transmissibility, *R*_{0} and initiation date of kuru had not been estimated previously, whereas the incubation period was previously estimated to have a mean of 10.3–13.2 years using back-calculation methods [28] and approximately 12 years based on the termination of transmission as well as the continued occurrence of infection [22]. Our estimate of 12 years with a credible interval of 11–15 years is consistent with both of these studies.

We extended our epidemic model to a suite of alternative models incorporating the human population genetics of kuru resistance to evaluate the proposed hypotheses of resistance mechanisms for the *PRNP* gene. Specifically, we compared predictions of genotype frequencies from these different models with human population genetic data sampled during the kuru epidemic. The only model consistent with the genotype data assumes that heterozygotes are resistant to infection and have a longer incubation period than homozygotes. This is also consistent with other studies that have hypothesized that genotype MM confers an increased infection rate and shorter incubation period than the other genotypes [34], as has been found in other prion diseases [12].

Our mathematical model is a simplification of the kuru dynamics within the East Highland region of Papua New Guinea. The incidence is known to have varied between different linguistic groups in the Eastern Highlands area, such as between the South Fore and the Northern Fore [35], within which kuru spread from village to village, some more affected than others [31]. Although there have been reports of kuru entering the Fore population around 1920, some Fore villages did not describe any kuru mortality until the 1940s, suggesting a possible travelling wave of disease. Owing to the limited spatial stratification of the kuru data, our mathematical model does not capture this temporal spread and is therefore an approximation of the dynamics throughout the entire region. The likely effect of assuming spatial homogeneity when transmission is actually spatially heterogeneous may be an elevation in the estimation of contact rate. Nonetheless, the cumulative amount of transmission would not be affected by such spatial heterogeneity. Consequently, if we assume that the true dynamics of kuru are more consistent with a travelling wave epidemic, the genotype frequencies observed are the average frequency of the total population, with each subpopulation contributing different frequencies at any one time. Furthermore, assuming that all individual genotypes behave according to the simple genotype model, the resistance hypotheses analysed here are valid. If the different villages that experienced the travelling wave at different times comprised individuals of different genotype frequencies prior to the epidemic, our conclusions would be affected. However, recent genotyping of young modern Fore and young individuals from the East Highland Papua New Guinea region that was not affected by kuru did not show a significant deviation from the 25/50/25% MM/MV/VV genotype frequency [16], suggesting that it was unlikely that there were differences in the genotype frequencies between villages initially. Additionally, there is evidence that the age of onset of kuru is strongly correlated with the probability that the patient is heterozygote [12,33]. These data suggest that heterozygote kuru patients, irrespective of the time or place of kuru infection, had longer incubation periods.

Our model does not account for differences in either contact or susceptibility between children and adults or women and adults. It has been established that women and children were affected predominantly at the start of the epidemic, and men towards the end of the epidemic, as men were only exposed to infection during childhood [35]. Sample sizes for genotype data in the later stages of the epidemic are sparse, leading to larger standard errors around the central estimates of the genotype frequencies. Thus, the conclusions drawn as a result of the frequency data from 1964 onwards are less reliable than those arising from frequency data prior to 1964.

Our results suggest that kuru may have imposed short-term balancing selection pressure on humans, where a fitness advantage for heterozygosity caused an increase in the frequency of heterozygotes in the population during the kuru epidemic. Although genome scans have indicated numerous potential targets of balancing selection [36,37], there are only a few examples of balancing selection imposed by a known pathogen in humans, including the major histocompatibility complex/human leucocyte antigen locus, the *β*-haemoglobin locus and the G6PD locus [15].

In summary, our study is the first to estimate the incubation period of kuru using a dynamic model and the first study to estimate the transmissibility of kuru and its reproductive number. Moreover, we also test the hypothesis of genetic susceptibility to kuru using genotype frequency data. Through this approach, we are able to evaluate the epidemiological mechanism of resistance and to demonstrate a signal of balancing selection, by which kuru-resistant heterozygotes were selected because of their increased fitness during the kuru epidemic. We have revealed in this study epidemiological mechanisms by which this balancing selection could have acted. These results have implications for understanding the resistance mechanisms of other prion or prion-like diseases.

## Acknowledgements

We thank Angelika Hofman for editing assistance and are grateful for the input of two anonymous reviewers. K.E.A. was funded by the Miriam Weston Trust.

- Received April 11, 2013.
- Accepted May 16, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.