Introduction

Plasmodium knowlesi is a zoonotic malaria parasite carried by simian reservoirs and transmitted to people through bites of infected Anopheles mosquitoes. Since the first report of hundreds of human P. knowlesi cases in Malaysian Borneo in 2004, reported P. knowlesi incidence has markedly increased across Southeast Asia1,2. Although increases in cases are partially due to improved detection with the application of molecular diagnostics, data indicate genuine increases in both the proportion and incidence of P. knowlesi3,4. In Malaysia, no indigenous cases of nonzoonotic malaria have been reported since 2018 and P. knowlesi is now the main cause of malaria in humans and can cause severe and fatal malaria when untreated5,6. These increases in reported incidence are strongly associated with deforestation, suggesting that landscape and land use change may be increasing contact between people, mosquitoes and macaque reservoirs of P. knowlesi7,8. P. knowlesi risk in humans is highly spatially heterogenous, with clusters of cases reported in households and villages9,10,11. This may reflect similar spatial variation in the force of infection from wildlife reservoir populations or increased adaptation to transmission between human hosts.

Many emerging pathogens have both zoonotic and nonzoonotic transmission pathways in which an initial spillover event from an animal reservoir may lead to stuttering chains of human-to-human transmission or larger outbreaks12. For P. knowlesi, genetic evidence suggests that transmission to humans remains primarily zoonotic, occurring through bites from vectors infected by wild non-human primate populations13. Evidence suggests P. knowlesi infections are typically chronic with low-level parasitemias in long-tailed macaques (Macaca fascicularis), the main reservoir14. The spatial distribution of human cases mirrors P. knowlesi infection rates in macaques, with the highest infection rates in humans and macaques reported in Malaysian Borneo15. P. knowlesi has a rapid 24 h replication rate with observed prepatent periods in humans ranging from 9 to 12 days16. Although limited data are available on the duration and frequency of human infectiousness to mosquitoes, the feasibility of human–mosquito–human transmission of P. knowlesi has been demonstrated once experimentally in a laboratory study16 but has not been formally described in natural settings. A recent systematic literature review concluded nonzoonotic P. knowlesi transmission was biologically plausible, but there was no empirical evidence that human-mosquito-human P. knowlesi transmission occurs naturally17. Notably, no previous studies had systematically assessed whether spatiotemporal patterns of P. knowlesi cases were consistent with nonzoonotic transmission.

Quantifying the relative contributions of zoonotic and nonzoonotic transmission is essential for designing effective surveillance and control. Transmission pathways are represented by two epidemiological parameters: the spillover rate, the rate at which a pathogen is transmitted from animals to humans, and the reproductive number (R), the number of secondary human cases resulting from an infectious human individual over the course of their infection18. Directly attributing human infections to zoonotic sources is challenging as the distributions and prevalence of wildlife reservoirs are largely unknown, and infection occurs indirectly via mosquito vectors. Alternatively, model-based inference methods can be used to infer likely transmission routes from routinely collected surveillance data on reported cases18,19,20,21. Diffusion network approaches have been used to quantify individual malaria case reproduction numbers (RC) by evaluating the probability of two cases being connected using the time of symptom onset, spatial locations of cases, and estimated serial interval (SI), i.e. the time from symptom onset in one case to symptom onset in the secondary case19,22,23,24. These methods also allow assessment of the likelihood of cases occurring from outside (unobserved) sources not represented in human surveillance data, such as introductions through zoonotic spillover18,21.

Here, we use a diffusion network approach within a Bayesian framework to estimate RC for P. knowlesi and nonzoonotic malaria species across Malaysia from 2012 to 2020. Due to an active malaria elimination programme, Malaysia has an exceptionally strong surveillance system, and detailed information is available for all cases. We use this surveillance dataset and parameters identified from the literature to estimate the SI for nonzoonotic P. knowlesi transmission. We then quantify estimates and uncertainty around RC for the zoonotic malaria P. knowlesi and nonzoonotic malaria P. falciparum and P. vivax. Models of P. knowlesi transmission were consistent with a large reservoir of unobserved infection sources, suggesting human cases are primarily driven by spillover and with limited subcritical onward transmission (RC < 1). In contrast, RC estimates for P. falciparum and P. vivax identified historical sustained human–mosquito–human transmission (i.e. RC > 1), with subsequent decreases in RC observed over time as the country approached elimination of transmission of these two nonzoonotic malaria parasites.

Results

Between 16 December 2011 and 3 January 2021, 32,635 malaria cases were reported to the Malaysian national malaria surveillance system. Of these, the majority (n = 26,093, 79.95%) were reported from the East Malaysian states of Sabah and Sarawak (Fig. 1). All suspected P. knowlesi cases were confirmed by molecular methods, with confirmed P. knowlesi accounting for over 70% (n = 23,143) of malaria cases nationally and no indigenous nonzoonotic malaria cases reported after 2018 (Fig. 2). Based on travel history to malaria-endemic areas within the past month and previous malaria diagnoses, cases are classified by the Malaysian Ministry of Health as indigenous, introduced, imported or reoccurrences (for P. vivax). Most malaria cases were classified as indigenous (n = 27,125/32,635, 83.12%), including over 99% of P. knowlesi cases. Over 80% (n = 18,691/23,143) of P. knowlesi cases were reported in men, with a median age of 38 years (interquartile range [IQR]: 26–49 years). Similarly, most nonzoonotic malaria cases were also reported in adult men; the median age was 31 (IQR: 21–43 years), and 84% (7945/9489) were male.

Fig. 1: Annual P. knowlesi incidence in Malaysia from 2012–2020.
figure 1

Annual parasite incidence (API) of P. knowlesi per district, API is presented as the annual number of malaria cases per 1000 people.

Full size image
Fig. 2: Monthy numbers of reported P. knowlesi cases and nonzoonotic malaria cases.
figure 2

Plots show the total numbers of reported cases for P. knowlesi and for nonzoonotic malaria species (P. falciparum, P. vivax, P. malariae, P. ovale) for East Malaysia (Sabah and Sarawak, top plot) and West Malaysia (bottom).

Full size image

Malaysia has high health system coverage with free treatment for malaria25. For available surveillance data, the median time between symptom onset and treatment was 5 days, with over 97% of individuals seeking treatment within 2 weeks (Fig. 3a). Out of 16,765 P. knowlesi cases screened for gametocytes by microscopy, 5.74% (95% confidence interval [CI]: 5.39–6.10%) had microscopically detectable gametocytes, consistent with other reports on the P. knowlesi gametocyte carriage in clinical cases26. Infectiousness to mosquitoes is primarily dependent on the density of malaria gametocytes, with infectious individuals predominantly having gametocyte densities high enough to be detected microscopically27,28. Proportions of infectious P. knowlesi cases increased with time since symptom onset; however, due to the relatively rapid reporting and treatment times, the mean duration of human infectiousness was only 1.645 days (range: 0–27 days) (Fig. 3b, Supplementary Information).

Fig. 3: Duration of infectiousness from surveillance data.
figure 3

The time of infectiousness was calculated from a time between symptom onset and treatment for P. knowlesi cases and b proportion and 95% CI of clinical P. knowlesi cases having microscopically observed gametocytes on presentation to health facility. Data are presented as the observed proportion of cases with microscopically observed gametocytes from 16,765 clinical P. knowlesi cases with the error bars representing 95% confidence intervals, the line is predicted probability; see Supplementary Fig. 3.

Full size image

We assessed the likelihood of cases being part of the same transmission chain based on species-specific SI parameters, dates and geographical locations of reported cases, and case classification as indigenous or non-indigenous (imported, introduced or reoccurrences). This modelling framework additionally included a threshold parameter (ε) that prevents the likelihood from becoming zero due to a complete break in the transmission chain and, therefore, epidemiologically represents the probability of infection from unobserved external sources needed to sustain transmission. We assessed ε ranges using mean prior values ranging from <0.001, a prior distribution centring around the assumption that nearly all cases are infected from human sources included in the surveillance data to 1, a prior distribution assuming that all human cases are infected from external zoonotic or imported sources (Supplementary Table 4)19. P. knowlesi models with ε values approaching 1 resulted in better fits, implying that all or nearly all human cases are infected from external sources, most likely to be zoonotic spillover from non-human primates in this context (Supplementary Tables 4 and 5). For each individual, their estimated number of secondary infections was calculated to give a time-varying point estimate of RC. An individual RC value < 1 indicates a low likelihood that the case is the source infection for any other case in the dataset. Mean RC values for P. knowlesi in East Malaysia from 2012 to 2020 were 0.0739 (95% CI: 0.0006–0.278) and 0.0147 (95% CI: 0.000003–0.0637) for West Malaysia. There were no estimates of RC greater than 1 in the best-fitting models, and no scenarios assessed that had mean RC values greater than 1, indicating no sustained transmission (Supplementary Information). As we estimated time-varying RC values for each individual, the estimates of all RC values as less than 1 also suggest there is no major overdispersion, excluding large super-spreading events. Although there were clear increasing trends in the numbers of reported P. knowlesi cases through the study period, this was not reflected in increases in RC estimates, and there was no evidence of deviation from zoonotic–mosquito–human transmission patterns.

In contrast, models for nonzoonotic malaria identified evidence for clear indigenous human-mosquito-human transmission chains, with 17% of P. falciparum cases and 24% of P. vivax cases with estimated RC values greater than 1 from 2012 to 2020. For P. falciparum, the mean RC value was 0.568 (95% CI: 0–1.97) for East Malaysia and 0.386 (95% CI: 0–2.04) for West Malaysia. The mean RC for P. vivax was 0.524 (95% CI: 0–2.436) and 0.230 (95%CI: 0, 1.820) for East and West Malaysia, respectively (Fig. 4). The probability of P. falciparum or P. vivax cases leading to onward transmission varied over space and time (Fig. 5). There were clear decreasing temporal trends in RC estimates for both nonzoonotic malaria species, with no RC values above 1 detected after 2018 when the last indigenous nonzoonotic malaria case was reported by the Malaysian Ministry of Health. Models clearly identified an outbreak of P. vivax reported in West Malaysia between 2016 and 2017, demonstrating the utility of these methods in pre-elimination settings29.

Fig. 4: Case reproductive numbers for P. falciparum, P. vivax and P. knowlesi.
figure 4

RC estimates per year for P. falciparum, P. vivax and P. knowlesi in East and West Malaysia. Data are presented as the median, first and third quartiles and range, with points more than 1.5 times the interquartile range plotted separately as outliers.

Full size image
Fig. 5: Mean posterior probability of RC (case reproductive number) estimates exceeding 1.
figure 5

The probability of RC exceeding 1 is shown for a P. falciparum and b P. vivax; P. knowlesi probabilities were not calculated as there were no estimates of RC greater than 1 in the best-fitting models, but village-level RC estimates can be found in Supplementary Fig. 3.

Full size image

Discussion

Quantifying rates of nonzoonotic transmission is critical to understanding the dynamics of emerging zoonotic diseases and identifying appropriate interventions. While many vector-borne and zoonotic diseases have emerged to become widespread in human populations (e.g. P. knowlesi, Lyme disease), the increase in human case numbers is not necessarily correlated with increased nonzoonotic transmission30,31. Here, we demonstrate the use of quantitative model-based approaches to identify key parameters and estimate individual reproductive numbers of the zoonotic and vector-borne disease P. knowlesi. Model results strongly suggest that P. knowlesi remains a primarily zoonotic disease driven by spillover, finding no evidence of sustained nonzoonotic P. knowlesi transmission. These results sharply contrast with RC estimates of the nonzoonotic malaria P. falciparum and P. vivax, where clear chains of human-mosquito–human transmission were identified. This highlights the utility of this method to assess the probability of nonzoonotic transmission using high-quality routine surveillance data and provides a critical tool to monitor changing transmission patterns.

A key determinant of the probability of nonzoonotic P. knowlesi transmission is whether humans are infectious. Using the largest dataset of gametocyte carriage rates in P. knowlesi clinical cases reported, consistent proportions of cases had microscopically observed gametocytes (5.74%), supporting the possibility of nonzoonotic P. knowlesi transmission. However, P. knowlesi gametocyte carriage rates in clinical cases upon presentation to health facilities are notably lower than for other nonzoonotic malaria species. For example, for P. vivax infections, which are biologically similar to P. knowlesi, up to 100% of clinical cases have microscopically observed gametocytes at presentation32. More extensive meta-analyses of gametocyte carriage in P. falciparum cases have also identified a wider range (12.1–32.2%) of microscopically observed gametocytes at rates typically higher than P. knowlesi cases reported within this dataset33. While this does suggest human P. knowlesi cases have the potential to be infectious, this parasite may be less adapted to human hosts in terms of gametocytogenesis. Critical questions remain about how infective P. knowlesi-infected humans are to mosquitoes, particularly considering the strong circadian rhythms of infectiousness to mosquitoes observed in P. knowlesi-infected non-human primates34.

Our model fits strongly suggests that P. knowlesi remains primarily driven by zoonotic spillover, finding no evidence of sustained nonzoonotic P. knowlesi transmission. For all of the best fitting models, there were no estimates of RC greater than 1, and in all scenarios assessed, there were no estimates suggesting RC is routinely over 1. There was additionally no clear evidence of changes to P. knowlesi RC over time that would suggest P. knowlesi may be adapting to transmission between humans. This contrasts with RC estimates for both P. falciparum and P. vivax, where RC estimates frequently exceed 1 and show clear decreasing temporal trends during an extensive malaria elimination campaign. This highlights the utility of these methods in monitoring changing transmission patterns of both emerging and endemic vector-borne diseases19,20.

While these results can shed insight into broad transmission patterns, this analysis is not without significant limitations. Model-based estimates alone cannot be used to conclusively prove that two cases are linked. This modelling approach estimated RC under different assumptions but did not reconstruct transmission networks linking cases. Further genetic, epidemiological and ecological data would be needed to prove transmission occurred between two individuals and rule out multiple spillover events occurring from zoonotic reservoirs. Additionally, while P. knowlesi results are consistent with a large reservoir of unobserved infections, this modelling framework does not explicitly attribute these to zoonotic sources. Other potential sources include asymptomatic or unreported infections. However, given the high coverage of the Malaysian surveillance system, infrequency and low parasite density of asymptomatic human P. knowlesi infections35,36 and widespread distribution of infected simian reservoirs15, unobserved infections would appear to be mostly from zoonotic sources. As higher rates of low-density, asymptomatic P. knowlesi infections have been reported in people residing near clinical cases, further studies are needed to identify whether peri-domestic or localised nonzoonotic P. knowlesi transmission is occurring but not detected by clinical surveillance systems35. Higher estimates of RC values could be used to guide targeted investigations at locations with the highest likelihood of nonzoonotic transmission. Future studies could additionally explore geographic variability in human–mosquito contact rates, human immunity, treatment-seeking behaviours and other factors impacting transmission probabilities and SI distributions.

Despite these limitations, results indicate there is no evidence of sustained nonzoonotic P. knowlesi transmission and demonstrate decreases in the transmission of P. falciparum and P. vivax in Malaysia throughout the past decade. This has important implications for the design of malaria surveillance and control programmes in areas with zoonotic malaria. Current measures, such as insecticide-treated net distribution and active case detection, have proven highly effective at controlling nonzoonotic malaria. While these measures may also be reducing nonzoonotic P. knowlesi transmission, new control strategies are urgently needed to mitigate the risks of spillover from wildlife hosts. With the inclusion of zoonotic malaria in global malaria elimination certification criteria, P. knowlesi poses a major threat to malaria elimination across Southeast Asia37. These methods provide a template for monitoring potential changes in transmission patterns of P. knowlesi and other zoonotic and vector-borne diseases.

Methods

Ethics

Ethical approval for this study was granted by the Medical Ethics and Research Committee of the Ministry of Health Malaysia (NMRR-21-470-58703) and the London School of Hygiene and Tropical Medicine Research Ethics Committee (Ref. 25654). The study protocol was registered under the National Medical Research Registry, Malaysia.

Surveillance data

National malaria surveillance data collected by the Malaysian Ministry of Health were assembled from 16 December 2011 to 3 January 2021. Data on malaria cases included demographic data (age, gender and occupation), case notification date to the Ministry of Health, symptom onset date, diagnosis date, method of detection, importation status, parasite species diagnosed, presence of gametocytes and health facility of detection. Symptom onset dates were excluded as unreliable if they occurred after diagnosis or notification dates or more than 30 days prior to diagnosis. Due to frequent misdiagnosis of P. knowlesi by microscopy, diagnosis of all presumptive P. knowlesi cases and imported malaria, cases were confirmed using species-specific molecular methods25. Additionally, surveillance data included geographical data on the case residence, probable location of the infection (state, district, village name and/or GPS coordinates) and travel history. All case locations were manually confirmed and geolocated to the nearest village centroid (Supplementary Information). Cases without accurate address data were excluded from the final analysis (n = 2117).

SI estimation

The SI, the time between symptom onset of primary and secondary malaria cases, requires estimating the duration of a series of sequential processes which need to occur in a human–mosquito–human transmission cycle38. While these intervals can be inferred from contract tracing data for directly transmitted diseases, SIs can only be estimated indirectly for vector-borne diseases39,40. SI durations have previously been estimated for P. falciparum and P. vivax19,40. However, there is a lack of empirical evidence on human–mosquito–human P. knowlesi transmission as this has only been observed once under experimental conditions16. The rapid replication cycle of P. knowlesi and weak evidence of adaptation to humans suggests the SI may differ from that of other nonzoonotic malaria parasites34.

To estimate the SI for human–mosquito–human P. knowlesi transmission, we adapted a quantitative model-based approach40. This models the SI as the sum of random variables representing sequential steps in the transmission cycle, including the prepatent period, human to mosquito transmission period, extrinsic incubation period, mosquito-to-human transmission period and infection-to-reporting periods. These were parameterised based on a combination of data from secondary literature and the Malaysian malaria surveillance dataset (Supplementary Information). Due to the typically low parasite densities and infrequency of asymptomatic human P. knowlesi infections in the region, we assumed asymptomatic infections did not contribute to transmission35,41. We also assumed individuals were not infectious after reporting due to Malaysia’s mandatory treatment and hospitalisation policy. Finally, no data were available on the duration or timing of infectiousness of P. knowlesi in humans, and the relatively short duration between symptoms onset and treatment within this population precluded estimation of natural recovery rates. While P. knowlesi gametocytes are highly synchronous and have strong circadian rhythms in macaques34,42, available evidence suggests this is not the case in human infections43.

Transmission likelihood and R C estimation

To estimate the joint likelihood of transmission between cases and a case having an unobserved source of infection, we used an adapted version of the NetRate algorithm to calculate RC for each case19,20,22,44. This uses a Bayesian framework to estimate the likelihood of two cases being within the transmission chain based on the geographical locations of cases, indigenous case status, time of symptom onset and SI parameters (Supplementary Information). Case reproduction numbers (RC) are then estimated by, for each case, summing the normalised likelihood of that case being the infector of each of the other cases. Non-indigenous cases and cases with a symptom onset date before the candidate case are not included as potential infectees. This approach additionally includes a provision for missing or unobserved sources of infections that can infect any observed individual with probability ε. These epsilon edges (ε) link cases to an external source when the likelihood of transmission between two cases is sufficiently low. Higher values of ε indicate a larger proportion of cases are assumed to be infected by external sources outside the surveillance data; in the case of P. knowlesi, this is mostly to be zoonotic reservoirs.

We fit separate models for P. knowlesi, P. falciparum and P. vivax for both East Malaysia (Borneo) and West Malaysia. From surveillance reports, most cases reported only local travel within the same village or subdistrict. To account for a high likelihood of local transmission, we used a fixed value for δ of 0.1, the parameter determining transmission likelihood based on Euclidean distance between cases (assessed in kilometres); this corresponds to most transmission events occurring within 10 km based on reported travel history (Supplementary Information). We fit shifted Rayleigh distributions to describe a prior distribution of possible SIs for P. knowlesi and nonzoonotic malaria. This distribution is defined by α, the shape of the SI distribution, and γ, the shifting parameter accounting for the minimum time between a bite from an infected mosquito and infectiousness in a human. For P. falciparum and P. vivax, we used previously published parameters to inform the SI distribution, with γ fixed at 15 days and normally distributed priors on α with a mean of 0.003 and a standard deviation of 0.1; this corresponds to a mean SI of 36 days (95% range: 21–60 days)20. For P. knowlesi, the γ was fixed at 14 days based on the modelled time to infectiousness fit from the Malaysian surveillance data (Supplementary Information). To account for uncertainty in the SI distribution, we used a normally distributed prior for α with a mean of 0.002 and standard deviation of 0.001, capturing the full range of plausible estimates of the P. knowlesi SI.

For each malaria species and location, surveillance data was input as a time-ordered series of malaria notification dates. As there was considerable uncertainty in the proportion of zoonotic transmission of P. knowlesi, we conducted a sensitivity analysis using normally distributed priors for ε with mean values of 0.01, 0.1 and 1. As ε values were bounded between 1e−11 and 1, prior estimates were truncated at 1, forming a half-normal estimate for higher values. Similar sensitivity analysis was conducted for P. falciparum and P. vivax, with mean values of priors for ε ranging from 0.0001 to 0.000001, reflecting the expected high sensitivity of malaria surveillance systems. We assessed model fit based on second-order Akaike’s Information Criteria (AIC) and visually inspected the posterior distributions of key parameters (γ, δ, ε) to select the best-fitting model. For all final models, we estimated RC by creating a matrix of potential infectors and infectees with the normalised likelihood of cases being connected. For each individual, the fractional number of secondary infections was calculated to give a time-varying point estimate of RC. Models were previously validated through simulation studies44, with code available at https://github.com/IzzyRou/spatial_rcs.

Spatial and temporal patterns of RC were explored for all species and locations. For P. falciparum and P. vivax, the probability of onward transmission (RC > 1) was assessed using statistical models with joint temporally and spatially structured random effects with Bayesian inference approximated using the Integrated Nested Laplace Approximation (INLA)45. We classified cases into dichotomous categories based on RC estimates. For each species, we fit a binomial model with a logit link to the number of cases with RC estimates > 1 out of the total cases per year detected within 5 km2 grid cells (Supplementary Information). The best-fitting model was selected using Deviance Information Criteria (DIC) and included a spatial effect and a temporally structured random walk model of order 246. The spatial effect was modelled as a Matern covariance function using an approximation to a stochastic partial differential equations approach. Posterior expectations and probabilities were estimated using 1000 samples. All data analysis was conducted in R statistical software v 4.1, and data were checked and cleaned using Quantum GIS v 3.24.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.