Estimates of the changing age-burden of Plasmodium falciparum malaria disease in sub-Saharan Africa

Estimating the changing burden of malaria disease remains difficult owing to limitations in health reporting systems. Here, we use a transmission model incorporating acquisition and loss of immunity to capture age-specific patterns of disease at different transmission intensities. The model is fitted to age-stratified data from 23 sites in Africa, and we then produce maps and estimates of disease burden. We estimate that in 2010 there were 252 (95% credible interval: 171–353) million cases of malaria in sub-Saharan Africa that active case finding would detect. However, only 34% (12–86%) of these cases would be observed through passive case detection. We estimate that the proportion of all cases of clinical malaria that are in under-fives varies from above 60% at high transmission to below 20% at low transmission. The focus of some interventions towards young children may need to be reconsidered, and should be informed by the current local transmission intensity.


Transmission model
The transmission model which we fitted to the prevalence and clinical incidence data was in most respects the same as a previously published model 11 , but is described here for completeness. At each point in time people can be in one of six states: susceptible ( S ), treated clinical disease (T ), untreated clinical disease ( D ), asymptomatic infection which may be detected by microscopy ( A ), sub-patent infection (U ) and protected by a period of prophylaxis from prior treatment ( P ). Latent liver stage infection is modelled as a delay in the force of infection of duration E d . People move between these states as shown in Figure 6. The partial differential equations for the human dynamics are as follows, with t representing time and a age, is the probability that clinical malaria is effectively treated, and each d is the mean duration in a state . Rather than numerically solve the partial differential equations, we stratified the infection states by age, and also by degree of exposure to mosquitoes. The equilibrium states for a given EIR can then be found by iteratively going through the age groups from youngest to oldest.

Heterogeneity in biting rates
We assume that each person has a relative biting rate  following a log-normal distribution between people with parameters 2 /2   and  , so that has a mean of 1. With a mean EIR for adults of 0  , the EIR  and force of infection  at age a are given by where b is the probability of infection if bitten by an infectious mosquito, and   and 0 a determine how the biting rate changes with age. To model this heterogeneity using a compartmental model, the infection states are stratified by exposure as well as age. Suppose that there are n exposure categories. Let 1 ,..., n xx and 1 ,..., n ww be the Gauss-Hermite abscissas and weights for integrating a function multiplied by a standard normal probability density 71 . We take a proportion i w of the population to be in exposure category i and set the relative biting rate in this category to be Then the EIR in the category i is

Immunity functions
We consider three points at which immunity may act:  a reduction in the probability of infection following an infectious challenge In our previously published model, the effect of blood-stage immunity was instead to reduce the duration of infection.
The three types of acquired immunity increase with exposure or wane as follows, starting from zero at birth: 1 Here, each u represents the time during which immunity cannot be boosted after a previous boost and each d governs the duration of immunity. The immunity functions are stratified by age and exposure to mosquitoes.

13
The probabilities of infection, detection and clinical disease are given by Hill functions, as these have a sufficiently flexible shape for a quantity which is assumed to decrease from its value with no immunity to a minimum value with high immunity.
The probability of infection is:  determine the probabilities that states A and U are detected by PCR, which are given by A q  and U q  respectively.

Infectiousness to mosquitoes
The lower probability of detection is assumed to be due to a lower parasite density, which also reduces onward infectivity. For a given set of parameters, including the EIR for each site, let the model-predicted incidence of clinical malaria in site j and age group i be ji  , and the person-time at risk and number of events be ji T and ji y respectively. We also include a parameter r for the method of case detection, with 1 r  for daily ACD and less than 1 for weekly ACD or PCD. We then assume that The rationale for this formulation of the age-specific random effects is that if the age groups 1,..., 4 k  coincided with the age groups in the study, then the Poisson-Gamma model would be equivalent to a negative binomial model (conditional upon the random effects u ). However, there was finer age stratification in some studies than others, and simply having a negative binomial likelihood for each age group i would give different statistical weight to a study depending on how finely the data were grouped.
The overall likelihood for clinical incidence data in a study is thus   where f is a log-normal probability density, g is a Gamma probability density and P is a Poisson probability. The random effects jk v can be analytically integrated out of this likelihood 72 , and we numerically integrated the random effects u out by adaptive Gauss-Hermite quadrature. Using numerical integration rather than sampling u in the MCMC fitting has a number of advantages: it reduces the number of quantities sampled; the EIRs for each site were sampled, and these are correlated with the random effects; and the model equilibrium only needs to be found once, and then in the numerical integration the same equilibrium value is used repeatedly. The accuracy of the integration was checked by increasing the number of quadrature points and ensuring that the change in the posterior parameter estimates was negligible. and invlogit its inverse function, and  is an over-dispersion parameter. We integrated the random effects w out of the likelihood by adaptive Gauss-Hermite quadrature. Data were stratified into age groups with boundaries 0, 2, 5, 10, 15, 20, 30, 40 years if they were presented in finer age groups than this: unlike clinical incidence, parasite prevalence does not have a sharp peak at young ages, hence less information is lost by grouping the data.
The infectivity parameters were fitted alongside the rest of the model. The data for infectivity consisted of three published studies in which mosquitoes were fed on a sample from the population of endemic areas, regardless of the volunteers' infection status 54-56 . In each study, numbers of mosquitoes feeding and infected were recorded by age group, and there was information on parasite prevalence. A Betabinomial likelihood was used with over-dispersion parameter I  .

Prior distributions and posterior estimates
The prior distribution for the log of the EIR in each of the datasets was taken as normally distributed with the standard deviation chosen according to whether or not the estimate was in the same year, and whether it was from the same village or town or from elsewhere in the same region (Supplementary Table S2). Data on treatment rates by country were obtained from the relevant Demographic Health Surveys (DHS) and Malaria Indicator Surveys (MIS) 41 . These surveys have data on the proportion of fever cases in under-fives which were treated with an antimalarial drug, which we took to be the treatment rate. As estimates were variable between surveys in the same country, we chose compromise values based on surveys near in time (Table 1). The prior estimate for the treatment rate was taken as 90% in the closely monitored cohorts in Dielmo and Ndiop, Senegal, and 5% for data from the Garki study. The prior distribution for the proportion of clinical cases effectively treated was a Beta distribution with mean at 75% of the DHS-derived value, and sum of Beta shape parameters = 25, where the factor 75% accounts for imperfect efficacy of pre-ACT treatments 73 .
Prior distributions for all other parameters in the model are similar to those previously reported 11 . Supplementary Table S1 shows the prior and posterior estimates for all the model parameters if the parameter was estimated. Some parameters were fixed based on values from the literature, and some shape parameters were given moderately informative prior distributions in order to ensure identifiability. In particular, the parameter governing heterogeneity in biting ( ) was fixed. This parameter has a substantial impact on the basic reproduction number and hence on the effectiveness of interventions. When we tried fitting the parameter, its posterior value differed greatly from the prior value, even though it is clear that the data does not contain genuine information about heterogeneity in biting. The observed and model-predicted probability of mosquito infection by age of the human population are shown for each study in Supplementary Figure S5.

Spatial Data Sources and Modelled Relationships
The methods for estimating the Africa-wide burden of clinical malaria in 2010 based on the parasite prevalence in 2 to 10 year-olds are summarised in the main text. Here we give further details.
To estimate the population by age across sub-Saharan Africa, we started with Landscan estimates of the population at 5km resolution in 2007 43 ; we then multiplied these by country-specific growth rates from the World Bank to obtain estimates of the population in 2010; and finally, used the estimates of the age distribution in each country from the UN to estimate the proportion of people in each 5km square and each five-year age group.
In each 5km square in sub-Saharan Africa we drew a sample of 10,000 prevalence values from the posterior distribution in that square: this distribution is from the model fitted by the Malaria Atlas Project (MAP) 42 . We then took a sample of 10,000 parameter sets from the posterior distribution of our fitted model, and for each pair of prevalence value and model parameter set we calculated the incidence that would be detected by different case detection methods, by five-year age groups, using the equilibrium model solutions. We then multiplied the incidence by the estimated number of people in each age group in that 5km square, to give a posterior distribution of the number of cases in each age group. This gives a posterior distribution of the proportion of cases by age at 5km resolution, or the numbers of cases by age group can be aggregated to give a posterior distribution of the Africa-wide burden and age-distribution of cases. Our model is not well-validated for prevalences above 80%, and needs implausibly high EIRs to reach prevalences above 85%. Hence when sampling from the posterior distribution of prevalence in each square, when the sampled value was above 80%, we set the clinical incidence to be the same as it would be at a prevalence of 80%.
There is wide uncertainty in the proportion of clinical cases in each age group across Africa (Supplementary Figure S3), due primarily to the uncertainty in the underlying prevalence estimates. If instead of the full posterior distribution, we use the posterior mean prevalence from the MAP model, so that the only uncertainty is in our model estimates, the proportion of cases in each age group is much more precisely estimated (Supplementary Figure S4).

Sensitivity analysis to immunity model
The relationship between clinical incidence and prevalence is somewhat sensitive to the choice of immunity model (Supplementary Figure S1B). The model presented in the main text (blue line) assumes that immunity against clinical disease increases with exposure, but with a delay after each boost to immunity during which immunity cannot be further boosted. This delay was included for biological plausibility. In this model clinical incidence always increases with prevalence. We fitted an alternative model in which immunity to disease is boosted at every exposure, so that in equation Error! Reference source not found. CA I increases as follows: Now the fitted incidence decreases at higher prevalences (red line). The model the main text incorporates a pre-erythrocytic/liver-stage immunity against infection and also models the effect of blood-stage immunity in reducing the probability of detection, with the choice of model based on the demonstrated partial efficacy of the RTS,S vaccine and on estimates of the detectability of parasite clones from genotyping studies 22,40 . Refitting with either of these removed also changes the incidence-prevalence curves (green and yellow lines). Incorporating immunity that reduces the duration of infection to a minimum duration of patent infection of 50 days made little difference to these relationships (not shown).

Sensitivity analysis to treatment rate
For model predictions and burden estimates, we assumed a treatment rate of 40% (the proportion of clinical cases treated with an effective drug), as the DHS and MIS surveys do not have data for every country, whereas we used country-specific data when fitting the model. A higher EIR is required to explain a given prevalence if there is a high treatment rate (Supplementary Figure S6a). Hence if the model is calibrated to match a given prevalence, there is a higher incidence if the treatment rate is higher (Supplementary Figure S6b). This results in the estimated Africa-wide burden increasing with the assumed treatment rate (Supplementary Table S3). However the estimated proportion of cases that occurs in under-fives is less sensitive to this assumption.

Sensitivity analysis to seasonal variation
Supplementary Figure S2 shows how the results are affected if there is seasonal variation in transmission, with the mosquito density varying over the year in proportion to the curves in Supplementary Figure S2d, with the same pattern repeating each year. The incidence and prevalence shown are the mean over the year. With the model presented in the main text, the incidence in under-fives is higher when there is marked seasonality than in a perennial transmission setting (Supplementary Figure S2a), whereas the proportion of cases in under-fives is less affected (Supplementary Figure S2b). In a model in which immunity to clinical disease increases in proportion to exposure (the model plotted with red dashed lines in Supplementary Figure S1), seasonality makes less of a difference to the incidence in under-fives (Supplementary Figure S2c), except at high prevalences when seasonal variation results in lower incidence for a given prevalence.