Modelling pathogen load dynamics to elucidate mechanistic determinants of host-Plasmodium falciparum interactions.

During infection, increasing pathogen load stimulates both protective and harmful aspects of the host response. The dynamics of this interaction are hard to quantify in humans, but doing so could improve understanding of the mechanisms of disease and protection. We sought to model the contributions of the parasite multiplication rate and host response to observed parasite load in individual subjects infected with Plasmodium falciparum malaria, using only data obtained at the time of clinical presentation, and then to identify their mechanistic correlates. We predicted higher parasite multiplication rates and lower host responsiveness in cases of severe malaria, with severe anaemia being more insidious than cerebral malaria. We predicted that parasite-growth inhibition was associated with platelet consumption, lower expression of CXCL10 and type 1 interferon-associated genes, but increased cathepsin G and matrix metallopeptidase 9 expression. We found that cathepsin G and matrix metallopeptidase 9 directly inhibit parasite invasion into erythrocytes. The parasite multiplication rate was associated with host iron availability and higher complement factor H levels, lower expression of gametocyte-associated genes but higher expression of translation-associated genes in the parasite. Our findings demonstrate the potential of using explicit modelling of pathogen load dynamics to deepen understanding of host-pathogen interactions and identify mechanistic correlates of protection.

Introduction

Improved methods are needed to identify mechanisms which protect against human infectious diseases in order to develop better vaccines and therapeutics1,2. Pathogen load is associated with the severity of many infections3, and is a consequence of how fast the pathogen can replicate, how long the infection has been ongoing, and the inhibition or killing of pathogen by the host response (Fig. 1a). The contribution of these factors varies within an individual over the course of infection, as well as between individuals. Identifying mechanistic correlates of the processes which determine pathogen load is likely to be more informative than identifying correlates of pathogen load per se. However, in humans the timing of infection is rarely known and treatment cannot usually be withheld to observe the natural dynamics of pathogen load and host response. Here we present an approach to estimate the latent determinants of parasite load dynamics. We use these estimates to better understand severe malaria phenotypes and to identify mechanisms inhibiting parasite growth and controlling parasite multiplication during Plasmodium falciparum malaria in Gambian children.

Fig 1

Estimating the dynamics of parasite load and host response in malaria.

(a) Conceptual model of determinants of parasite load. (b) Schematic of relationships between parasite load, multiplication rate (m), P, and parasite growth inhibition (PGI) derived from the longitudinal malariatherapy dataset. (c) Correlation matrix for P, parasite biomass, duration of illness and TNF concentrations in 2000 simulated Gambian children (Spearman correlation, LOWESS fit lines). (d) Performance in simulated subjects of the best models to predict ln P and m, compared with predictions made using individual variables only. Boxes show median and interquartile range, whiskers extend 1.5-times the interquartile range or to limit of range, n=100 simulated datasets (each of 139 subjects). (e-i) Comparisons of parasite biomass (e), TNF (f), duration of illness (g), predicted m (h), predicted P (i), in 139 Gambian children with uncomplicated (UM, n=64) or severe malaria (SM1, prostration, n=36; SM2, any combination of cerebral malaria, hyperlactatemia or severe anemia, n=39). Box and whiskers as in d; P for Kruskal-Wallis (above plots) and Mann-Whitney tests (UM vs SM2, within plot). (j, k) Correlation of predicted m (j) or P (k) with age, P for Spearman correlation, n=139.

Results

Estimating determinants of parasite load and host response dynamics in humans

To estimate the determinants of parasite load dynamics in naturally-infected malaria patients we calibrated a statistical prediction model using outputs from a mechanistic simulation which combined information from two datasets. A historical dataset of the longitudinal course of untreated infection in 97 patients who were deliberately inoculated with P. falciparum as a treatment for neurosyphilis (malariatherapy dataset) (Supplementary Fig.1) was used as a reference for changes in parasite load over time4. A dataset from 139 naturally infected Gambian children with malaria (Gambian dataset, Supplementary Table 1, Supplementary Dataset 1) was used for subsequent discovery of the determinants of parasite load dynamics. We used an existing mathematical model for the malariatherapy data (the Dietz model4), which estimated latent variables thought to determine changes in parasite load over time in each individual, and modified the model to better represent the features of the Gambian dataset. We used the modified model to simulate a large number of in-silico Gambian patients, with all latent variables and course of infection fully known, and then trained a statistical model to learn from these simulations the relationships between variables available in the real Gambian patient data and the unobservable, latent variables.

In the models4, the increase in parasite load up to the first peak is determined by two individual-specific latent variables (Fig. 1b, see Methods): the within-host multiplication rate, m, which is the initial rate of increase in parasite load before any constraint by the host response; and the parasite load required to stimulate a host response that reduces parasite growth by 50%, P,4. When m, P, and parasite load are known, parasite growth inhibition (PGI) by the host response can be calculated (see Methods). We allowed rescaling of P values between the malariatherapy and Gambian datasets, and incorporated plasma Tumour Necrosis Factors (TNF) concentrations as an indicator of the protective host response5,6, using a maximum-likelihood approach (see Methods and Supplementary Fig. 2). These modifications resulted in higher P values in the Gambian population than malariatherapy subjects, consistent with epidemiological data showing higher fever thresholds in P. falciparum infected children than in adults7. Other model assumptions and definitions are shown in Supplementary Table 2.

To test whether combination of a mechanistic simulation model with statistical learning of the relationships between latent and directly observable variables was better at predicting the determinants of parasite load than using observable variables alone, we simulated 2000 Gambian children with malaria with known values of m, P parasite biomass, duration of illness and plasma TNF (Fig. 1c and Supplementary Fig. 3) and then fit general additive models (GAMs) to predict values of m and P for individual children (Supplementary Table 3). The resulting models produced more accurate predictions of of m and P than using individual variables alone (Fig. 1d).

Next we used the GAMs to predict values of P and m for each of the 139 individuals in the Gambian dataset (Fig. 1e-k, Supplementary Fig. 4). Children with the most severe manifestations of malaria (SM2) had the highest parasite load, TNF, predicted m, and predicted P values, intermediate values were seen in those with prostration as the only manifestation of severe disease (SM1), and values were lowest in uncomplicated malaria (UM), whilst duration of illness did not differ significantly by clinical phenotype (Fig. 1e-i). These observations suggest that high parasite load and severe disease are most likely in individuals with either fast replicating parasites (high m) or less immune responsiveness (high P).

Since age can be a major determinant of malaria severity and naturally acquired immunity8, we examined whether age was associated with m or P. Age was not significantly correlated with m but was significantly negatively correlated with P (Fig. 1j,k). This implies little age-related acquisition of constitutive resistance (for example, naturally-acquired antibody-mediated immunity) in these children, as might be expected from the relatively low malaria transmission in this region of The Gambia9. However, these data also indicate that a lower parasite load would be needed to provoke an equivalent host response in older individuals without significant naturally acquired immunity.

Predicting severe malaria phenotype from within-host dynamics

We next asked whether individual estimates of m and P could be used to predict pathophysiological features malaria which had not been used in our model derivation. Severe malarial anemia (hemoglobin concentration <5g/dL), is most common in the youngest children in high transmission settings, but rare in lower transmission settings such as Greater Banjul region of The Gambia, where cerebral malaria was relatively more common10. Severe malarial anemia is characterised by a higher parasite biomass10–12, lower levels of both TNF and interleukin-10 (IL-10), but an elevated ratio of TNF:IL-1013,14 when compared to cerebral malaria. In our Gambian subjects, hemoglobin concentration could be predicted from estimated P, m and age; IL-10 concentration could be predicted from m and Pc (Supplementary Table 4, Fig. 2a-b). We simulated a population of Gambian infants, selected those predicted to have hemoglobin <5 g/dL, and compared their characteristics to real Gambian subjects with cerebral malaria. The simulated severe anemia cases had lower m but similar P, higher parasite biomass and longer duration of illness than the cerebral malaria patients (Fig 2c-f). Both TNF and IL-10 concentrations were predicted to be lower in severe anemia than in cerebral malaria (Fig 2g-h), whereas the TNF:IL-10 ratio was predicted to be higher in severe anemia (Fig 2i), supporting the biological plausibility of relationships defined in our model and illuminating a potential explanation for these distinct severe malaria phenotypes.

Fig 2

Contribution of parasite load dynamics to severe malaria phenotype.

(a, b) Comparison of predicted and actual hemoglobin (a, n=136) and IL-10 (b, n=139) concentrations in the Gambian children. Pearson correlation, shaded region, 95% CI of regression line. (c-i) Comparisons of m, P, parasite biomass, days of illness, plasma TNF, plasma IL-10, and plasma TNF:IL-10 ratio, in Gambian children with cerebral malaria (CM, n=12) and simulated Gambian infants with severe anemia (SA, n=24). Boxes show median and interquartile range, whiskers extend 1.5-times the interquartile range or to limit of range.

Estimating parasite growth inhibition reveals the protective effect of platelets

The role of the host response in restricting parasite load is often unclear in human malaria because the strongest host responses are often seen in patients with the highest parasite loads and most severe disease15,16. For example platelets directly inhibit parasite growth16,17, and the reduction in platelet count typically seen in malaria is partly a consequence of the protective mechanism of platelet adhesion to infected red cells16. However the reduction in platelet count is greatest in individuals with the highest parasite load and most severe disease18, which seems counterintuitive if the low platelet counts indicate parasite killing. In Gambian children, estimated PGI did not differ significantly by clinical phenotype (Fig. 3a) indicating that the components of the host response which restrain parasite growth are similarly activated in severe and uncomplicated disease groups at the time of hospital presentation, but implying that this response developed too late to prevent high parasite load in the severe cases. Subjects with severe disease had the lowest platelet counts (Fig. 3b and Supplementary Table 1) and the highest parasite loads (Fig. 1d), but the protective role of platelets was evident through the significant (P=0.0001) correlation with PGI (Fig. 3c). Thus considering differences between individuals in observed parasite load and host response as part of a dynamic rather than static process can resolve counterintuitive associations.

Fig 3

The protective effect of platelets is revealed by estimating parasite growth inhibition.

(a,b) Comparisons of PGI (a) and platelet count (b) in 139 Gambian children with uncomplicated (UM, n=64) or severe malaria (SM1, prostration, n=36; SM2, any combination of cerebral malaria, hyperlactatemia or severe anemia, n=39 (platelet data missing from 4 subjects)). (c) Correlation between platelet count and PGI (n=135) shows that low platelet count is associated with greater parasite growth inhibition. Boxes show median and interquartile range, whiskers extend 1.5-times the interquartile range or to limit of range; P for Kruskal-Wallis (above plots) test (a, b) and for Spearman correlation (c).

Predicting mechanistic correlates of parasite growth inhibition

To determine whether our model-derived estimates could be used to discover aspects of host-parasite interaction we sought to identify mechanistic correlates of protection and susceptibility. We analysed human whole blood gene expression, with gene signature-based deconvolution to adjust for leukocyte-mixture19, from samples of 24 children at the time of presentation (13 with UM, 11 with SM, Supplementary Table 5). Of 11702 detected human genes, 51 were significantly correlated (26 positively, 25 negatively) with estimated PGI after adjustment for false discovery rate (Benjamini-Hochberg adjusted P<0.05, Fig. 4a, Supplementary Table 6). We reasoned that genes positively correlated with PGI should be enriched for effector mechanisms which act to reduce parasite load, whilst genes negatively correlated with PGI should be enriched for mechanisms which favour increase in parasite load. Eight of these genes were also correlated with parasite biomass and three with TNF (Supplementary Table 6).

Fig. 4

Transcriptional correlates of parasite growth inhibition

(a) Volcano plot showing association between gene expression and parasite growth inhibition after adjustment for leukocyte mixture in a linear model. Log fold change (log FC) is the coefficient of log adjusted gene expression vs. parasite growth inhibition. Positive log FC indicates that increasing gene expression is associated with increasing parasite growth inhibition. Negative log FC indicates that increasing gene expression is associated with decreasing parasite growth inhibition. P calculated using two-sided likelihood ratio test, adjusted for multiple testing using the Benjamini-Hochberg method: false discovery rate adjusted P <0.05 (FDR) is considered significant (above dashed line, colored circles). The 10 significant genes with greatest positive and negative log FC are labelled. Analyses based on data from n=24 subjects. (b,c) Primary networks derived from the genes significantly associated with PGI, with positive (b, n=26) and negative (c, n=25) log FC.

Genes positively correlated with PGI (Fig 4a) showed limited canonical pathway enrichments (Supplementary Table 7) but 16 (62%) were linked together in a network around extracellular signal-regulated kinases ERK1/2 and AKT serine/threonine kinase (Fig. 4b). These kinases integrate cellular inflammatory and metabolic responses to control innate defence mechanisms such as cytokine secretion, phagocytosis and degranulation20,21. The 25 genes negatively correlated with PGI were strongly enriched in immune response pathways (Supplementary Table 7). Network analysis showed 15 (60%) of the negatively correlated genes were linked through a network focussed around type 1 interferon (Fig. 4c), consistent with observations that sustained type 1 interferon signalling is associated with higher parasitemia in mice22–25 and potentially in humans22,26. C-X-C motif chemokine ligand 10 (CXCL10, also known as IFN-γ-inducible protein of 10 kDa, IP-10) expression had the greatest log-fold change of the genes negatively correlated with PGI (Fig. 4c), consistent with findings that CXCL10 deletion and neutralisation decrease parasite load in mice27.

We investigated whether associations with PGI were dependent on assumptions we made about the true severity rate in Gambian children, which we assumed to be 5% based on published data in other settings28,29. Varying this to credible extremes of 1% and 10% and repeating the process of calibration between datasets, fitting of models to predict m and P, and estimating new values for PGI, resulted in little difference in the genes identified as significantly associated with PGI, or the significance of individual genes (Supplementary Table 8).

Cathepsin G and MMP9 directly inhibit parasite growth

The 26 genes positively correlated with PGI have not, to our knowledge, previously been described as having anti-parasitic effects so we sought direct biological evidence, focussing on two encoding secreted proteins as the best candidates: CTSG (cathepsin G) and MMP9 (matrix metallopeptidase 9, also known as matrix metalloproteinase 9 and gelatinase B), which both encode neutrophil granule proteins30. We tested whether these proteases could inhibit parasite growth in vitro. Cathepsin G and MMP9 both inhibited growth of P. falciparum 3D7 strain (Fig. 5a). Addition of cathepsin G to schizont cultures produced a dramatic reduction in invasion of new erythrocytes, and pretreatment of erythrocytes with cathepsin G before adding them to schizont cultures produced a similar reduction in their invasion (Fig. 5b), indicating that cathepsin G acts primarily on the erythrocyte. Addition of MMP9 to schizont cultures produced a more modest reduction, whilst pretreatment of erythrocytes did not reduce invasion, implying that MMP9 likely acts against schizonts or free merozoites rather than preventing invasion at the erythrocyte surface (Fig. 5b).

Fig. 5

Effects of cathepsin G and MMP9 on parasite growth and expression of erythrocyte invasion receptors

(a) Effect of cathepsin G (18μg/mL, n=5) and MMP9 (16μg/mL, n=3) or no treatment (n=8) on in vitro growth of P. falciparum 3D7 (n are biological replicates, results representative of two independent experiments). (b) Effect of cathepsin G (18μg/mL) and MMP9 (18μg/mL) on erythrocyte invasion of P. falciparum 3D7 when added directly to schizonts and donor red cells, or when pre-incubated (PT) with donor red cells before washing and adding to schizonts (n=3 biological replicates per condition, representative of two independent experiments). (a, b) Show mean (95% CI) and P for two-sided unpaired t-test. (c) Cathepsin G and MMP9 concentrations in plasma from healthy donor whole blood (n=8) unstimulated or stimulated with 1μM PMA, and from Gambian children with P. falciparum malaria (n=34). Bars show median, P for two-sided Wilcoxon matched pairs test. (d-e) Dose effects on growth inhibition by MMP9 against P. falciparum 3D7 (d), and invasion inhibition by cathepsin G pre-treatment against four parasite strains (e) (n=3 biological replicates per dose, mean (95% CI) and P for linear trend, each result representative of two independent experiments). (f) Additive effect of Cathepsin G 1μg/mL and MMP9 1μg/mL against P. falciparum 3D7 invasion (n=4 biological replicates per condition, mean (95% CI) and P for ANOVA, representative of three independent experiments). (g) Dose response for erythrocyte surface receptor cleavage by cathepsin G (n=3 biological replicates per dose, mean +/- standard error, asymmetric 5-parameter logistic regression fit lines, representative of two experiments). (h) Effect of PMA stimulation of healthy donor (n=8) whole blood on erythrocyte surface receptor expression assessed by fluorescence intensity (P for two-sided Wilcoxon matched pairs test). (i) Comparison of proportion of erythrocytes with detectable receptor expression in acute (day 0) and convalescent (day 28) samples from Gambian children with malaria (n=6, P for one-sided Wilcoxon matched pairs test).

In order to identify biologically relevant concentrations of cathepsin G and MMP9 we measured their concentrations in whole blood from healthy donors, before and after stimulating degranulation, and in plasma from children with malaria at the time of clinical presentation (Fig. 5c). Local concentrations which might occur in vivo, adjacent to degranulating neutrophils, could be at least an order of magnitude higher31. MMP9 is also known to be released from other cell types in response to P. falciparum, including vascular endothelial cells32. MMP9 dose-dependently inhibited parasite growth over a physiological range of concentrations (Fig. 5d). Similarly, parasite invasion was dose-dependently inhibited by cathepsin G pre-treatment of erythrocytes, with similar effects in each of four parasite strains with different invasion phenotypes33 (Fig. 5e). Combined treatment with low doses of MMP9 and cathepsin G – in the range detected in patient plasma – showed an additive effect (Fig. 5f).

Cathepsin G has previously been reported to cleave red cell surface glycophorins34, therefore we asked whether it might also cleave other RBC surface proteins which are used as invasion receptors by P. falciparum. Consistent with its broad inhibition of parasite invasion, cathepsin G dose-dependently cleaved the majority of P. falciparum invasion receptors including glycophorins A, B, and C, CD147 (basigin), CD108 (semaphorin 7A), and complement receptor 1 (CR1), but not CD55 (DAF) (Fig. 5g). MMP9 did not cleave any of these surface receptors (Supplementary Fig. 5). PMA stimulation of healthy donor whole blood recapitulated the loss of erythrocyte surface glycophorins A and B, CD108 and CD147 in all donors, decreased glycophorin C expression in 6 of 8 healthy donors, but did not consistently reduce CR1 (Fig. 5h) (as might be expected from the dose-response curves, Fig 5g). In samples from Gambian children on the day of presentation with P. falciparum malaria, the proportions of erythrocytes with detectable expression of glycophorins A and B and CD147 were significantly lower than in convalescent samples (28 days after treatment), and there was a trend to lower expression of CD108 and glycophorin C (Fig. 5i). These results would be consistent with cleavage of these surface molecules in vivo during acute infection. The variable expression seen at day 28 (Fig. 5i) may indicate the persistence of modified erythrocytes in the circulation. The importance of glycophorins and basigin in RBC invasion and genetic susceptibility to severe malaria is well established36–38, and so it is highly likely that the cleavage of these erythrocyte receptors by cathepsin G would have a protective effect in vivo.

Host and parasite factors associated with parasite multiplication rate

In our model, m is influenced by constitutive host and parasite factors but independent of any parasite load-dependent responses. We sought to confirm associations with two constitutive host factors known to influence parasite growth: iron39 and complement factor H (FH)40,41 (Supplementary Dataset 1). Since we did not have premorbid blood samples we used convalescent blood as a proxy for pre-infection status, with samples collected 28 days after treatment when the host response was quiescent (median C-reactive protein 1.1mcg/mL (IQR 0.5-5.1, n=70), similar to healthy West African population levels42).

Iron deficiency is protective against malaria43 and reduces parasite multiplication in vitro. Consistent with this, levels of hepcidin (a regulator of iron metabolism and marker of iron sufficiency or deficiency44) were significantly correlated with m (rs=0.21, P=0.049) in 92 children who had not received blood transfusion.

FH is a constitutive negative regulator of complement activation which protects host cells from complement mediated lysis45 but many pathogens including P. falciparum have evolved FH binding proteins to benefit from this protection40,41. FH protects blood-stage parasites from complement mediated killing in vitro and higher plasma levels are associated with susceptibility and severity of malaria46. In the 14 children with residual day 28 plasma available, FH correlated with m (rs=0.75, P=0.002), providing further support that the quantitative estimates from our model exhibit expected relationships with known determinants of parasite growth.

We investigated whether we could identify any parasite processes associated with m, through correlation with parasite gene expression. Of 3704 parasite genes detected by RNA-Seq, adjusted for developmental stage distribution19, no individual genes passed the FDR adjusted P-value threshold of <0.05. Therefore we used weighted gene correlation network analysis to reduce dimensionality47, generating 17 modules of co-expressed parasite genes. Module eigengene values19 of two modules correlated with m (unadjusted Spearman correlation P<0.05); their hub-genes were PF3D7_1136000 (a conserved Plasmodium protein of unknown function) and PF3D7_1238300 (putative pre-mRNA-splicing factor CWC22). The PF3D7_1136000 module was negatively correlated (rs=-0.5, P=0.01) with m and contained 140 genes with greatest gene ontology enrichment in microtubule-based movement (Supplementary Tables 9 & 10). The PF3D7_1136000 module genes have high tolerance to insertional mutagenesis (Fig. 6a) and high parasite fitness following mutation (Fig. 6b), characteristics of winning mutants in competitive growth assays48, supporting the concept that lower expression of these genes may favour more rapid growth. 77 (55%) of the genes in this module exhibit greatest expression during gametocyte development49, consistent with the concept that increased sexual-stage commitment results in reduced asexual replication50. In contrast, the PF3D7_1238300 module was positively correlated with m (rs=0.46, P=0.03), and contained 45 genes enriched in translation functions (Supplementary Tables 9 & 10), plausible determinants of m, with mutagenesis tolerance typical of essential genes (Fig. 6a,b). Parasite genes differentially expressed between severe and uncomplicated malaria cases19 were highly over-represented in this module (16 of 45 (36%), P=1.2x10-8, Fisher exact test).

Figure 6

Parasite gene expression modules associated with predicted multiplication rate.

(a,b) Violin plots showing comparison of mutation insertion scores (a) and mutation fitness scores (b) between modules associated with predicted multiplication rate (PF3D7_1136000, n=138 genes; PF3D7_1238300 n=42 genes) and all other genes (n=3421). (Violin plots indicate distribution of data (kernel density estimates) and median (red circle); P for comparison between each module and all other genes using a two-sided Mann-Whitney test).

Discussion

Using a model-based approach to estimate the within-host dynamics of pathogen load and its determinants in human infection we could estimate the relative contributions of parasite multiplication and host response to parasite load measured at a single point in time, and we have used these predictions to identify mechanistic determinants of parasite load in malaria. Our approach is based on clearly defined assumptions, but as with any attempt to model complex biology, alternative approaches are possible. We cannot, at present, propagate uncertainty throughout the sequential stages of the model fitting, prediction of parameter estimates in individual subjects, and association of these parameter estimates with real variables. However, estimating the dynamics of parasite load allows us to make inferences about disease biology and mechanisms associated with PGI which could not have been made using only direct measurements. Our mechanistic validation suggests that the relative estimates of latent variables are accurate enough to be useful, providing proof-of -principle that pathogen load dynamics can be estimated in humans. This approach could be refined and expanded to identify additional genetic and serological determinants of pathogen load dynamics. The latter should be identified prospectively, since use of convalescent samples may introduce confounding.

Parasite load is only one of the factors associated with severe malaria and its interpretation is dependent on epidemiological context10,15,29. Variations in the host response, naturally acquired immunity, and the expression of P. falciparum erythrocyte membrane protein 1 (PfEMP1) variants are also important determinants of severity and of disease phenotype10,15. We have previously suggested that variation in the dynamics of parasite load may explain why cerebral malaria and severe anaemia occur with parasites expressing the same PfEMP1 variants10, and our model-based approach predicted that slower growth and longer duration of illness may distinguish severe anemia from cerebral malaria.

The importance of pathogen load and the dynamic nature of host-pathogen interactions are often omitted from studies of life-threatening infectious diseases in humans3. Much of our understanding of the host-pathogen interactions comes from comparisons between individuals at the point of clinical presentation, despite the fact that they may be at different stages in the dynamic process of infection. This can result in seemingly paradoxical observations such as high levels of TNF or low levels of platelets associated with severe malaria15,16, whilst evidence also indicates that TNF and platelets mediate defense against malaria parasites5,6,15–17. Considering the dynamic nature of the host-parasite interaction may explain these paradoxes and identify protective mechanisms more efficiently.

We identified several mechanisms which might be considered as prototypes for host-directed therapy in malaria. Inhibition of type-1 interferon or CXCL10 signalling with inhibitory antibodies or small molecules might be strategies to enhance control of parasite load. The therapeutic potentials of cathepsin G and MMP9 may be counterbalanced by risk of collateral tissue damage, but selective targeting of receptors on the erythrocyte surface may be a useful paradigm for both treatment and prevention of malaria.

Our approach could be applied to some other infectious diseases in which pathogen load can be measured and for which we do not have effective treatments, including emerging viral infections like Ebola, and possibly highly resistant bacterial pathogens, for which host-directed therapies may life-saving2.

Methods

Subjects and laboratory assays

We used data from all of the malariatherapy patients reported by Dietz et al.4 and from all 139 Gambian subjects reported in our previous studies11,51,52 who had all of the following data available: age, parasite biomass estimate, plasma TNF concentration, duration of illness and severity of illness. No subjects were excluded after this selection, and all available data was included in analyses, with the exception that one outlier was excluded from parasite gene expression analysis. As described previously11,51,52, Gambian children (<16 years old) were recruited with parental consent from three peri-urban health centres in the Greater Banjul region, from August 2007 through January 2011 as part of a study approved by the Gambia Government/MRC Laboratories Joint Ethics Committee, and the Ethics Committee of the London School of Hygiene and Tropical Medicine. P. falciparum malaria was defined by compatible clinical symptoms in the presence of ≥5000 asexual parasites/μL blood, and any children suspected or proven to have bacterial co-infection were excluded. Severe malaria was specifically defined by the presence of prostration (SM1) or any combination of three potentially overlapping syndromes (cerebral malaria (CM), severe anemia (SA, hemoglobin <5 g/dL), and hyperlactatemia (blood lactate >5 mmol/L) - collectively SM2)11,51–53. Clinical laboratory assays, measurements of plasma TNF and IL-10 by Luminex, measurements of gene expression by RT-PCR, and estimation of total parasite biomass from PfHRP2 ELISA have been previously described11,52. Subject-level data from these Gambian children is available as Supplementary Dataset 1.

Statistical analyses

Statistical analyses were undertaken using the R statistical software54 and GraphPad Prism (GraphPad Software, Inc.). Directly measured continuous variables were compared between groups using unpaired or paired student’s t-test (when normally distributed) and the Mann-Whitney or Wilcoxon matched pairs tests (when normal distribution could not be assumed), and ANOVA or Kruskal-Wallis test for comparison across multiple groups. Associations between measured variables and latent variables were assessed using generalised additive models (GAM55, with the R package “mgcv”); the generalised cross-validation score and explained variance were used to select the best GAM once all model terms had significant effects (P<0.05). It was not possible to propagate uncertainty estimates through all stages from model development, calibration to the Gambian data, and prediction of latent variables in individual subjects, and so statistical analyses of latent variable were undertaken using their predicted values without any measure of uncertainty, and using non-parametric methods. Correlations between predicted values of latent variables and measured variables were done using Spearman correlation.

All hypothesis tests were two-sided with alpha = 0.05 unless specifically stated otherwise. One-sided testing was only used when justified by small sample size and a strong a priori hypothesis for the direction of effect. We did not adjust for multiple hypothesis testing, except in the case of gene expression analyses where false-discovery rate was controlled using the Benjamini-Hochberg method. Dose-response curves were fitted using asymmetrical sigmoidal five-parameter logistic equation in GraphPad Prism.

Model relating parasite multiplication, host response and parasite load

A process-based, stochastic simulation model was devised to reproduce the clinical data collected from the Gambian children. This was achieved by combining the information in the Gambian data with a model describing the first wave of parasitemia in non-immune adults who were deliberately infected with P. falciparum malaria to treat neurosyphilis (“malariatherapy”)4. These malariatherapy data, from the pre-antibiotic era, are the main source of information on the within-host dynamics and between-host variation in the course of parasitemia in untreated malaria infections. The model of Dietz et al.4 was modified and extended in order to be applied to the Gambian data, and we made the assumption that the Gambian children presented to hospital prior to the first peak of parasitemia.

Model of ascending parasitemia in malariatherapy subjects

The model relates parasite density after each 2-day asexual blood stage cycle (P(t+2)) to the parasite density at the end of the previous cycle (P(t)) by the following equation: The host-specific parasite multiplication rate, m, is a property of both parasite and host, allowing for growth-inhibition by constitutive factors; the proportion of parasites that will survive the effects of the density-dependent host response in the present cycle is S: where P is the host-specific parasite load threshold at which the host response is strong enough to inhibit 50% of parasite growth in that cycle. Parasite growth inhibition (PGI(t)) is defined as 1-S(t).

Consistent with the original Dietz model, P(0) was set to 0.003 parasites/μl4.

The original Dietz model included an additional parameter, S, to help describe the decline in parasitemia after the peak of the first wave. S is the proportion of isogenic parasites surviving an additional density- and time-dependent host response, which might represent adaptive immunity (4). Estimates of the range of values of S in the Dietz dataset and model were used when simulating data but since this parameter has little influence on parasite densities prior to the peak it was not used to make subsequent predictions of m and P in individual Gambian subjects.

At the explicit request of Klaus Dietz and Louis Molineaux, we hereby communicate the following correction regarding their assertion that the malariatherapy patients had not received any treatment4: it was later found that 47 of these patients had indeed received subcurative treatment, and that those patients had significantly higher parasite densities. This is unlikely to influence our analysis, because treatment would only be provided when malariatherapy patients became very unwell, presumably at maximum parasitemia, whereas we assume that most patients with naturally acquired infection likely present prior to the peak parasitemia that might occur in the absence of treatment.

Fitting of the malariatherapy model to data from Gambian children

Individual-level parameter estimates for the malariatherapy dataset were kindly provided by Klaus Dietz. The logarithms of these 97 estimates of m and Pc were well described by a multivariate normal distribution, providing a quantitative description of inter-individual variation in the dynamics of the first wave of parasitemia. In order to use the Dietz model to simulate the Gambian data, a number of modifications and extensions were made. Some of these required estimation of additional parameters by comparing the model simulations with the Gambian data. Dietz et al. provided a statistical description of the parasite density at which first fever occurred (the “fever threshold”) in the form of the distribution of the ratio of threshold density to peak parasitemia. The median density at first fever was at 1.4% of peak density. We introduced the assumption that the onset of fever occurs at a particular threshold value of S, because fever is dependent on the production of cytokines like interleukin-6 and TNF, both components of the host response. This constitutes a process-based model for the onset of fever rather than a purely statistical one. Because individuals differ in their response to parasite load (captured through variation in P), this results in variation of parasite densities at first fever but ignores any potential variation among individuals with respect to magnitude of host response necessary to generate fever. The host response threshold for the onset of fever S = 0.86 was determined as the value of S calculated at 1.4% of the peak density of a simulated individual with the median parameter values. This yielded a distribution of fever ratios similar to the one described by Dietz et al.4, albeit with less variation.

To simulate the time between onset of fever and clinical presentation we made use of the self-reported duration of symptoms in the Gambian data. The model which was most consistent with these values assumed a gamma-distributed duration of symptoms in non-severe cases, and a possibility to present earlier in the case of more severe disease. Since parasite biomass is related to likelihood of having severe malaria11,12,56 the probability of early presentation on any day after onset of fever was set proportional to the (density-dependent) probability of having severe disease on that day. Scale (ζ) and shape (κ) parameters of the gamma distribution as well as the factor (ξ) for determining the probability of early presentation were estimated from the Gambian data.

We assumed that TNF production τ(t) increases monotonically with density dependent host response (1-S) and represented this relationship using a heuristic function of the form with free parameters a, b, λ* and γ estimated from the Gambian data.

The Gambian children had on average higher parasite densities than the malariatherapy patients, which led to a bad fit of the original model to the Gambian data. This was overcome by introducing the assumption that the Gambian children had a different range of values of P to the adult malariatherapy patients. A factor π was therefore estimated by which the ln P value from the Dietz model was multiplied. This led to overall higher parasite densities upon presentation. However, our model uses parasite biomass and its relationship with disease severity to predict day of presentation, and there is an interaction between the mean ln P and the variation in ln P, as well as the proportion of severe malaria in the simulated Gambian population. Based on the relatively low malaria transmission in the Banjul area of The Gambia, we assumed that severe cases (defined by the presence of any of: prostration, hyperlactatemia, severe anemia or cerebral malaria) were over-represented by hospital-based recruitment and that in an unselected population of children of similar age to those in our dataset only approximately 5% of all malaria infections would be severe28,29. Therefore we estimated a factor δ by which the variance of ln P should be multiplied such that both rate of severity as well as the distribution of parasite biomass matched well after fitting our simulation to the Gambian data.

The free parameters ζ, κ, ξ, a, b, λ*, γ, π and δ (Supplementary table 11), together summarized as θ, were estimated by fitting model simulations to the information on TNF, parasite density, and duration of symptoms, for any given candidate parameterization, a total of 139 clinically presenting individuals were simulated from the model, which corresponds to the size of the Gambian dataset. An objective function L(θ) was calculated, and a simulated annealing algorithm (provided by the “optim” function in R) determined the value for θ which maximizes this function. The log-likelihood L (θ) was comprised of three separate objectives. The first objective represented the log-probability that the frequency of severe cases in the simulation was equal to an assumed 5%, employing a binomial likelihood, given the actual number of severe cases sampled in 139 simulated individuals. The second objective considered the overlap between the bivariate distribution of ln parasite density vs. ln TNF concentration in the simulated data compared to the Gambian dataset. An approximate numerical value for this partial log-likelihood was obtained as the log probability of the Gambian data (density and TNF) given a two-dimensional kernel density estimate of the simulation output as a likelihood model. Kernel density estimates were obtained using the “kde2d” function in the “MASS” package in R. In this calculation, the TNF/density data points of severe or prostrated Gambian patients entered the partial likelihood with a weight of 1/11, to account for the oversampling of severe cases in the Gambian data. The third objective concerned the two-dimensional distribution of log density and duration since first fever. This partial log-likelihood was obtained using the same kernel-based approach described above, with weights of 1/11 for severe and prostrated cases. The overall log-likelihood L (θ) was calculated as a weighted sum of the three partial log-likelihoods, with the log-probability of having the desired true severity rate weighted with a factor of 68, which ensured similar magnitude of the three partial log-likelihoods at the optimum.

The results of the fitting algorithm were visually confirmed to yield a good overlap of the joint distributions of density and biomass, the duration of symptoms, TNF and biomass between simulation and the Gambian children. Approximate confidence intervals for the parameter estimates were determined by employing a 2nd degree polynomial to estimate the curvature of the maximum simulated likelihood surface in the vicinity of the parameter point estimate, assuming independence of parameters.

As in the original model of Dietz et al.4, peripheral parasite densities were used to determine the dynamic changes in parasitemia, implying a correlation between peripheral densities and total parasite biomass. Total parasite biomass per kg was calculated from the predicted parasite density by the equation 70,000 x 1.09 x predicted parasite density in parasites/μL, as has been determined previously for uncomplicated malaria cases in the Gambian dataset11.

Deterministic relationships between observable and latent variables

The range of values of m and ln P in a simulated population of 2000 patients were determined and each divided into 50 equally spaced increments in order to generate 2500 possible combinations of m and ln P for which all model outcomes were determined in order to visualize their relationships. For the purpose of this analysis, the time-dependent adaptive immune response parameters (which comprise S) were set for all subjects at their respective population median values. The model of Dietz et al. makes use of discrete 2 day time intervals4, corresponding to the duration of the intraerythrocytic cycle in a highly synchronised infection. However, naturally acquired infections are rarely this synchronous and the time since infection of our Gambian patients is an unknown continuous variable. In order to cope with this we assumed that the relationship between predicted outcome variables (parasite biomass, duration of illness and TNF concentration) and explanatory variables (m and P) could be approximated by smoothed GAM. We used the GAM to estimate values of m, P and parasite growth inhibition (PGI, 1-S) in the Gambian children, based on their known total parasite biomass, duration of symptoms and TNF concentration.

Predicting severe anemia and IL-10 concentrations

We used the data from the Gambian children to predict hemoglobin and IL-10 concentrations as continuous variables, using GLM with predicted P, predicted m, and age as explanatory variables. We then simulated a population of 50,000 1-year olds with malaria, allowing for normal variation in baseline hemoglobin concentration57, and adjusting P values according to a linear relationship between predicted ln P and age in the Gambian children. To predict the occurrence of severe anemia, we calculated the proportion of subjects estimated to have hemoglobin <5g/dL, and for these we calculated IL-10 concentrations as a continuous outcome.

RNA-sequencing and data analysis

We used RNA-sequencing data from all 24 subjects who were included in our previously reported study19 and had data to allow estimation of parasite growth inhibition and multiplication rate. RNA extraction, library preparation, sequencing and downstream analysis, including adjustment for leukocyte and parasite developmental stage mixture, have all been previously described19.

The association of gene expression with m and PGI was determined using a generalized linear model approach in edgeR, allowing adjustment for leukocyte and parasite developmental stage mixture. Coefficients and P-values were calculated for the relationships between adjusted log gene expression and PGI for all detected genes. False discovery rate (FDR) was then computed using the Benjamini-Hochberg approach and FDR below 0.05 was considered to be significant in the initial analysis. FDRs between 0 and 0.1 were considered to indicate consistent findings when comparing associations obtained under different model assumptions. Gene ontology (GO) terms were obtained from Bioconductor packages “org.Hs.eg.db” and “org.Pf.plasmo.db”. Fisher’s exact test was used to identify significantly over-represented GO terms from gene lists. The background gene sets consisted of all expressed genes detected in the data set. Enrichment analysis for biological process terms was carried out using the "goana()" function in edgeR. Using groups of genes significantly positively or negatively correlated with PGI, Ingenuity Pathway Analysis (Qiagen) was used to identify networks of genes functionally linked by regulators, interactions or downstream effects, which were visualized as radial plots centered around the most connected network member. The weighted gene co-expression network analysis (WGCNA) tool47 was used to construct modules of highly co-expressed parasite genes, based on analysis of 23 samples (sample HL_478 was removed as an outlier in parasite RNA-seq analysis) as described previously19. Module eigengene values for each subject were correlated with predicted m, using Spearman correlation.

Parasite culture, growth and invasion assays

P. falciparum 3D7 strain was used in continuous culture for all of the experiments unless otherwise stated. Asexual blood stage parasites were cultured in human blood group A red cells, obtained from the National Blood Service, at 1-5% hematocrit, 37°C, 5% CO2 and low oxygen (1% or 5%) as described previously58,59. Growth medium was RPMI-1640 (without L-glutamine, with HEPES) (Sigma) supplemented with 5 g/L Albumax II (Invitrogen), 147 μM hypoxanthine, 2 mM L-glutamine, and 10 mM D-glucose. Parasite developmental stage synchronization was performed using 5% D-sorbitol to obtain ring stage parasites or Percoll gradients for schizont stage enrichment58,60. For growth assays, schizonts were mixed at <1% parasitemia with uninfected erythrocytes at 2% final hematocrit. Cathepsin G (Abcam) or recombinant active MMP9 (Enzo) were added for 72 hour incubation to allow two replication cycles. Growth under each condition was calculated relative to the average growth in untreated samples. Invasion assays were performed by adding parasites synchronised at the schizont stage to target erythrocytes and incubating for 24 hours. Cathepsin G and MMP9 were either pre-incubated with the target cells overnight followed by four washes with RPMI to completely remove them, or they were added directly to the culture of schizonts with target erythrocytes for 24 hours. The same protocol was followed for other P. falciparum strains except Dd2, for which magnetic purification was used to purify schizonts61. For combined treatments, cathepsin G was added to target erythrocytes and MMP9 was added at the same time as schizonts.

Flow cytometry for parasitemia and erythrocyte surface receptor expression

Flow cytometry was performed using a BD LSR Fortessa machine and analysis was conducted using FlowJo v10 (TreeStar Inc.), and gating strategies are show in Supplementary Figure 5. To assess parasitemia, 1μl of sample at 50% hematocrit was stained with Hoechst 33342 (Sigma) and dihydroethidium (Sigma) and then fixed with 2% paraformaldehyde (PFA) before flow cytometry as previously described62. Erythrocyte surface receptor expression was assessed by median fluorescence intensity of erythrocytes labelled with monoclonal antibodies or by comparison with isotype control antibodies (Supplementary Table 12). Briefly, erythrocytes were washed twice before resuspending at 50% haematocrit, of which 1-2μl was stained in 100μl of antibody cocktail in FACS buffer (2% fetal bovine serum, 0.01% sodium azide in PBS) for 30 minutes in the dark on ice. Samples were washed twice in FACS buffer and then fixed in 300μl FACS buffer with 2% paraformaldehyde. Surface receptor loss was calculated from the difference between the treated and untreated sample median fluorescent intensities after the isotype control antibody fluorescence had been subtracted.

Whole blood stimulation and Cathepsin G and MMP9 ELISA

Whole blood was collected from 8 healthy adult donors and plated at 25% hematocrit, and incubated overnight with or without 1μM PMA (Sigma). Supernatant was collected to perform Cathepsin G (CTSG ELISA Kit-Human, Aviva Systems Biology) and MMP9 (Legend Max Human MMP-9, Biolegend) ELISAs, and erythrocytes were collected for assessment of surface receptor expression. The same ELISA kits were used to measure cathepsin G and MMP9 in acute (day 0) plasma samples from children with malaria.

C-reactive protein, Hepcidin, and complement Factor H ELISA

Using plasma samples collected 28 days after infection, CRP was measured using the human Simple Step ELISA kit (Abcam) and hepcidin concentration was measured in subjects who had not received blood transfusion using the Hepcidin-25 bioactive ELISA kit (DRG), both according to the manufacturer’s instructions, with duplicate measurements when sufficient plasma was available. Complement Factor H assays were performed using an in-house ELISA as described63.