Emulator-based Bayesian optimization for efficient multi-objective calibration of an individual-based model of malaria ====================================================================================================================== * Theresa Reiker * Monica Golumbeanu * Andrew Shattock * Lydia Burgert * Thomas A. Smith * Sarah Filippi * Ewan Cameron * Melissa A. Penny ## ABSTRACT Individual-based models have become important tools in the global battle against infectious diseases, yet model complexity can make calibration to biological and epidemiological data challenging. We propose a using a Bayesian optimization framework employing Gaussian process or machine learning emulator functions to calibrate a complex malaria transmission simulator. We demonstrate our approach by optimizing over a high-dimensional parameter space with respect to a portfolio of multiple fitting objectives built from datasets capturing the natural history of malaria transmission and disease progression. Our approach quickly outperforms previous calibrations, yielding an improved final goodness of fit. Per-objective parameter importance and sensitivity diagnostics provided by our approach offer epidemiological insights and enhance trust in predictions through greater interpretability. ## INTRODUCTION ### Individual-based models of infectious diseases Over the last century, mathematical modelling has become an important tool to analyze and understand disease- and intervention-dynamics for many infectious diseases. Individual-based models (IBMs), where each person is simulated as an autonomous agent, are now widely used. These mathematical models capture heterogeneous characteristics and behaviors of individuals, and are often stochastic in nature. This bottom-up approach of simulating individuals and transmission events enables detailed, robust and realistic predictions on population epidemic trajectories as well as the impact of interventions such as vaccines or new drugs (*1, 2*). Going beyond simpler (compartmental) models to capture stochasticity and heterogeneity in populations, disease progression, and transmission, IBMs can additionally account for contact networks, individual care seeking behavior, immunity effects, or within-human dynamics (*1-3*). As such, well-developed IBMs provide opportunities for experimentation under relatively naturalistic conditions without expensive clinical or population studies. Prominent recent examples of the use of IBMs include assessing the benefit of travel restrictions during the Ebola outbreak 2014–2016 (*4*) and guiding the public health response to the Covid-19 pandemic in multiple countries (*5*). IBMs have also been applied to tuberculosis (*6*), influenza (*7*), dengue (*8*), and many other infectious diseases (*2*). Within the field of malaria, several IBMs have been developed over the last 15 years and have been used to support understanding disease and mosquito dynamics (*9-11*), predict the public health impact or carry out economic analyses of (new) interventions (*12-15*); and investigate drug resistance (*16*). Many have had wide-reaching impact, influencing WHO policy recommendations (*12, 17-19*) or strategies of national malaria control programs (*20*). ### Calibration caveats and the curse of dimensionality For model predictions to be meaningful, modelers need to ensure their models accurately capture abstractions of the real world. The potential complexity and realism of IBMs often come at the cost of long simulation times and potentially large numbers of input parameters, whose exact values are often unknown. Parameters may be unknown because they represent derived mathematical quantities that cannot be directly measured or require elaborate, costly experiments (for example shape parameters in decay functions (*21*)), because the data required to derive them in isolation is incomplete or accompanied by inherent biases, or because they interact with other parameters. Calibrating IBMs poses a complex high-dimensional optimization problem and thus algorithm-based calibration is required to find a parameter set that ensures realistic model behavior, capturing the biological and epidemiological relationships of interest. Local optima may exist in the potentially highly irregular, high-dimensional goodness-of-fit surface, making iterative, purely sampling-based algorithms (e.g. Particle Swarm Optimization or extensions of Newton-Raphson) inefficient and, in light of finite runtimes and computational resources, unlikely to find global optima. Additionally, the *curse of dimensionality* means the number of evaluations of the model scales exponentially with the number of dimensions (*22*). As an example, for the model discussed in this paper, a 23-dimensional parameter space at a sampling resolution of one sample per 10 percentile cell in each dimension, would yield 10*number of dimensions* = 1023 cells. This is larger than number of stars in the observable Universe (of order 1022 (*23*)). Furthermore, most calibrations are not towards one objective or dataset. For multi-objective fitting, each parameter set requires the evaluation of multiple outputs and thus multiple simulations to ensure that all outcomes of interest are captured (in the model discussed here epidemiological outcomes such as prevalence, incidence, or mortality patterns). In this study, we applied a new approach to calibrate a well-established and used IBM of malaria dynamics called *OpenMalaria*. Malaria IBMs in particular are often highly complex (e.g. containing multiple sub-modules and many parameters), consider a two-host system influenced by seasonal dynamics, and often account for multifaceted within-host dynamics. OpenMalaria features within-host parasite dynamics, the progression of clinical disease, development of immunity, individual care seeking behavior, vector dynamics and pharmaceutical and non-pharmaceutical antimalarial interventions at vector and human level ([https://github.com/SwissTPH/openmalaria.wiki.git](https://github.com/SwissTPH/openmalaria.wiki.git)) (*3, 21, 24*). Previously, the model was calibrated using an asynchronous genetic algorithm (GA) to fit 23 parameters to 11 objectives representing different epidemiological outcomes, including age-specific prevalence and incidence patterns, age-specific mortality rates and hospitalization rates (*3, 21, 24*) (see Supplementary Texts 1 and 2 for details on the calibration objectives and data). However, the sampling-based nature and sequential function evaluations of GAs can be too slow for high-dimensional problems in irregular spaces where only a limited number of function evaluations are possible and valleys of neutral or lower fitness may be difficult to cross (*25*), (*26*). Other solutions to fit similarly detailed IBMs of malaria employ a combination of directly extracting parameter values from the literature where information is available, and fitting the remainder using multi-stage, modular Bayesian Markov Chain Monte Carlo (MCMC)-based methods (*27-32*). For these models, multiple fitting objectives are often not addressed simultaneously. Rather, to our knowledge, most other malaria IBMs are divided into functional modules (such as the human transmissibility model, within-host parasite dynamics model, and the mosquito or vector model), which are assumed to be influenced by only a limited number of parameters each. The modules are then fit independently and in a sequential manner (*28-32*). Modular approaches reduce the dimensionality of the problem, allowing for the use of relatively straightforward MCMC algorithms. However, these struggle with efficiency in high dimensions as their Markovian nature requires many sequential function evaluations (104–107 even for simple models), driving up computing time and computational requirements (*33*). Additionally, whilst allowing for the generation of posterior probability distributions of the parameters (*31*), the modular nature makes sequential approaches generally unable to account for interdependencies between parameters assigned to different modules and how their co-variation may affect disease dynamics. ### Emulators and Bayesian Optimization Progress in recent years on numerical methods for supervised, regularized learning of smooth functions from discrete training data allows us to revisit calibration of detailed mathematical models using Bayesian methods for global optimization (*34*). Current state-of-the art calibration approaches for stochastic simulators are often based around Kennedy and O’Hagan’s approach (*35*) (KOH), where a posterior distribution for the calibration parameters is derived through a two-layer Bayesian approach involving cascade of surrogates (usually GPs) (*36*). A first GP is used to model the systematic deviation between the simulator and the real process it represents, while a second GP is used to emulate the simulator (*37*). However, this approach is computationally intense when scaling to high-dimensional input spaces and multi-objective optimization. A fully Bayesian KOH approach is likely computationally heavy (*37*) for the efficient calibration of detailed malaria simulators like OpenMalaria. Single-layer Bayesian optimization with Gaussian processes (GPs) on the other hand have gained popularity as an efficient approach to tackle expensive optimization problems, for example in hyperparameter search problems in machine learning (*38, 39*). Assuming that the parameter-solution space exhibits a modest degree of regularity, a prior distribution is defined over a computationally expensive objective function by the means of a light-weight probabilistic emulator such as a Gaussian process. The constructed emulator is sequentially refined by adaptively sampling the next training points based on acquisition functions derived from the posterior distribution. The trained emulator model is used to make predictions over the objective functions from the input space with minimum evaluation of the expensive *true* (simulator) function. Purely sampling-based iterative approaches (like genetic algorithms) are usually limited to drawing sparse random samples from proposals located nearby existing samples in the parameter space. In contrast, the use of predictive emulators permits exploration of the entire parameter space at higher resolution. This increases the chances of finding the true global optimum of the complex objective function in question and avoiding local optima. Here, we use a single-layer Bayesian optimization approach to solve the multidimensional, multi-objective calibration of OpenMalaria (Fig. 1). Employing this single-layer Bayesian approach further allows for the direct comparison to previous calibration attempts for OpenMalaria as the objective functions are retained. We prove the strength and versatility of our approach by optimizing its 23 input parameters using real-world data on 11 epidemiological outcomes in parallel. To emulate the solution space, we explore and compare two prior distributions, namely a GP emulator and a *superlearning* algorithm in form of a Gaussian process stacked generalization (GPSG) emulator. We first use a GP emulator to emulate the solution space. Whilst GP emulators provide flexibility whilst retaining relative simplicity (*39*) and have been used previously as priors in Bayesian optimization (*38*), stacked generalization algorithms have not. They provide a potentially attractive alternative as they have been shown to outperform GPs and other machine learning algorithms in capturing complex spaces (*14, 40*). The stacked generalization algorithm (*40*) builds on the idea of creating ensemble predictions from multiple learning algorithms (*level 0 learners*). The cross-validated predictions of the level 0 learners are incorporated into a general learning system (level 1 meta-learner). This allows for the combination of memory-efficient and probabilistic algorithms in order to reduce computational time, whilst retaining probabilistic elements required for adaptive sampling. Here, we showcase the efficiency and speed of the Bayesian optimization calibration scheme and propose a novel modus operandi to parameterize complex mathematical models that harvests recent computational developments and is scalable to high dimensions in multi-objective calibration. ![Fig. 1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F1.medium.gif) [Fig. 1.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F1) Fig. 1. Overview of model calibration approaches by Bayesian optimization using Gaussian process and machine learning emulators. **A. General Framework.** The input parameter space is initially sampled in a space-filling manner, generating the initial core parameter sets (initialization). For each candidate set, simulations are performed with the model, mirroring the studies that yielded the calibration data. The deviation between simulation results and data is assessed, yielding goodness of fit scores for each parameter set. An emulator (C or D) is trained to capture the relationship between parameter sets and goodness of fit and used to generate out-of-sample predictions. Based on these, the most promising additional parameter sets are chosen (adaptive sampling by means of an acquisition function), evaluated, and added to the training set of simulations. Training and adaptive sampling are repeated until the emulator converges and a decision on the parameter set yielding the best fit is made. **B. Acquisition Function**. The acquisition function is used to determine new parameter space locations. Thus, ***θ*** is a vector of input parameters (23-dimensional for the model described here) to be evaluated during adaptive sampling. It incorporates both predictive uncertainty of the emulator and proximity to the minimum. **C. Gaussian process emulator**. A heteroscedastic Gaussian process is used to generate predictions on the loss functions, ![Graphic][1], for each input parameter set ***θ*. D. Gaussian process stacked generalization emulator**. Three machine learning algorithms (level 0 learners: bilayer neural net, multivariate adaptive regression splines and random forest) are used to generate predictions on the individual objective loss functions ![Graphic][2] and ![Graphic][3] (collectively ![Graphic][4]) at locations ***θ***. These predictions are inputs to a heteroscedastic (level 1 learner) which is used to generate the stacked learner predictions ![Graphic][5] and derive predictions on the overall goodness of fit ![Graphic][6]. ## RESULTS ### Calibration workflow The developed model calibration workflow approach is summarized in Fig. 1A. In brief, goodness of fit scores were first derived for randomly generated, initial parameter sets. The goodness of fit scores were defined as a weighted sum of the loss functions for each of 11 fitting objectives. These span various epidemiological measures capturing the complexity and heterogeneity of the malaria transmission dynamics, including the age-prevalence and age-incidence relationships, and are informed by a multitude of observational studies (see methods and Supplementary Text 2). Next, GP and GPSG emulators were trained on the obtained set of scores and used to approximate the relationship between parameter sets and goodness of fit for each objective. After initial investigation of different machine learning algorithms, the GPSG was constructed using a bilayer neural net, multivariate adaptive regression splines and random forest as *level 0* learners and a heteroscedastic Gaussian process as *level 1* learner (Fig. 1C-D, see methods and supplement). Using a lower confidence bound acquisition function based on the emulators’ point and uncertainty predictions for proposed new candidate parameter sets, the most promising sets were chosen. These parameter sets were simulated and added to the database of simulations for the next iteration of the algorithm. At the next iteration, the emulators are re-trained on the new simulation database and re-evaluated (Fig. 1B). This iterative process of simulation, training and emulation was repeated until a memory limit of 1024GB was hit. Approximately 130,000 simulations were completed in up to this point. ### Algorithm performance by iteration and time and convergence Both emulators adequately captured the input-output relationship of the calculated loss-functions from the simulator, with better accuracy when close to minimal values of the weighted sum of the loss functions, *F* (Fig. 2A). This is sufficient as the aim of both emulators within the Bayesian optimization framework is to find minimal loss function values rather than an overall optimal predictive performance for all outcome values. Examples of truth vs predicted estimates on a 10% holdout set are provided in Fig. 2A (additional plots for all objectives can be found in Supplementary Fig.s S2-S5). A *satisfactory fit* of the simulator was previously defined by a loss function value of *F* = 73.2 (*21*). The *previous best* model fit derived using the GA had a weighted sum of the loss functions of *F* = 63.7 (*21*). *Satisfactory fit* was achieved by our approach in the first iteration of the GPSG-based Bayesian optimization algorithm (GPSG-BO), and after six iterations for the GP-based algorithm (GP-BO) (Fig. 2B). The *current best* fit was approximately retrieved after six iterations for the GPSG-BO algorithm and after nine iterations for GP-BO, and was improved by both algorithms after ten iterations (returning final values *F* = 58.3 for GP-BO and 59.6 for GPSG-BO). This shows that the Bayesian optimization approach with either of our emulators very quickly achieves a better simulator fit than obtained with a classical GA approach that was previously employed to calibrate OpenMalaria. Of the two emulators, the GP approach finds a parameter set associated with a better overall accuracy and the GPSG reaches *satisfactory* values faster (both in terms of iterations and time). A likely explanation for this is that the GPSG-BO is unable to propagate its full predictive variance into the acquisition function. Only uncertainty stemming from the level 1 probabilistic learner (GP) is therefore captured in the final prediction. This leads to underestimation of the full predictive variance, and a bias towards exploitation in the early stages of the GPSG-BO algorithm (as illustrated by early narrow sampling, see Supplementary Fig.s S6-7). ![Fig. 2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F2.medium.gif) [Fig. 2.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F2) Fig. 2. Emulator performance. **A. Example of emulator predictions vs true values on a 10% holdout set.** Predictions are shown for the final iteration of each optimization (iteration 30 for GP-BO and iteration 23 for GPSF-BO). Here, emulator performances are shown for objective 4 (the age-dependent multiplicity of infection,*f*4) and the weighted sum *F*. Plots for all other objectives are provided in the supplement. GP = Gaussian process emulator, GPSG = Gaussian process stacked generalization emulator **B. Convergence**. Weighted sum of loss functions over 11 objectives associated with the current best fit parameter set by time in seconds. Satisfactory fit of OpenMalaria refers to a weighted sum of loss functions value of 73.2 (as defined previously (*21*)). Previous best fit for OpenMalaria was achieved by the genetic algorithm had a loss function value of 63.7. Our new approach yields a fit of 58.2 for GP-BO in iteration 21 within in 1.02*6 seconds (∼12 days) and 59.6 for GPSG-BO in iteration 10 in 6.00e5 seconds (∼7 days). GP-BO = Gaussian process emulator Bayesian optimization, GPSG-BO = Gaussian process stacked generalization emulator Bayesian optimization). **C. Example log prior parameter distributions and posterior estimates**. The most influential parameters on the weighted sum of the loss functions are shown here (see Fig. 3C). All other plots can be found in the supplement. The posterior estimates for GP-BO and GPSG-BO are shown in relation to those previously derived through optimization using a genetic algorithm (GA-O) Fig. 2C shows examples of the posterior estimates returned by the optimization algorithms in context of the log prior distributions for the parameters with the greatest effects on *F* (see also Fig. 3C). All algorithms return parameter values within the same range and (apart from parameter 4), clearly distinct from the prior mean. The fact that highly similar parameter values are identified by multiple algorithms strengthens confidence in the final parameter sets yielded by the algorithms. ![Fig. 3.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F3.medium.gif) [Fig. 3.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F3) Fig. 3. **A. Multiplicity of infection by age.** Comparison of simulator goodness of fit for objective 4, the age-specific multiplicity of infection (number of genetically distinct parasite strains concurrently present in one host). Simulations were carried out for the same random seed for all parameterizations and for a population size of N=5,000. **B. Simulated epidemiological relationship between transmission intensity (entomological inoculation rate, EIR) and P. falciparum prevalence (PfPR****2-10****)**. Simulated epidemiological relationship between the transmission intensity (EIR in number of infectious bites per person per year) and infection prevalence in individuals aged 2-10 years (PfPR2-10) under the parameterizations achieved by the different optimization algorithms. Lines show the mean across 100 random seed simulations for a simulated population size N=10,000 and the shaded area shows the 95% confidence interval. **C. Parameter effects on the objective variance**. Using the GP emulator, a global sensitivity analysis (Sobol analysis) was conducted. The tile shading shows the total effect indices for all objective functions and parameters grouped by function. SEN= Senegal, TZN = Tanzania. ### Optimal Goodness of Fit The best fit parameter sets yielded by our approach are provided in the supplement (Table S2). Importantly, after ten iterations of the GPSG-BO algorithm (approximately 7 days), and 20 iterations for the GP-BO algorithm (approximately 12 days), both approaches yielded similar values of the 11 objective loss functions, along with similar weighted total loss function values, and qualitatively similar visual fits and predicted trends to the data (Fig. 3A-B and supplement). We found this to be an unexpectedly fast result of the two algorithms. Details of the algorithm’s best fits to the disease and epidemiological data are shown in Supplementary Fig.s S8-S18. Overall, several objectives had visual and reduced loss-function improvements, for example to the objective on the multiplicity of infection (Fig. 3A). ### Impact / Parameter sensitivity analysis & External validation An additional benefit of using emulators is the ability to understand the outcome’s dependence on and sensitivity to the input parameters. To identify the most influential parameters for each of the 11 fitting objectives, we used the GP emulator trained on all available training simulation results from the optimization process (R2=0.53 [objective 7] - 0.92 [objective 3]) to conduct a global sensitivity analysis by variance decomposition (here via Sobol analysis (*41*)). Fig. 3C shows Sobol total effect indices quantifying the importance of individual parameters and describing each parameter’s contributions to the outcome variance for each objective. Our results indicate that most objectives are influenced by multiple parameters from different groups, albeit to varying degrees, thus highlighting the importance of simultaneous multi-objective fitting. Clusters of influential parameters can be observed for most objectives; for example, parameters associated with incidence of acute disease influence clinical incidence and pyrogenic threshold objectives. Some parameters have strong influence on multiple objectives, such as parameter 4, the critical value of cumulative number of infections and influences immunity acquisition; and parameter 10, a factor required to determine the pyrogenic threshold, which we find to be a key parameter determining infections progressing to clinical illness. ### Algorithm validation In order to test if our algorithms can recover a known solution, the final parameter sets for both approaches were used to generate synthetic field data sets, and our approaches were subsequently applied to recover the known parameter set. For the GP, 13 of the 23 parameters were recovered (Supplementary Fig. S19A). Those not recovered largely represented parameters to which the weighted loss function was found to be insensitive (Fig. 3C). Thus, rather than showing a shortcoming of the calibration algorithm, this suggests a potential for dimensionality reduction of the simulator and re-evaluation of its structure. ### Comparison of key epidemiological relationships and implications for predictions The new parameterizations for OpenMalaria were further explored to assess key epidemiological relationships, in an approach similar multiple-model comparison in Penny et al. 2016 (*12*). We examined incidence and prevalence of disease, as well as incidence of mortality for multiple archetypical settings, considering a range of perennial and seasonal transmission intensity and patterns. The results are presented in Fig. 3B and Supplementary Fig.s S20-30. The new parameterizations result in increased predicted incidence of severe episodes and decreased prevalence for all transmission intensities (thus also slightly modifying the prevalence-incidence relationship). While we found that the overall implications for the other simulated epidemiological relationship were small, the differences in predictions for severe disease may carry important implications for public health decision making. We conclude that our new parameterizations do not fundamentally bring into question previous research conducted using OpenMalaria, but we do suggest re-evaluation of adverse downstream events such as severe disease and mortality. ## DISCUSSION Calibrating individual-based models can be challenging as many techniques struggle with high dimensionality, or become infeasible with long model simulation times and multiple calibration objectives. However, ensuring adequate model fit to key data is vital, as this impacts the weighting, we should give model predictions in the public health decision making process. The Bayesian optimization approaches presented here provide fast solutions to calibrating individual-based models while improving model accuracy, and by extension prediction accuracy. Using a Bayesian optimization approach, we calibrated a detailed simulator of malaria transmission and epidemiology dynamics with 23 input parameters simultaneously to 11 epidemiological outcomes, including age-incidence and -prevalence patterns. The use of a probabilistic emulator to predict goodness-of-fit, rather than conducting sparse sampling, allows for cheap evaluation of the simulator at many locations and increases our confidence that the final parameter set represents a global optimum. Our approach provides a fast calibration whilst also providing a better fit compared to the previous parameterization. We are further able to define formal endpoints to assess calibration alongside *visual confirmation* of goodness of fit (*21, 28*), such as the emulator’s predictive variance approaching the observed simulator variance. The emulator’s ability to quantify the input stochasticity of the simulator also enables simulation at small population sizes, contributing to fast overall computation times. Despite the demonstrated strong performance of stacked generalization in other contexts such as geospatial mapping (*14, 40, 42-45*), we found that using a *superlearning* emulator for Bayesian optimization was not superior to traditional GP-based methods. In our context using GPSG sped up convergence of the algorithm, but both approaches, GP and GPSG, led to equally good fits. Each approach does however, have different properties with context-dependent benefits: The dimensionality reduction provided by GPSG approaches may lead to computational savings depending on the *level 0* and *level 1* learners. At the same time, only level 1 learner uncertainty is propagated into the final’ ’ predictions, which affects the efficacy of adaptive sampling and may lead to overly exploitative behavior, where sampling close to the point estimate of the predicted optimum is overemphasized, rather than exploring the entire parameter space (see supplements S2 and S3 on selected points). On the other hand, exploration/exploitation trade-offs for traditional GP-BO algorithms have long been examined and *no regret* solutions have been developed (*46*). The methodology presented here constitutes a highly flexible framework for individual based model calibration and aligns with the recent literature on using emulation in combination with stochastic computer simulation experiments of infectious diseases (*47*). Both algorithms can be applied to other parameterization and optimization problems in disease modelling and also in other modelling fields, such as physical or mobility and transport models. Furthermore, in the GPSG approach, additional or alternative level 0 can be easily incorporated. Possible extensions to our approach include combination with methods to adaptively reduce the input space for constrained optimization problems (*48*), or other emulators may be chosen depending on the application. For example, homoscedastic GPs, which are faster than the heteroscedastic approach presented here, may be sufficient for many applications (but not for our IBM in which heteroscedastic was required due to the stochastic nature of the model). Alternatively, the computational power required by neural net algorithms scales only linearly (compared with a nominal cubic scaling for GPs) with the sample size, and we envisage wide applications for neural net-based Bayesian optimization in high dimensions. In our example, the bilayer neural net algorithm completed training and prediction within seconds whilst maintaining very high predictive performance. Unfortunately, estimating the uncertainty required for good acquisition functions is difficult in neural networks, but solutions are being developed (*39, 49*). These promising approaches should be explored as they become more widely available in high-level programming languages. With the increased availability of code libraries and algorithms, Bayesian optimization with a range of emulators is also becoming easier to implement. The probabilistic, emulator-based calibration approach is accompanied by many benefits, including relatively quick global sensitivity analysis. As explored in this work, GP-based methods are easily coupled with sensitivity analyses, which provide detailed insights into a model’s structural dependencies and the sensitivity of its goodness of fit to the input parameters. To the best of our knowledge, no other individual-based model calibration study has addressed this. In the case of malaria models, we have shown the interdependence of all OpenMalaria model components and a relative lack of modularity. In particular, within-host immunity-related parameters were shown to influence all fitting objectives, including downstream events such as severe disease and mortality when an infection progresses to clinical disease. Thus, calibrating within-host immunity in the absence of key epidemiology and population outcomes can lead to suboptimal calibration and ultimate failure of the model to adequately capture disease biology and epidemiology. We have employed a different approach to calibrating OpenMalaria compared with previous methods but reach broadly similar comparisons to the natural history of disease. We also attainted a slightly improved but similar goodness of fit, the main benefit being improved fitting times and the ability to measure parameter importance. Given the high number of influential parameters for each epidemiological objective in our parameter importance investigations, and the overlap between parameter-objective associations, we argue that, where possible, multi-objective fitting should be preferred over purely sequential approaches. Our approach confirms that using a parallel approach to parameterization rather than a modular, sequential, one captures the joint effects of all parameters and ensures that all outcomes are simultaneously accounted for. To the best of our knowledge, no model of malaria transmission of comparable complexity and a comparable number of fitting objectives was simultaneously calibrated to all its fitting objectives. Disregarding the joint influence of *all* parameters on the simulated outcomes may negatively impact the accuracy of model predictions, in particular on policy-relevant outcomes of severe disease and mortality. Despite providing relatively fast calibration towards a better fitting parameter set, several limitations remain in our work. We have not systematically tested that a global optimum has been reached in our new approach, but assume it is close to a global minimum for the current loss-functions defined, as further iterations did not yield changes, and both the GP and GPSG achieved similar weighted loss function and parameter sets. We aimed to improve the algorithm to calibrate detailed IBM, but we did not incorporate new data, which will be important moving forward as our parameter importance and validation analysis highlights several key epidemiological outcomes on severe disease and mortality are sensitive to results. The key limitations of Bayesian optimization, particularly when using a Gaussian process emulator, are the high computational requirements in terms of memory and parallel computing nodes due to increasing runtimes and cubically scaling memory requirements of GPs. For this reason, we opted to not employ fully Bayesian KOH methods, which would double the number of GPs that would need to be run. Yet, memory limits may be reached before the predictive variance approached its limit. Furthermore, we chose an acquisition function with high probability to be *no regret* (*46*), but this likely overemphasizes exploration in the early stages of the algorithm considering the dimensionality of the problem and finite runtime. We opted here for pure exploitation every 5 iterations, but a more formal optimization of the acquisition function should be explored. The GPSG approach presented here can partially alleviate this challenge, depending on the choice of learning algorithms, but the iterative nature and need for many simulations remain. Memory- and time-saving extensions are thus worth exploring, such as incorporating GPU computing or adaptively constraining the prior parameter space, dimensionality reduction, or addressing alternative acquisition functions. Additionally, as with all calibration methodologies, many choices are left to the user, such as the size of the initial set of simulations, the number of points added per iteration, or the number of replicates simulated at each location. There is no general solution to this as the optimal choices are highly dependent on the problem at hand, and we did not aim to optimize these. Performance might be optimized further through a formal analysis of all these variables, however the methodology here is already fast, effective, and highly generalizable to different types of simulation models and associated optimization problems. Improving the loss-functions or employing alternative *Pareto front* efficiency algorithms was not the focus of our current study but would be a natural extension of our work, as would be alternative approaches to the weighting of objectives, which remains a subjective component of multi-objective optimization problems (*50*). A model’s calibration to known input data forms the backbone of its predictions. The workflow presented here provides great advances in the calibration of detailed mathematical models of infectious diseases such as IBMs. Provided sufficient calibration data to determine goodness-of-fit, our approach is easily adaptable to any agent-based model and could become the new modus operandi for multi-objective, high-dimensional calibration of stochastic disease simulators. ## METHODS ### Preparation of calibration data and simulation experiments Disease transmission models generally have two types of parameter inputs: core parameters, inherent to the disease and determining how its natural history is captured, and simulation options characterizing the specific setting and the interventions in place (Fig. 1A in the main manuscript). The simulation options specify the simulation context such as population demographics, transmission intensity, seasonality patterns and interventions and typically vary depending on the simulation experiment. In contrast, the core parameters determine how its epidemiology and aetiopathogenesis are captured. These include parameters for the description of immunity (e.g. decay of maternal protection), or for defining clinical severe episodes (e.g. parasitemia threshold). To inform the estimation of core parameters, epidemiological data on the natural history of malaria extracted from published literature and collated in previous calibrations of OpenMalaria(*3, 21, 24*) were re-used in this calibration round. These include demographic data such as age-stratified numbers of host individuals which are used to derive a range of epidemiological outcomes such as age-specific prevalence and incidence patterns, mortality rates and hospitalization rates. Site-specific OpenMalaria simulations were prepared, representing the studies that yielded these epidemiological data in terms of transmission intensity, seasonal patterns, vector species, intervention history, case management, and diagnostics (*24*). The mirroring of field study characteristics in the simulation options ensured that any deviation between simulation outputs and data could be attributed to the core parameters. Age-stratified simulation outputs to match to the data include numbers of host individuals, patent infections, and administered treatments. A summary of the data is provided in the Supplementary text 2. ### General Bayesian Optimization framework with emulators In our proposed Bayesian optimization framework (Fig. 1) we evaluated the deviation between simulation outputs and the epidemiological data by training probabilistic emulator functions that approximate the relationship between core parameter sets and goodness of fit. To test the optimization approach in this study we considered the original goodness of fit metrics for OpenMalaria detailed in (*21*) and in Supplementary Text 2, which uses either Residual Sum of Squares (RSS) or negative log-likelihood functions depending on the epidemiological data for each objective (*21, 24*). The objective function to be optimized is a weighted sum of the individual objectives’ loss functions. We adopted a Bayesian optimization framework where a probabilistic emulator function is constructed to make predictions over the loss functions for each objective from the input space, with a minimum amount of evaluations of the (computationally expensive) simulator. We compared two emulation approaches. Firstly, a heteroskedastic Gaussian process (GP) emulator and secondly a stacked generalization emulator (*40*). For approach 1 (*GP-BO*), we fitted a heteroskedastic Gaussian process with the input noise modelled as another Gaussian process (*51*) with a Matérn 5/2 kernel to account for the high variability in the parameter space (Fig. 1C) (*38, 52*). For approach 2 (*GPSG-BO*), we selected a two-layer neural network (*53-55*), multivariate adaptive regression splines (*56*), and a random forest algorithm (*57, 58*) as level 0 learners. With each iteration of the algorithm, the training was extended using adaptive sampling based on an acquisition function (*lower confidence bound*) that accounts for uncertainty and predicted proximity to the optimum of proposed locations (Fig. 1B). As the emulator performance improves (as assessed by its predictive performance on the test set) we gain confidence in the currently predicted optimum. ### Malaria transmission and disease simulator We applied our novel calibration approach to OpenMalaria ([https://github.com/SwissTPH/openmalaria](https://github.com/SwissTPH/openmalaria)), an open-source modelling platform of malaria epidemiology and control. It features several related individual-based stochastic models of *P. falciparum* malaria transmission and control. Overall, the OpenMalaria IBM consists of a model of malaria in humans linked to a model of malaria in mosquitoes and accounts for individual level heterogeneity in humans (in exposure, immunity, and clinical progression) as well as aspects of vector ecology (e.g. seasonality and the mosquito feeding cycle). Stochasticity is featured by including between- and within-host stochastic variation in parasite densities with downstream effects on immunity (*24*). OpenMalaria further includes aspects of the health system context (e.g. treatment seeking behavior and standard of care) (*3, 24*) with additional probabilistic elements such as treatment seeking probabilities or the option for stochastic results of diagnostic tests. An ensemble of OpenMalaria model alternative variants is available defined by different assumptions about immunity decay, within-host dynamics, heterogeneity of transmission, along with more detailed sub-models that track parasite genetics, and pharmacokinetic and pharmacodynamics. The models allow for the simulation of interventions, such as the distribution of insecticide treated nets (ITNs), vaccines, or reactive case detection (*59, 60*), in comparatively realistic settings. Full details of the model and the history of calibration can be found in the original publications (*3, 21, 24*) and are summarized in Supplementary Texts 1 and 2. In our application, we use the term *simulator* to refer to the OpenMalaria base model variant (*21*). ### Calibrating OpenMalaria: loss functions and general approach #### Aim Let ***f*** (***θ***) denote a vector of loss functions obtained by calculating the goodness of fit between simulation outputs and the real data (full details of loss function can be found in supplementary Text 2). In order to ensure a good fit of the model, we aim to find the parameter set ***θ*** that achieves the minimum of the weighted sum of 11 loss functions (corresponding to the 10 fitting objectives) ![Graphic][7], where *f**i*(***θ***) is the value of objective function *i* at ***θ*** and *w**i* is the weight assigned to objective function *i*: ![Formula][8] The weights are kept consistent with previous rounds of calibration and chosen such that different epidemiological quantities contributed approximately equally to *F*(***θ***) (see Supplementary Text 2). #### Step 1: Initialization Let *D* = 23 denote the number of dimensions of the input parameter space **Θ** and *W* = 11 the number of objective functions *f**i* (*θ*) *i* = 1, …, 11. Prior distributions consistent with previous fitting runs (*21*) were placed on the input parameters. As each parameter is measured in different units, we sampled from the *D* -dimensional unit cube **Θ** and converted these to quantiles of the prior distributions (*21*) (Supplementary Text 2and Supplementary Fig. S6). Previous research suggests that in high-dimensional spaces quasi-Monte Carlo (qMC) sampling outperforms random or Latin Hypercube designs for most function types and leads to faster rates of convergence (*61, 62*). We therefore used Sobol sequences to sample 1,000 initial locations from **Θ**. The GP can account for input stochasticity of the simulator. For each sample, we simulated 2 random seeds at a population size of 10,000 individuals. Additionally, 100 simulations were run at the centroid location of the unit cube to gain information on the simulator noise. Using small noisy simulations with small populations speeds up the fitting as the noisy simulations are less computational expensive than larger population runs. Replicates were used to detect signals in noisy settings and estimate the pure simulation variance (*51*). The 2000 unique locations were randomly split into a training set (90%) and a test set (10%). All simulator realizations at the centroid were added to the training set. #### Step 2: Emulation ##### 2.1: Emulator Training Each emulator type for each objective function was trained in parallel to learn the relationships between the normalized input space **Θ**, and the log-transform of the objective functions *f* (***θ***). In each dimension *d* ∈ *D*, the mean *ε**d* and standard deviation *σ**d* of the training set were recorded, *d* = 1, …, 23. ##### 2.2 Posterior prediction We randomly sampled 500,000 test locations in **Θ** from a multivariate normal distribution with mean ***θ******opt*** and covariance matrix **Σ**, where ***θ******opt*** is the location of the current best location and **Σ** is determined based on previously all sampled locations, and scaled each dimension to mean *ε**d* and standard deviation *σ**d*. The trained emulators were used to make predictions ![Graphic][9] of the objective functions ***F*** (***θ***) at the test locations. Mean estimates, standard deviations, and nugget terms were recorded. The full predictive variance at each location ***θ*** ∈ **Θ** corresponds to the sum of the standard deviation and nugget terms. From this, we derived the weighted sum ![Graphic][10], using weights *w* consistent with previous fitting runs (Smith 2012) with greater weighting for further downstream objectives. The predicted weighted loss function at location ***θ*** was denoted ![Graphic][11] with a predicted mean ![Graphic][12] and variance ![Graphic][13]. Every 15 iterations, we increase the test location sample size to 5 Million to achieve denser predictions. #### Step 3: Acquisition We chose the lower confidence bound (LCB) acquisition function to guide the search of the global minimum (*63*). Lower acquisition corresponds to *potentially* low values of the weighted objective function, either because of a low mean prediction value or large uncertainty (*64*). From the prediction set at iteration *t*, we sample without replacement 250 new locations ![Graphic][14], with the hyperparameter *v* = 1 and ![Graphic][15], where *T**t* is the number of previous unique realisations of the simulator at iteration *t*, and *δ* = 0.01 is a hyperparameter (*46*). We choose this method as with high probability it is *no regret (46, 64)*. With increasing iterations, confidence bound-based methods naturally transition from mainly exploration to exploitation of the current estimated minimum. In addition to this, we force exploitation every 10 iterations by setting *τ**t* = 0). #### Step 4: Simulate The simulator was evaluated at locations identified in step 3 and the realisations were added to the training set. Steps 2-4 were run iteratively. The Euclidian distance between locations of current best realisations was recorded. #### Step 5: Convergence Convergence was defined as no improvement in the best realisation, argmin***F******F***. ### Emulator definition We compared two emulation approaches. Firstly, a heteroskedastic GP emulator and secondly a stacked generalization emulator (*40*) using a two-layer neural net, multivariate adaptive regression splines (MARS) and a random forest as level 0 learners and a heteroskedastic GP as level 1 learner: #### 1.1.1 Heteroskedastic Gaussian Process (hetGP) (*65*) We fitted a Gaussian process with the input noise modelled as another Gaussian process (*51*). After initial exploration of different kernels, we chose a Matérn 5/2 kernel to account for the high variability in the parameter space. A Matérn 3/2 correlation function was also tested performed equally. Each time the model was built (for each objective at each iteration), its likelihood was compared to that of a homoscedastic Gaussian process and the latter was chosen if its likelihood was higher. This resulted in a highly flexible approach, choosing the best option for the current task. #### 1.1.2 Gaussian Process Stacked Generalization (GPSG) Stacked generalization was first proposed by Wolpert 1992 (*40*) and builds on the idea of creating ensemble predictions from multiple learning algorithms (level 0 learners). In *superlearning*, the cross-validated predictions of the level 0 learners are fed into a level 1 meta-learner. We compared the 10-fold cross-validated predictive performance of twelve machine learning algorithms on the test set. All algorithms were accessed through the mlr package in R (*66*). We compared two neural network algorithms (brnn (*54*) for a two layer neural network and nnet for a single-hidden-layer neural network (*67*)), five regression algorithms (cvglmnet (*68*) for a generalised linear model with LASSO or Elasticnet Regularization and 10-fold cross validated lambda, glmboost (*69*) for a boosted generalized linear model, glmnet (*68*) for a regular GLM with Lasso or Elasticnet Regularisation, mars for multivariate adaptive regression splines (*70*), and cubist for rule-and instance-based regression modelling (*71*)), three random forest algorithms (randomForest (*58*), randomForestSRC (*72*) and ranger (*73*)), and a tree-like node harvesting algorithm (nodeHarvest (*74*)). Extreme gradient boosting and support vector regression were also tested but excluded from the comparison due to its long runtime. Their performance was compared with regards to runtimes, and correlation coefficients between predictions on the test set and the true values. Based on these, we selected the two-layer neural network (brnn) (*55*), multivariate adaptive regression spline (mars) (*70*), and random forest (randomForest) (*58*) algorithms. This ensemble of machine learning models constituted the level 0 learners and was fitted to the initialization set. Out-of-sample predictions from a 10-fold cross validation of each observation were used to fit the level 1 heteroskedastic Gaussian process. As in approach 1, we opted for a Matérn 5/2 kernel and retained the option of changing to a homoscedastic model where necessary. ### Emulator performance We ascertained that both emulators captured the input-output relationship of the simulator by tracking the correlation between true values ***f*** and predicted values ![Graphic][16]on the holdout set of 10% of initial simulations with each iteration (truth vs predicted R2 0.51-0.89 for GP vs 0.37-0.77 for GPSG after initialization, see supplement S1). Transition from exploration to exploitation during adaptive sampling was tracked by recording the distribution of points selected during adaptive sampling in each iteration (Supplementary Fig.s S2 and S3). ### Sensitivity analysis A global sensitivity analysis was conducted on a heteroskedastic GP model with Matérn 5/2 kernel that was trained on all training simulation outputs (n=5,400) from the fitting process. We used the Jansen method of Monte Carlo estimation of Sobol’ sensitivity indices for variance decomposition (*75, 76*) with 20 000 sample points and 1000 bootstrap replicates. Sobol’ indices were calculated for all loss functions ***f*** as well as for their weighted sum ***F*** and in all dimensions. Whilst keeping the number of sample points to as low as possible for computational reasons, we ascertained that first-order indices summed to 1 and total effects >1. We further ensured that the overall results of the Sobol’ analysis were consistent with the results of other global sensitivity analyses, namely the relative parameter importance derived from training a random forest (Supplementary Fig. S32). ### Synthetic data validation Synthetic field data was generated by forward simulation using the final parameter sets from each optimization process. The two optimization algorithms were run anew using the respectively generated synthetic data to calculate the goodness of fit statistics. The parameter sets retrieved by the validation were compared against the parameterization yielded by the optimization process. ### Epidemiological outcome comparison We conducted a small experiment to compare key epidemiological outcomes from the new parameterizations with the original model and that detail in a four malaria model comparison in Penny et al. 2016 (*12*). We simulated malaria in archetypical transmission and seasonality settings using the different parameterizations. The experiments were set up in a full-factorial fashion, considering the simulation options described in Table 1. Monitored outcomes were the incidence of uncomplicated, severe disease, hospitalizations, and indirect and direct malaria mortality over time and by age, prevalence over time and by age, the prevalence-incidence relationship, and the EIR-prevalence relationship. Simulations were conducted for a population of 10,000 individuals over 10 years. View this table: [Table 1:](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/T1) Table 1: Full experimental design in setting archetypes. Experiments were run at 36% probability that an infected individual receives effective care within 14 days. ### Software Consistent with previous calibration work, we used OpenMalaria version 35, an open-source simulator written in C++ and further detailed in full in the supplement, OpenMalaria wiki ([https://github.com/SwissTPH/openmalaria/wiki](https://github.com/SwissTPH/openmalaria/wiki)) or in the original publications (*3, 21, 24*). Calibration was performed using R 3.6.0. For the machine learning processes, all algorithms were accessed through the mlr package version 2.17.0(*66*). The heteroskedastic Gaussian process utilised the hetGP package under version 1.1.2(*65*). The sensitivity analysis was conducted using the soboljansen function of the sensitivity package version 1.21.0 in R (*77*). All algorithms were adapted to the operating system (CentOS 7.5.1804) and computational resources available at the University of Basel Center for Scientific Computing, SciCORE, which uses a Slurm queueing system. The full algorithm code is available on GitHub ([https://github.com/reikth/BayesOpt_Calibration](https://github.com/reikth/BayesOpt_Calibration)) and can be easily adapted to calibrate any simulation model. The number of input parameters and objective functions are flexible. Thus, to adapt the code to other simulators, code should be updated to run the respective model simulator, and tailored to user’s operating system. Further requirements to adapt the workflow are sufficient calibration data, and a per-objective goodness-of-fit metric. ## Data Availability All code is available on git: [https://github.com/reikth/BayesOpt\_Calibration](https://github.com/reikth/BayesOpt_Calibration) ## DATA AND CODE AVAILABILITY Code is publicly available on GitHub under [https://github.com/reikth/BayesOpt\_Calibration](https://github.com/reikth/BayesOpt_Calibration) and all calibration data was detailed in (*24*) and further available from the researchers on request. ## FUNDING The work was funded by the Swiss National Science Foundation through SNSF Professorship of MAP (PP00P3_170702) supporting MAP, MG, and LB. TR was supported by Bill & Melinda Gates Foundation Project OPP1032350. EC’s research is supported by funding from the Bill and Melinda Gates Foundation to Curtin University (Opportunity ID: OPP1197730). ## AUTHOR CONTRIBUTIONS MAP and EC conceived the study. Algorithm development by EC, TR, MAP, and SF with implementation and preparation for sharing on GitHub by TR. Loss functions by MAP and TAS. Sensitivity analysis by TR with inputs from MG and LB. First draft was written by TR and MAP, all authors contributed to writing and interpretation of results and approved the final manuscript. ## COMPETING INTERESTS Authors declare no competing interests. ## SUPPLEMENTARY MATERIALS Supplementary Texts 1 and 2 Supplementary Figures S1-S32 Supplementary Tables S1-S6 ## Supplementary Information for ## 1 SUPPLEMENTARY TEXT 1: MALARIA TRANSMISSION MODEL ### 1.1 Main features We test our calibration algorithm on OpenMalaria, an individual-based model of malaria dynamics. To provide context of the model’s structure and the role of the fitted parameters (see supplementary text 1), we here briefly describe its main features and key equations. This description is adapted from that provided in Smith et al. 2012 (*1*) and Smith et al. 2006 (*2*). Full details of all model components can be found in *The American Journal of Tropical Medicine and Hygiene*, Volume 75, Issue 2 Supplement (2006). OpenMalaria features discrete individual-based stochastic simulations of malaria in humans in 5-day time steps. Every infection and individual are characterized by a set of continuous state variables, namely, parasite densities, infection durations, and immune status. Key processes and relationships regarding the course of infection simulated by model include the attenuation of inoculations, acquired pre-erythrocytic immunity, acquired blood-stage immunity, morbidity (acute and severe) and mortality (malaria-specific and indirect), anemia, and the infection of vectors as a function of parasite densities in the human. Other model components include a vector model and a case management system. All individual components have previously been well documented (*1, 2*). A visual summary of the model with references to further details on each component is provided in Fig. S1. In our current recalibration only the original (base) model variant is used to test our new approach (*1*). Parameters estimated during the calibration process are highlighted and summarized in Table S1 at the end of this section. Other parameter values were drawn from the literature or were calibrated to separate data: for example, the empirical parasite density model of Maire et al. 2006 (*3*) was calibrated to malariatherapy (*4*) data and not recalibrated at the population level. ![Figure S1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F4.medium.gif) [Figure S1.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F4) Figure S1. Visual summary of OpenMalaria with references to original publications on the model components. Adapted from Smith et al. 2006, Fig.3 (*2*).References from top to bottom and left to right (Attenuation of inoculations (*5*)),(Acquired pre-erythrocytic immunity (*5*)),(Infection of vectors [(*6*), (*7*)], Acquired blood-stage immunity (*3*), Anemia (*8*)),(Indirect mortality (neonatal) (*9*), Acute malaria morbidity (*10*), severe malaria morbidity (*11*)),(Indirect mortality excluding neonatal (*11*), Malaria specific mortality (*11*)) ### 1.2 Infection of the human host The seasonal pattern of entomological inoculation rate (EIR) determines seasonal pattern of transmission and thus the parasite densities in the individual, modified by natural or acquired immunity and interventions (*2*). #### 1.2.1 Differential feeding by mosquitoes depending on body surface area In the base model, the expected number of entomological inoculations experienced by individual *i* of age *a* at time *t* is ![Formula][17] where *E**max*(*t*) refers to the annual entomological inoculation rate (EIR) computed from human bait collections on adults and *A*(), is the individual’s availability to mosquitoes, assumed to be proportional to average body surface area, depending only on age *a*(*i, t*). *A*(*a*(*i,t*)) increases with age up to age 20 years where it reaches a value of *A**max* (the average body surface of people ≥ 20 years old in the same population). The biting rate in relation to human weight is based on data from The Gambia published by Port and others (*12*), where the proportion of mosquitoes that had fed on a host were analyzed in relation to the host’s contribution to the total biomass and surface area of people sleeping in one mosquito net (*5*). ### 1.2.2 Control of pre-erythrocytic stages The number of infective bites received per unit time for each individual *i*, adjusted by age, is given by Eq. 1 above. A survival function *S*(*i, t*) defines the probability that the progeny of an inoculation survives to give rise to a patent blood stage infection, i.e. the proportion of inoculations that result in infections or the susceptibility of individual *i* at time *t*. The force of infection is modelled as ![Formula][18] where *E**a*(*i, t*) is the expected number of entomological inoculations endured by individual *i* at time *t*, adjusted for age and individual factors, and the number of infections *h*(*i, t*) acquired by individual ! in five-day time step *t*, follows a Poisson distribution: ![Formula][19] The susceptibility of individual *i* at time *t, S*(*i, t*)is defined as: ![Formula][20] where ![Graphic][21]and *S*∞ are constants representing the lower limit of success probability of inoculations in immune individuals, critical value of cumulative number of entomologic inoculations, critical value of *E**a*(*i, t*), steepness of relationship between success of inoculation and *X**p*(*i, t*), and, the lower limit of success probability of inoculations at high where *E**a*(*i, t*), respectively. Here ![Formula][22] *S*∞ and *E*∗ are fixed to *S*∞ = 0.049, and *E*∗ = 0.032 inoculations/person-night and are detailed in (*5*). #### 1.2.3 Course of infection in the human host The model for each individual infection N in host ! comprises a time series of parasite densities. The base model for infection within humans is described in Maire et al. 2006 (*3*). In brief, the duration of each infection, *τ**max* is sampled from ![Formula][23] parameterised against malaria therapy data (*4*) and detailed in Maire et al. 2006 (*3*). In the absence of previous exposure or concurrent infections, the log density of infection *j* in host *i* at each time point, *τ* = 0,1, …, *τ**max* (*i*, N) is normally distributed with expectation ![Formula][24] Where *y**G*(*τ, τ**max*) is taken from a statistical description of parasite densities in malariatherapy patients and *d*(*i*) describes between-host variation with a log-normal distribution with variance ![Graphic][25]. We consider the possibility of multiple concurrent infections in the same individual at the same time. Exposure to asexual blood stages is measured by ![Formula][26] where *Y*(*i, τ*) is the total parasite density of individual *i* at time *τ* and *y*(*i, j, τ*) is the density of infection N in individual ! at time *τ* and ![Formula][27] In the presence of previous exposure and co-infection, the expected log density for each concurrent infection is then: ![Formula][28] where *M*(*i, t*)is the total multiplicity of infection of in individual *i* at time *t*, and ![Formula][29] where ![Graphic][30] (note that a continuous time approximation to this is given in the original publications (*3, 5*) and hence measures the cumulative parasite load. Furthermore ![Formula][31] where,![Graphic][32], the number of inoculations since birth, excluding the one under consideration, which measures the diversity of inocula experienced by the host up to the time point under consideration. ![Formula][33] which measures the effect of maternal immunity. ![Graphic][34], and *α**m* are all constants estimated in the fitting process. These constants are described in Table S1, or further in Maire et al. 2006 (*3*). Variation within individuals described as ![Graphic][35], where ![Formula][36] and ![Graphic][37] and ![Graphic][38]are constants, described in Table S1. The simulated density of infection *j* in individual *i* at time *τ* is then drawn from a normal distribution: ![Formula][39] The total density of all infections in individual *i* at time t is then the sum of the densities of concurrent infections *j* ![Formula][40] #### 1.2.4 Infectivity of the human host The model infectivity of the human host is described in Ross 2006 where infectivity of individual I at time t is given by the distributed lag model: ![Formula][41] where *t* is in 5-day units and ![Formula][42] where ![Graphic][43] are constants representing contributions of past infections to gametocyte densities (detailed in Table S1), and to be calibrated at the population level. We define ![Formula][44] where Φ is the cumulative normal distribution, ![Graphic][45] is the density of female gametocytes necessary for infection of the mosquito, and ![Graphic][46] is constant (depending on the blood meal volume, gametocyte viability and system variability). Thus, the proportion of mosquitoes infected by individual *i* at time *t* is defined as ![Formula][47] and the probability of a mosquito becoming infected during any feed is ![Formula][48] where *η* is a constant scale factor and to be calibrated. We define ![Graphic][49] as the value of *k**u*(*t*) in the simulation of an equilibrium scenario to which an intervention has been applied. Let ![Graphic][50] be the corresponding entomologic inoculation rate. ![Graphic][51] and ![Graphic][52] are the corresponding values for the intervention scenario. Then ![Formula][53] where *l**v* corresponds to the duration of the sporogenic cycle in the vector, which we approximate with two time steps (10 days). ![Graphic][54] is the total vectorial capacity) ### 1.3 Morbidity In order to simulate the clinical state of individual *i* at time *t*, for each five-day time step 5 independent samples from the simulated parasite density distribution are drawn for each concurrent infection *j*. #### 1.3.1 Acute morbidity (uncomplicated clinical cases) The model for an episode of acute morbidity was originally described in (*10*) and occurs in individual *I* at time *t* with probability ![Formula][55] where *Y*∗is the pyrogenic threshold and *Y**max* is the maximum density of five daily densities sampled during the five-day interval *t*. The pyrogenic threshold changes over time following ![Formula][56] where *f*1(*Y*(*i, t*)) is a function describing the relationship between accrual of tolerance and the parasite density *Y*(*i, y*); *f*2(*Y*∗(*i, t*)) describes the saturation of this accrual process at high values of *Y*∗ and ![Graphic][57] determines the decay threshold with first-order kinetics, ensuring that the parasite tolerance is short-lived (*10*). Here *f*1(*Y* (*i, t*)) is defined to ensure that the stimulus is not directly proportional to *Y* but rather that it asymptotically reaches a maximum at high values of *Y*: ![Formula][58] At high values of *Y*∗, a higher parasite load is required to achieve the same increase: ![Formula][59] Thus, the pyrogenic threshold *Y*∗is defined to follow ![Formula][60] and the initial condition ![Graphic][61] at the birth of the host, where ![Graphic][62] and ![Graphic][63] are targets of the calibration, and are defined in Table S1. #### 1.3.2 Severe disease The model for severe disease was described in Ross et al 2006 (*11*) and two different classes of severe episodes are considered by the model, *B*1and *B*2. *P**B*1(*i, t*) is the probability that an acute episode (*A*) is of class *B*1 and ![Formula][64] where ![Graphic][65] is a constant to be calibrated and *H*(*i, t*) is the clinical status of individual *i* at time *t*. Class *B*2 of severe malaria episodes occurs when an otherwise uncomplicated episode coincides with some other insult, which occurs with risk ![Formula][66] where *F* is the limiting value of *F*(*a*(*i, t*)) at birth and ![Graphic][67] is the age at which it is halved and both are to be calibrated. The probability that individual *i* experiences an episode belonging to class *B*2 at time *t*, conditional on there being a clinical episode at that time is ![Formula][68] The age ant time specific risk of severe malaria morbidity conditional on a clinical episode is then given by ![Formula][69] #### 1.3.3 Mortality Malaria deaths in hospital are a random sample of admitted severe malaria cases, with age-dependent sampling fraction *Q**h*(*a*), the hospital case fatality rate, derived from the data of Reyburn et al (2004) (*13*). The original model was described in Ross et al. 2006 (*11*). The severe malaria case fatality in the community for age group *a, Q**c*(*a*) is estimated as ![Formula][70] where ϕ1 the estimated odds ratio for death in the community compared to death in in-patients is an age-independent constant to be calibrated and *Q**h*(*a*) is the hospital case fatality rate. The total malaria mortality is the sum of the hospital and community malaria deaths. The risk of neonatal mortality attributable to malaria (death in class *D*1) in first pregnancies is set equal to 0.3 *μ**PG* where ![Formula][71] where *x**PG* is related to *x**MG*, the prevalence in simulated individuals 20-24 ears of age via ![Formula][72] and ![Graphic][73] and ![Graphic][74] are constants to be calibrated and are detailed in Table S1. An indirect death in class *D*2 is provoked at time *t*, conditional on there being a clinical episode at that time with probability *P**D*2(*i, t*) where ![Formula][75] and ![Formula][76] where *Q**D* is limiting value of ![Graphic][77] at birth and ![Graphic][78] is a constant to be calibrated. Deaths in class *D*2 occur 30 days (six time steps) after the provoking episode. View this table: [Table S1:](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/T2) Table S1: Names and details of OpenMalaria core parameters. GA-O = Genetic algorithm optimization, GP-BO = Gaussian process-based Bayesian optimization, GPSG-BO = Gaussian process stacked generalization-based Bayesian optimization ## 2 SUPPLEMENTARY TEXT 2: CALIBRATION APPROACH AND DATA SUMMARY A comprehensive epidemiological calibration dataset was collated in order to parameterize OpenMalaria. This calibration dataset covers a total of eleven different epidemiological relationships (or objectives for fitting) that span important aspects of the natural history of malaria. Data were collated from different settings (see Table S2 for summary) and were detailed in the original model descriptions (*2, 10*) and a later parameterization (*1*). A total of 61 simulation scenarios were setup to parameterize OpenMalaria, constructed to simulate the study surveys and study sites that yielded the calibration dataset. The study site observations were replicated in OpenMalaria by reproducing the timing of the surveys and their endpoints (such as prevalence and incidence) and matching simulation options to the setting with regards to transmission intensity and seasonality, vector species, treatment seeking behavior and anti-malarial interventions. The objectives and data are further detailed below. The parameter estimation process is a multi-objective optimization problem with each of the epidemiological quantities in Table S2 representing one objective. The aim of the optimization is to find a parameter set that maximizes the goodness of fit by minimizing a loss statistic computed as the weighted sum of the loss functions for each objective. Building a weighted average reduces the multiple loss terms to a single overall loss statistic, defined as: ![Formula][79] where *f**ij*(***θ***) is the loss function for parameter vector ***θ***, epidemiological quantity *i* and dataset *j*, and the weights *w**i* were chosen so that different epidemiological quantities contribute approximately equally to *F*(***θ***). For the current calibration, we utilised the loss functions from Smith et al. 2012 (*1*), the loss function *f**i*(***θ***) for each objective *i* use either (negative) log-likelihoods or Residual Sum of Squares (RSS) with an unknown minimum. We did not update these loss-functions in order to compare to our previous approaches. The likelihood functions are given by ![Formula][80] where the observed values are *x*1, …, *x**n* and the model parameters ***θ***. In practice, it is easier to work with the log likelihood, namely ![Formula][81] The loss functions *f**i*(***θ***) used for each objective are detailed in the following sections. ### 2.1 Objectives: Epidemiological data and loss functions Below we described each fitting objective in terms of the data (setting, surveys, observations, references) along with the associated loss function and original references. Table S3 provides an overview of the 61 simulation scenarios used for calibration, and which objective they contribute to. #### 2.1.1 Age pattern of incidence after intervention ##### 2.1.1.1 Data The data used for the calibration of objective 1 (Age pattern of incidence) consists of eight cross-sectional surveys of infection rates by age and EIR in Matsari village, capturing 12 age groups each. Matsari village was monitored entomologically for four years (Nov 1970 - Nov1973) during the Garki Project and multiple anti-malaria interventions were administered (*14*). From October 1970 to March 1972 (the baseline / pre-intervention phase), eight cross-sectional malariologic surveys of the whole village population and intensive entomologic surveillance (human bait collection of mosquitoes and dissections of the mosquito salivary glands for sporozoites) were carried out. The latter was used to estimate a baseline transmission intensity of 67 inoculations per person per year (EIR) and to derive seasonal transmission patterns. Mid-1972 marked the beginning of the intervention phase, during which an additional eight surveys were carried out at 10-week intervals (surveys 9-16). During this time, indoor residual spraying with Propoxur was carried out comprehensively in the village, along with mass treatment of the population with Sulfadoxine-pyrimethamine at 10 week-intervals immediately after assessment of individuals’ parasitologic status. The experimental setup is summarised in Fig. 3 of Smith et al 2006 (*5*). Incidence data (number of patent infections and number of hosts by age) from surveys 9-16 was used for our calibration. ###### Sites and scenario numbers Matsari, Nigeria (30) ###### Original reference detailing data and model fits *Smith TA, Maire N, Dietz K, Killeen GF, Vounatsou P et al. Relationship between the entomological inoculation rate and the force of infection for Plasmodium falciparum malaria. Am J Trop Med Hyg. Volume 75, No. 2 Supplement. 2006 (5)* ##### 2.1.1.2 Loss function: Binomial Log Likelihood We denote the Binomial log likelihood for this objective to be ![Formula][82] where *a* is the number of age groups,. the number of surveys, *p**j,k* the scenario data number of parasite positive hosts and *H**j,k* the scenario data number of hosts for age group *k* and survey *j*. Parameter ![Graphic][83] is associated with the model predictions and is given by ![Formula][84] where ![Graphic][85] are the predicted number of parasite positive hosts and ![Graphic][86] the predicted number of hosts for age group *k* and survey *j*. #### 2.1.2 Age patterns of prevalence ##### 2.1.2.1 Data The data used for the calibration of objective 2 (age-patterns of prevalence) consists of six cross-sectional malariology surveys conducted in the Rafin Marke, Matsari, Sugungum villages in Nigeria 1970-1972 (12 age groups each, part of the Garki Project during the pre-intervention period) (*14*), Navrongo in Ghana 2000 (12 age groups) (*15*) and Namawala 1990-1991 (*16*) and Idete in Tanzania (11 and 6 age groups, respectively) 1992-1993 (*17*). In all study sites, annual transmission intensity (EIR) and seasonal patterns were assessed using light trap or human night bait collections and dissections of the salivary glands (see Fig. 2 in Maire et al. 2006 (*3*)). In all sites except Idete, the health system at the time of the surveys treated only a small proportion of the clinical malaria episodes. In the Idete, the village dispensary was assumed to treat approximately 64% of clinical malaria (based on the published literature). During simulation, prevalence was defined by comparing each predicted parasite density with the limit of detection used in the actual study. ###### Sites and scenario numbers Sugungum, Nigeria (24); Rafin-Marke, Nigeria (28); Matsari, Nigeria (29); Idete, Tanzania (31); Navrongo, Ghana (34); Namawala, Tanzania (35) ###### Original reference detailing data and model fits *Maire N, Smith TA, Ross A, Owusu-Agyei S, Dietz K, et al. A model for natural immunity to asexual blood stages of Plasmodium falciparum malaria in endemic areas. Am J Trop Med Hyg. Volume 75, No. 2 Supplement. 2006 (3)* ##### 2.1.2.2 Loss function: Binomial Log Likelihood We denote the binomial log likelihood for each scenario of this objective to be ![Formula][87] where *a* is the number of age groups, *s* the number of surveys, *p**j,k* the scenario data number of parasite positive hosts and *H**j,k* the scenario data number of hosts for age group k and survey j. Parameter *p**j,k* is associated with the model predictions and is given by ![Formula][88] where ![Graphic][89] are the predicted number of parasite positive hosts and ![Graphic][90] the predicted number of hosts for age group *k* and survey *j*. #### 2.1.3 Age patterns of parasite density ##### 2.1.3.1 Data The same data sources as for objective 2 (age pattern of prevalence) were used for calibration of objective 3 (age pattern of parasite density). Parasite densities in sites that were part of the Garki project (Sugungum, Rafin-Make and Matsari, Nigeria) were recorded by scanning a predetermined number of microscope fields on the thick blood film and recording how many had one or more asexual parasites visible. These were converted to numbers of parasites visible by assuming Poisson distribution for the number of parasites per field and a blood volume of 0.5 mm3 per 200 fields. In the other studies (Idete and Namawala, Tanzania and Navrongo, Ghana), parasites were counted against leukocytes and converted to nominal parasites/microliter assuming the usual standard of 8,000 leukocytes/microliter. The biases in density estimates resulting from these different techniques were accounted for by multiplying the observed parasite densities with constant values estimated for Garki (*v*) and non-Garki (*v*1) studies to rescale them to the values in malariatherapy patients (*18*). ###### Sites and scenario numbers Sugungum, Nigeria (pre-intervention, 24); Rafin-Marke, Nigeria (pre-intervention, 28); Matsari, Nigeria (pre-intervention, 29); Idete, Tanzania (31); Navrongo, Ghana (34); Namawala, Tanzania (35) ###### Original reference detailing data and model fits *Maire N, Smith TA, Ross A, Owusu-Agyei S, Dietz K, et al. A model for natural immunity to asexual blood stages of Plasmodium falciparum malaria in endemic areas. Am J Trop Med Hyg. Volume 75, No. 2 Supplement. 2006* (*3*) ##### 2.1.3.2 Loss function: Log-normal log likelihood For objective 3 (age pattern of parasite densities) we denote the log-Normal log likelihood for each scenario to be ![Formula][91] where *n* is the number of observations in the data set, *ρ* = exp(−0.5 log(2 *π*)), a constant from the log-normal likelihood, RSS is the residual sum of squares given by ![Formula][92] and *σ* is the standard deviation given by ![Formula][93] Here, *v* is the appropriate density bias, which is a fitting parameter, *a* is the number of age groups, *s* is the number of surveys, *P**j,k* the scenario number of parasite positive hosts, and *Y**j,k* the sum of the log densities, ![Graphic][94] the predicted number of parasite positive hosts and ![Graphic][95] the predicted sum of the log densities for age group *k* and survey *j*. The density bias are fitting parameters *v* and *v*1. #### 2.1.4 Age pattern of number of concurrent infections ##### 2.1.4.1 Data For objective 4 (age pattern of number of concurrent infections), the dataset from Navrongo, Ghana (also used in the calibration of objectives 2 and 3) was used to calibrate to the total numbers of distinct parasite infections in one individual in each age group, and at each survey. Distinct infections were detected by polymerase chain reaction-restriction fragment length polymorphism in the sampled individuals. ###### Sites and scenario numbers Navrongo, Ghana (34) ###### Original reference detailing data and model fits *Maire N, Smith TA, Ross A, Owusu-Agyei S, Dietz K, et al. A model for natural immunity to asexual blood stages of Plasmodium falciparum malaria in endemic areas. Am J Trop Med Hyg. Volume 75, No. 2 Supplement. 2006 (3)* ##### 2.1.4.2 Loss function: Poisson Log Likelihood Assuming that both the data and the simulations are Poisson distributed about the correct value and thereby also allowing for over-dispersion, we denote the Poisson log likelihood for each scenario to be for the objective of age pattern of number of concurrent infections to be ![Formula][96] where *a* is the number of age groups, *s* the number of surveys, *Pn**j,k* the scenario data total patent infections for age group *k* and survey *j*. Parameter *λ* *j,k* is associated with the model predictions and is given by ![Formula][97] where ![Graphic][98] are the predicted total of patent infections and ![Graphic][99] the predicted number of hosts for age group *k* and survey *j* and *H**j,k* is the scenario data number of hosts for age group *k* and survey *j*. #### 2.1.5 Age pattern of incidence of clinical malaria ##### 2.1.5.1 Data Two distinct datasets representing three study sites (Table S4) were used for the calibration of objective 5 and objective 6 (age pattern of incidence of clinical malaria). For Objective 5, the dataset contains data on the age pattern of clinical episodes in the villages of Ndiop and Dielmo in Senegal (*19, 20*). During the study period of July 1990 - June 1992, the village populations were visited daily to detect and treat any clinical malaria attacks with quinine. Cases were detected by reporting of symptoms (fever) during daily active case detection and subsequent thick blood smear microscopy. Only symptomatic individuals (axillary temperature ≥ 38.0°C or rectal temperature ≥ 38.5°C). Due to the active case detection and rapid treatment all symptomatic episodes are assumed to be effectively treated in these villages during the study period. No effective treatment of clinical malaria was assumed prior to the study period. The annual patterns of transmission were replicated as reported by Charlwood et al (1998) (*21*). A proportion *P**t* =35.75% are assumed to be treated effectively in Idete. As all individuals reporting to the village dispensary were treated presumptively with chloroquine, this proportion corresponds to the proportion of episodes reported to the village dispensary. ###### Sites and scenario numbers Ndiop, Senegal (232), Dielmo, Senegal (233) ###### Original reference detailing data and model fits *Smith TA, Ross A, Maire N, Rogier C, Trape J-F et al. An epidemiologic model of the incidence of acute illness in Plasmodium falciparum malaria. Am J Trop Med Hyg. Volume 75, No. 2 Supplement. 2006 (10)* ##### 2.1.5.2 Loss function: RSS-biased We denote a loss function based on biased residual sum of squares: ![Formula][100] where *a* is the number of age groups, *s* the number of surveys, *s*1 the initial survey number, and *R* is the residual given by ![Formula][101] where *I**j,k* is the observed recorded incidence rate, ![Graphic][102] are the predicted total cases (severe and uncomplicated), ![Graphic][103] the predicted number of hosts for age group *k* and survey *j* and *μ* is a bias related to the scenario. For scenarios 232 and 233 (representing Ndiop and Dielmo, Senegal) this bias is *μ* = 5 indicating the duration in years for which episodes are collected. For scenario 49 in Objective 6 (Idete, Tanzania) the bias is *μ* = 0.357459 and represents the proportion of episodes reported to the village dispensary. #### 2.1.6 Age pattern of incidence of clinical malaria: infants ##### 2.1.6.1 Data Objective 6 (age pattern of incidence of clinical malaria in infants) is informed by a dataset on incidence that contains passive case detection data on the age-incidence in infants recorded at the health centre in Idete, Tanzania from June 1993-October 1994 (*17*). The annual patterns of transmission were replicated as reported by Charlwood et al (1998) (*21*). ###### Sites and scenario numbers Idete, Tanzania (49)) ###### Original reference *Smith TA, Ross A, Maire N, Rogier C, Trape J-F et al. An epidemiologic model of the incidence of acute illness in Plasmodium falciparum malaria. Am J Trop Med Hyg. Volume 75, No. 2 Supplement. 2006 (10)* ##### 2.1.6.2 Loss function: RSS-biased The loss function for Objective 6 is the same as Objective 5. For scenario 49 (Idete, Tanzania) the bias is *μ* = 0.357459 and represents the proportion of episodes reported to the village dispensary. #### 2.1.7 Age pattern of threshold parasite density for clinical attacks ##### 2.1.7.1 Data Objective 7 (Age pattern of threshold parasite density for clinical attacks), uses the dataset from Dielmo, Senegal (see objective 5) for calibration. The pyrogenic threshold in the (OpenMalaria) predictions is output as the sum of the log threshold values across age groups. The pyrogenic threshold per age group is given as the parasite:leucocyte ratio for recorded incidence of disease. To adjust these densities to the same scale as that used in fitting the simulation model to other datasets, the parasite:leukocyte ratios were multiplied by a factor of 1,416 to give a notional density in parasites/microliter of blood. This number was derived as follows: Parasites were counted against leukocytes and converted to nominal parasites/microliter assuming the usual (though biased) standard of 8,000 leukocytes/microliter. The biases in density estimates resulting from these different techniques was accounted for by multiplying the observed parasite densities with constant values estimated for Garki (*v*) and non-Garki (*v*1) studies to rescale them to the values in malariatherapy patients (*18*).The value 1416 comes from ![Formula][104] where the original *v*1 ≈ 0.18. ###### Sites and scenario numbers Dielmo, Senegal (234) ###### Original reference detailing data and model fits *Smith TA, Ross A, Maire N, Rogier C, Trape J-F et al. An epidemiologic model of the incidence of acute illness in Plasmodium falciparum malaria. Am J Trop Med Hyg. Volume 75, No. 2 Supplement. 2006 (10)* ##### 2.1.7.2 Loss function: RSS-biased (log) For the objective 7 (Age pattern of threshold parasite density for clinical attacks) we denote a residual sum of squares loss function given by (13) with ![Formula][105] where *Y*∗ is the observed pyrogenic threshold, ![Graphic][106] are the predicted sum log pyrogenic threshold, ![Graphic][107] the predicted number of hosts for age group *k* and survey *j* and is a bias related to the scenario. Here, this bias is related to the log parasite/leucocyte ratio and thus *μ* = 1/(8000*v*1) where *v*1 is the non-Garki density bias. #### 2.1.8 Hospitalization rate in relation to prevalence in children ##### 2.1.8.1 Data Data on the relative incidence of severe malaria-related morbidity and mortality in children <9 years old across different transmission intensities were originally collated by Marsh and Snow (1999) (*22*) (Table 4). Data measurements per age group were available as the relative risk of severe disease compared to age group 1 and the proportion/prevalence of severe episodes. A total of 26 entries on the relationship between severe malaria hospital admission rates and *P. falciparum* prevalence were used to calibrate objective 8 (Hospitalisation rate in relation to prevalence in children), each represented in a separate simulation scenario, with one observation per scenario. These are summarised in Table S5. To obtain a continuous function relating hospital incidence rates to prevalence, linear interpolation between data points was performed. To convert hospital incidence rates to community severe malaria incidence, the hospital admission rates was divided by the assumed proportion of severe episodes representing to hospital (48%). There was assumed to be no effective treatment of uncomplicated malaria episodes or malaria mortality. ###### Sites and scenario numbers Bo, Sierra Leone (501); Niakhar, Senegal (502), Farafenni, The Gambia (503); Areas I-V, The Gambia (504-508); Gihanga, Burundi (509); Katumba, Burundi (510); Karangasso, Burkina Faso (511); Kilifi North, Kenya (512); Manhica, Mozambique (514); Namawala, Tanzania (515); Navrongo, Ghana (516); Saradidi, Kenya (517); Yombo, Tanzania (518); Ziniare, Burkina Faso (519); Matsari, Nigeria (520); ITC control, Burkina Faso (521); Mlomp, Senegal (522); Ganvie, Benin (523); Kilifi Town, Kenya (524); Chonyi, Kenya (525); Bandafassi, Senegal (526); Kongodjan, Burkina Faso (527) ###### Original reference detailing data and model fits *Ross A, Maire N, Molineaux L and Smith TA. An epidemiologic model of severe morbidity and mortality caused by Plasmodium falciparum. Am J Trop Med Hyg. Volume 75, No. 2 Supplement. 2006 (11)* ##### 2.1.8.2 Loss function: squared deviation The loss function is denoted as the log of residual sum of squares ![Formula][108] where *a**s* is the access to treatment of severe cases (0.48, estimated in base model), ![Graphic][109] is the scenario predicted rate of severe episodes per 1000 person-years for age group *k* = 1 (0-9 years), and parameter ![Graphic][110] is the interpolated observed rate of severe episodes per 1000 person year given by ![Formula][111] where ![Graphic][112] is the predicted prevalence summed over all surveys, *P**u* and *P**l* are the observed prevalences above and below the predicted prevalence ![Graphic][113], respectively and *R**u* and *R**l* are the corresponding severe episode rates to the observed prevalences. The predicted prevalence is given by ![Formula][114] where ![Graphic][115] is the total number of parasite positive predicted and ![Graphic][116] are the total number of hosts (division by 24 to give mean values). The predicted rate of episodes per 1000 person year is given by ![Formula][117] where ![Graphic][118] is the number of severe cases predicted and with division by 2 to convert to from 2 years to 1 year and the division by 24 to give mean number of hosts. #### 2.1.9 Age pattern of hospitalization: severe malaria ##### 2.1.9.1 Data For objective 9 (Age pattern of hospitalisation), a subset of the data collated by Marsh and Snow (1999) (*22*) (see objective 8) is used. Detailed age-specific severe hospital admission rates were available for 5 of the sites (Table S6). The patterns of incidence by age were summarised by age in 1-4 and 5-9 year-old children and compared with 1-11 month old infants by calculating the relative risk. Of the five sites, four were selected for fitting objective 9 based on the predicted prevalence. Baku, The Gambia was excluded as the very low (2%) prevalence here could not be matched. ###### Sites and scenario number(s) Area V, The Gambia (158); Saradidi, Kenya (167); Ganvie, Benin (173); Bandafassi, Senegal (176) ###### Original reference detailing data and model fits *Ross A, Maire N, Molineaux L and Smith TA. An epidemiologic model of severe morbidity and mortality caused by Plasmodium falciparum. Am J Trop Med Hyg. Volume 75, No. 2 Supplement. 2006 (11)* ##### 2.1.9.2 Loss function: Residual sums of squares for relative risk We denote a loss function based on residual sum of squares: ![Formula][119] where *RR**k* is the relative risk of severe episode for age group *k* compared to age group 1 and ![Graphic][120] is the predictive relative risk for age group k compared to age group 1. The predicted relative risk is given by ![Formula][121] where ![Graphic][122] is the number of severe cases predicted for age group *k* and ![Graphic][123] the total number of hosts for age group *k*. #### 2.1.10 Malaria specific mortality in children (< 5 years old) ##### 2.1.10.1 Data For objective 10 (Malaria specific mortality in children (<5 years old)), a subset of the data collated by Marsh and Snow (1999) (*22*) (see objective 8) was used (*24*). Mortality data were derived from verbal autopsy studies in sites with prospective demographic surveillance and were adjusted for the effect of malaria transmission intensity on the sensitivity and specificity of the cause of death determination. The odds ratio for death of a case in the community relative to that in hospital was estimated by fitting to the malaria-specific mortality rates in children less than five years of age assuming the published hospital case fatality rate. Nine sites for which both malaria-specific mortality rates and seasonal transmission patterns were available were included for calibration. There is one observation per study site and simulation scenario, and predicted values are for one survey at the end of 2 years. ###### Sites and scenario number(s) Bo, Sierra Leone (301); Niakhar, Senegal (302); Farafenni, The Gambia (303); Kilifi North, Kenya (312); Navrongo, Ghana (316); Saradidi, Kenya (317); Yombo, Tanzania (318); Bandafassi, Senegal (326); Kongodjan, Burkina Faso (327) ###### Original reference detailing data and model fits *Ross A, Maire N, Molineaux L and Smith TA. An epidemiologic model of severe morbidity and mortality caused by Plasmodium falciparum. Am J Trop Med Hyg. Volume 75, No. 2 Supplement. 2006* (*11*) ##### 2.1.10.2 Loss function: Residual sums of squares For objective 10 on Malaria specific mortality in children, the loss function minimizes the log sum of squares ![Formula][124] where *DMR*1 is the observed direct mortality rate for age group 1 (0-5 years) and ![Graphic][125] is the predicted direct mortality rate for age group 1. The predicted direct mortality rate is given by ![Formula][126] where ![Graphic][127] is the number of direct malaria deaths cases predicted for age group 1 and ![Graphic][128] the total number of predicted hosts for age group 1. The division by 2 is to convert to yearly rate as the survey was conducted at the end of 2 years. #### 2.1.11 Indirect malaria infant mortality rate ##### 2.1.11.1 Data For objective 11 (indirect malaria infant mortality rate), a subset of the data collated by Marsh and Snow (1999) (*22*) (see objective 8) was used. These constitute a library of sites for which entomologic data were collected at least monthly and all-cause infant mortality rates (IMR) were available. There is one observation per scenario: all cause infant mortality rate (returned as a single number over whole intervention period). ###### Sites and scenario number(s) Bo, Sierra Leone (401); Niakhar, Senegal (402); Area V, The Gambia (408); Karangasso, Burkina Faso (411); Manhica, Mozambique (414); Namawala, Tanzania (415); Navrongo, Ghana (416); Saradidi, Kanya (417); Yombo, Tanzania (418); Mlomp, Senegal (422); Bandafassi, Senegal (426) ###### Original reference detailing data and model fits *Ross A, Maire N, Molineaux L and Smith TA. An epidemiologic model of severe morbidity and mortality caused by Plasmodium falciparum. Am J Trop Med Hyg. Volume 75, No. 2 Supplement. 2006 (11)* ##### 2.1.11.2 Loss function: Residual sums of squares The loss function minimises the log sum of squares: ![Formula][129] where *iDMR*1 the observed indirect mortality rate for age group 1 and ![Graphic][130] is the predicted indirect mortality rate for age group 1. ### 2.2 Tables S3-S4 View this table: [Table S2:](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/T3) Table S2: Epidemiological quantities and data sources used for parameterizing models. (a) Some scenarios are used to predict several outcomes, so the total of this column does not equal the total of 61 scenarios involved in fitting the models. (b) The number of data points is the sum over all scenarios and simulated survey periods of the number of age groups into which the data were disaggregated for comparison with the model predictions. (c) In relation to the EIR specified as a seasonal pattern. (d) Model predictions for this objective are compared with linear interpolations between the field data points. *The number of data points is the sum over all scenarios and simulated survey periods of the number of age groups into which the data were disaggregated for comparison with the model predictions. Table adapted from Table S1 in Smith et al 2012 **(*1*)**. View this table: [Table S3:](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/T4) Table S3: Calibration data for objectives 2-4, age patterns of prevalence, parasite densities, and multiplicity of infection View this table: [Table S4:](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/T5) Table S4: summary of study data set for objective 5: Age pattern of incidence of clinical malaria. View this table: [Table S5.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/T6) Table S5. Settings used for calibrating the incidence of severe malaria. (Adapte from Table1 from Ross et al. 2006 (*11*)) View this table: [Table S6:](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/T7) Table S6: Age-specific period prevalence rates* of severe malaria, severe malaria, severe malaria anaemia and acute respiratory-tract infections from five communities in The Gambia and Kenya. (Adapted from Table 2 from Snow et al 1997 (*23*)) ## 3 EMULATOR PERFORMANCE ![Figure S2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F5.medium.gif) [Figure S2.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F5) Figure S2. GP emulator performance. Emulator predictions vs true values on a holdout set compromising 10% of initial samples in iteration 1. w.sum is the weighted sum *F*, of the 11 objectives. Here, predictions are generated as the weighted sum of individual objective predictions. ![Figure S3.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F6.medium.gif) [Figure S3.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F6) Figure S3. GP emulator performance. Emulator predictions vs true values on a holdout set compromising 10% of initial samples in iteration 30 (final iteration). w.sum is the weighted sum *F*, of the 11 objectives. Here, predictions are generated as the weighted sum of individual objective predictions. ![Figure S4.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F7.medium.gif) [Figure S4.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F7) Figure S4. GPSG emulator performance. Emulator predictions vs true values on a holdout set compromising 10% of initial samples in iteration 1. w.sum is the weighted sum *F*, of the 11 objectives. Here, predictions are generated as the weighted sum of individual objective predictions. ![Figure S5.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F8.medium.gif) [Figure S5.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F8) Figure S5. GPSG emulator performance. Emulator predictions vs true values on a holdout set compromising 10% of initial samples in iteration 23 (final iteration). w.sum is the weighted sum *F*, of the 11 objectives. Here, predictions are generated as the weighted sum of individual objective predictions. ## 4 ADAPTIVE SAMPLING: SELECTED POINTS ### 4.1 GP-BO ![Figure S6.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F9.medium.gif) [Figure S6.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F9) Figure S6. GP-BO sampling behavior. Values in each dimension of the points sampled during adaptive sampling of GP-BO algorithm in iterations 1,10, 20, and 30. ### 4.2 GPSG-BO ![Figure S7.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F10.medium.gif) [Figure S7.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F10) Figure S7. GPSG-BO sampling behavior. Values in each dimension of the points sampled during adaptive sampling of GPSG-BO algorithm in iterations 1,10, 20, and 23. ## 5 OPENMALARIA: FINAL SIMULATOR FIT ![Figure S8.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F11.medium.gif) [Figure S8.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F11) Figure S8. Objective 1: Age pattern of prevalence in Matsari, Nigeria during the intervention. Final simulator fit using the parameter sets yielded using GP-BO and GPSG-BO compared to the previous parameterization (derived using optimization with a genetic algorithm, GA-O). ![Figure S9.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F12.medium.gif) [Figure S9.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F12) Figure S9. Objective 2: Age pattern of prevalence. Final simulator fit using the parameter sets yielded using GP-BO and GPSG-BO compared to the previous parameterization (derived using optimization with a genetic algorithm, GA-O). ![Figure S10.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F13.medium.gif) [Figure S10.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F13) Figure S10. Objective 3: Age pattern of parasite densities (geometric mean). Final simulator fit using the parameter sets yielded using GP-BO and GPSG-BO compared to the previous parameterization (derived using optimization with a genetic algorithm, GA-O). ![Figure S11.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F14.medium.gif) [Figure S11.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F14) Figure S11. Objective 4: Age pattern of number of concurrent infections. Final simulator fit using the parameter sets yielded using GP-BO and GPSG-BO compared to the previous parameterization (derived using optimization with a genetic algorithm, GA-O). ![Figure S12.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F15.medium.gif) [Figure S12.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F15) Figure S12. Objective 5: Age pattern of incidence of clinical malaria in Dielmo and Ndiop, Senegal. Final simulator fit using the parameter sets yielded using GP-BO and GPSG-BO compared to the previous parameterization (derived using optimization with a genetic algorithm, GA-O). ![Figure S13.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F16.medium.gif) [Figure S13.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F16) Figure S13. Objective 6: Age pattern of incidence of clinical malaria in Idete, Tanzania. Final simulator fit using the parameter sets yielded using GP-BO and GPSG-BO compared to the previous parameterization (derived using optimization with a genetic algorithm, GA-O). ![Figure S14.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F17.medium.gif) [Figure S14.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F17) Figure S14. Objective 6. Age pattern of threshold parasite density for clinical attacks. Final simulator fit using the parameter sets yielded using GP-BO and GPSG-BO compared to the previous parameterization (derived using optimization with a genetic algorithm, GA-O). ![Figure S15.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F18.medium.gif) [Figure S15.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F18) Figure S15. Objective 7: Hospitalization rate in relation to prevalence in children. Final simulator fit using the parameter sets yielded using GP-BO and GPSG-BO compared to the previous parameterization (derived using optimization with a genetic algorithm, GA-O). ![Figure S16.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F19.medium.gif) [Figure S16.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F19) Figure S16. Objective 8. Age pattern of hospitalization. Final simulator fit using the parameter sets yielded using GP-BO and GPSG-BO compared to the previous parameterization (derived using optimization with a genetic algorithm, GA-O). ![Figure S17.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F20.medium.gif) [Figure S17.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F20) Figure S17. Objective 9: Direct mortality in children <5 years old. Final simulator fit using the parameter sets yielded using GP-BO and GPSG-BO compared to the previous parameterization (derived using optimization with a genetic algorithm, GA-O). ![Figure S18.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F21.medium.gif) [Figure S18.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F21) Figure S18. Objective 10: All-cause infant mortality rate. Final simulator fit using the parameter sets yielded using GP-BO and GPSG-BO compared to the previous parameterization (derived using optimization with a genetic algorithm, GA-O). ## 6 VALIDATION ![Figure S19.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F22.medium.gif) [Figure S19.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F22) Figure S19. Data recovery validation of posterior estimates. Prior distributions of each parameter and parameter value identified by the optimization algorithm. The final parameter set was used to generate synthetic field data by simulating each of the 61 scenarios with the respective core parameter sets. The simulation outputs were reformatted to match the original field data, generating a synthetic field data set. The optimization with both algorithms was repeated using this synthetic field data. The plot shows the best parameter values in each dimension identified at the end of the validation optimization compared to the values identified in the original optimization. The grey area shows the prior distribution. **A. GP-BO validation. B. GPSG-BO validation** ## 7 OPENMALARIA SIMULATED EPIDEMIOLOGY ![Figure S20.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F23.medium.gif) [Figure S20.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F23) Figure S20. Seasonal pattern assumed for subsequent analyses. The monthly transmission intensity is equivalent to the annual transmission intensity (EIR) scaled by these values and forced to sum to the annual EIR. ![Figure S21.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F24.medium.gif) [Figure S21.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F24) Figure S21. Relationship between EIR and *Pf*PR2-10 under three parameterizations. Solid lines show medians and shaded regions show 95% credible intervals. EIR denotes the entomological inoculation rate. ![Figure S22.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F25.medium.gif) [Figure S22.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F25) Figure S22. Yearly incidence of clinical (uncomplicated) malaria as a function of PfPR2-10 displayed by parameterization and age group. Clinical incidence is presented in terms of the yearly number of events per person. We assume a probability of effective treatment within 14 days of uncomplicated malaria of 36% ![Figure S23.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F26.medium.gif) [Figure S23.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F26) Figure S23. Yearly incidence of total severe malaria as a function of PfPR2-10, displayed by parameterization and age group. Incidence is presented in terms of the yearly number of events in a population of 1000 individuals. It is assumed that 48% of severe malaria cases seek official care at a heath care facility (hospital). We assume a probability of effective treatment within 14 days of uncomplicated malaria of 36% ![Figure S24.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F27.medium.gif) [Figure S24.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F27) Figure S24. Yearly number of malaria-related deaths as a function of PfPR2-10, displayed by parameterization and age group. Malaria mortality incidence is presented in terms of the yearly number of deaths in a population of 1000 individuals. For the OpenMalaria model both deaths directly attributed to malaria (dotted curve) and all deaths associated with malaria (including both deaths directly attributable to malaria and those associated with comorbidities) are shown (full line). See Box S1.2 for definitions of deaths attributable to malaria in the models ![Figure S25.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F28.medium.gif) [Figure S25.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F28) Figure S25. Yearly incidence of clinical malaria in a seasonal transmission setting as a function of age, displayed by transmission intensity (PfPR2-10) and parameterization. Clinical incidence is presented in terms of the yearly number of events per person. The PfPR2-10 categories include simulated prevalences of 2.5-3.5%, 9-10%, 28-32%, and 47-53% labeled as 3%, 10%, 30%, and 50%, respectively. ![Figure S26.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F29.medium.gif) [Figure S26.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F29) Figure S26. Yearly incidence of clinical malaria in a perennial transmission setting as a function of age, displayed by transmission intensity (PfPR2-10) and parameterization. Clinical incidence is presented in terms of the yearly number of events per person. The PfPR2-10 categories include simulated prevalences of 2.5-3.5%, 9-10%, 28-32%, and 47-53% labeled as 3%, 10%, 30%, and 50%, respectively ![Figure S27.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F30.medium.gif) [Figure S27.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F30) Figure S27. Yearly incidence of total severe malaria in a seasonal transmission setting as a function of age, displayed by transmission intensity (PfPR2-10) and parameterization. Incidence is presented in terms of the yearly number of events per 1000 person-years. It is assumed that 48% of severe malaria cases seek official care at a heath care facility (hospital). The PfPR2-10 categories include simulated prevalences of 2.5-3.5%, 9-10%, 28-32%, and 47-53% labeled as 3%, 10%, 30%, and 50%, respectively ![Figure S28.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F31.medium.gif) [Figure S28.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F31) Figure S28. Yearly incidence of total severe malaria in a perennial transmission setting as a function of age, displayed by transmission intensity (PfPR2-10) and parameterization. Incidence is presented in terms of the yearly number of events per 1000 person-years. It is assumed that 48% of severe malaria cases seek official care at a heath care facility (hospital). The PfPR2-10 categories include simulated prevalences of 2.5-3.5%, 9-10%, 28-32%, and 47-53% labeled as 3%, 10%, 30%, and 50%, respectively ![Figure S29.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F32.medium.gif) [Figure S29.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F32) Figure S29. Yearly incidence of malaria-related deaths in a seasonal transmission setting as a function of age, displayed by transmission intensity (PfPR2-10) and parameterization. Malaria mortality incidence is presented in terms of the yearly number of deaths in a population of 1000 individuals. The dashed estimates represent direct malaria deaths, and the solid all malaria deaths (including those attributable to co-morbidities). The PfPR2-10 categories include simulated prevalences of 2.5-3.5%, 9-10%, 28-32%, and 47-53% labeled as 3%, 10%, 30%, and 50%, respectively ![Figure S30.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F33.medium.gif) [Figure S30.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F33) Figure S30. Yearly incidence of malaria-related deaths in a perennial transmission setting as a function of age, displayed by transmission intensity (PfPR2-10) and parameterization. Malaria mortality incidence is presented in terms of the yearly number of deaths in a population of 1000 individuals. The dashed estimates represent direct malaria deaths, and the solid all malaria deaths (including those attributable to co-morbidities). The PfPR2-10 categories include simulated prevalences of 2.5-3.5%, 9-10%, 28-32%, and 47-53% labeled as 3%, 10%, 30%, and 50%, respectively ## 8 LOG PRIOR DISTRIBUTIONS ![Figure S31.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F34.medium.gif) [Figure S31.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F34) Figure S31. Log prior distributions and final posterior estimates. Prior distributions of each parameter and final parameter values identified by each optimization algorithm (GP-BO and GPSG-BO) and compared to the current parameterization (derived using a genetic algorithm, GA). ## 9 RANGER IMPORTANCE ![Figure S32.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/07/08/2021.01.27.21250484/F35.medium.gif) [Figure S32.](http://medrxiv.org/content/early/2021/07/08/2021.01.27.21250484/F35) Figure S32. Random forest importance. Estimated parameter importance indices for all parameters and objectives. The indices were calculated using the ranger random forest package in R. ## ACKNOWLEDGMENTS We acknowledge and thank our colleagues in the Swiss TPH Disease Modelling unit. Calculations were performed at sciCORE ([http://scicore.unibas.ch/](http://scicore.unibas.ch/)) scientific computing core facility at University of Basel ## Footnotes * We updated literature review and methodology references. Additionally we changed to title * Received January 27, 2021. * Revision received July 7, 2021. * Accepted July 8, 2021. * © 2021, Posted by Cold Spring Harbor Laboratory This pre-print is available under a Creative Commons License (Attribution-NonCommercial 4.0 International), CC BY-NC 4.0, as described at [http://creativecommons.org/licenses/by-nc/4.0/](http://creativecommons.org/licenses/by-nc/4.0/) ## REFERENCES 1. 1. D. L. DeAngelis, V. Grimm, Individual-based models in ecology after four decades. F1000prime reports 6, (2014). 2. 2. L. Willem, F. Verelst, J. Bilcke, N. Hens, P. Beutels, Lessons from a decade of individual-based models for infectious disease transmission: a systematic review (2006-2015). BMC Infect Dis 17, 612 (2017). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/s12879-017-2699-8&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) 3. 3. T. Smith et al., Towards a comprehensive simulation model of malaria epidemiology and control. Parasitology 135, 1507–1516 (2008). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1017/S0031182008000371&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=18694530&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000261548800004&link_type=ISI) 4. 4. M. F. Gomes et al., Assessing the international spreading risk associated with the 2014 west african ebola outbreak. PLoS Curr 6, (2014). 5. 5. N. Ferguson et al., Report 9: Impact of non-pharmaceutical interventions (NPIs) to reduce COVID19 mortality and healthcare demand. Imperial College London 10, 77482 (2020). 6. 6. T. Cohen et al., Are survey-based estimates of the burden of drug resistant TB too low? Insight from a simulation study. PLoS One 3, e2363 (2008). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pone.0002363&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=18523659&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) 7. 7. M. E. Halloran et al., Modeling targeted layered containment of an influenza pandemic in the United States. Proc Natl Acad Sci U S A 105, 4639–4644 (2008). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NDoicG5hcyI7czo1OiJyZXNpZCI7czoxMToiMTA1LzEyLzQ2MzkiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMS8wNy8wOC8yMDIxLjAxLjI3LjIxMjUwNDg0LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 8. 8. T. A. Perkins et al., An agent-based model of dengue virus transmission shows how uncertainty about breakthrough infections influences vaccination impact projections. PLoS computational biology 15, e1006710 (2019). 9. 9. N. Chitnis, D. Hardy, T. Smith, A periodically-forced mathematical model for the seasonal dynamics of malaria in mosquitoes. Bull Math Biol 74, 1098–1124 (2012). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s11538-011-9710-0&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22218880&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) 10. 10. E. Cameron et al., Defining the relationship between infection prevalence and clinical incidence of Plasmodium falciparum malaria. Nat Commun 6, 8170 (2015). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ncomms9170&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=26348689&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) 11. 11. P. A. Eckhoff, Malaria parasite diversity and transmission intensity affect development of parasitological immunity in a mathematical model. Malar J 11, 419 (2012). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/1475-2875-11-419&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=23241282&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) 12. 12. M. A. Penny et al., Public health impact and cost-effectiveness of the RTS,S/AS01 malaria vaccine: a systematic comparison of predictions from four mathematical models. Lancet 387, 367–375 (2016). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S0140-6736(15)00725-4&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=26549466&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) 13. 13. H. C. Slater, P. G. Walker, T. Bousema, L. C. Okell, A. C. Ghani, The potential impact of adding ivermectin to a mass treatment intervention to reduce malaria transmission: a modelling study. J Infect Dis 210, 1972–1980 (2014). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/infdis/jiu351&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24951826&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) 14. 14. S. Bhatt et al., Improved prediction accuracy for disease risk mapping using Gaussian process stacked generalization. Journal of The Royal Society Interface 14, 20170520 (2017). 15. 15. P. Winskill, P. G. Walker, J. T. Griffin, A. C. Ghani, Modelling the cost-effectiveness of introducing the RTS,S malaria vaccine relative to scaling up other malaria interventions in sub-Saharan Africa. BMJ Glob Health 2, e000090 (2017). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NToiYm1qZ2giO3M6NToicmVzaWQiO3M6MTE6IjIvMS9lMDAwMDkwIjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjEvMDcvMDgvMjAyMS4wMS4yNy4yMTI1MDQ4NC5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 16. 16. T. D. Nguyen et al., Optimum population-level use of artemisinin combination therapies: a modelling study. Lancet Glob Health 3, e758–766 (2015). 17. 17. O. J. Brady et al., Role of mass drug administration in elimination of Plasmodium falciparum malaria: a consensus modelling study. Lancet Glob Health 5, e680–e687 (2017). 18. 18. W. H. Organization, Malaria vaccine: WHO position paper–January 2016. Weekly Epidemiological Record= Relevé épidémiologique hebdomadaire 91, 33–52 (2016). 19. 19. L. Okell et al., in Malaria Policy Advisory Committee meeting. (2015), pp. 16–18. 20. 20. M. Runge et al., Simulating the council-specific impact of anti-malaria interventions: a tool to support malaria strategic planning in Tanzania. PloS one 15, e0228469 (2020). 21. 21. T. Smith et al., Ensemble modeling of the likely public health impact of a pre-erythrocytic malaria vaccine. PLoS Med 9, e1001157 (2012). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pmed.1001157&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22272189&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) 22. 22. R. E. Bellman, Dynamic programming. (Princeton University Press, ed. 6, 1957). 23. 23. A. Craig, in BBC News. (2003), vol. 2020. 24. 24. T. Smith et al., Mathematical modeling of the impact of malaria vaccines on the clinical epidemiology and natural history of Plasmodium falciparum malaria: Overview. The American journal of tropical medicine and hygiene 75, 1–10 (2006). [FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiRlVMTCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czo2OiI3NS8xLzEiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMS8wNy8wOC8yMDIxLjAxLjI3LjIxMjUwNDg0LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 25. 25. D. E. Goldberg, Genetic algorithms in search. Optimization, and MachineLearning, (1989). 26. 26. P. S. Oliveto, T. Paixão, J. Pérez Heredia, D. Sudholt, B. Trubenová, in Proceedings of the Genetic and Evolutionary Computation Conference 2016. (2016), pp. 1163–1170. 27. 27. C. M. Hazelbag, J. Dushoff, E. M. Dominic, Z. E. Mthombothi, W. Delva, Calibration of individual-based models to epidemiological data: A systematic review. PLoS Comput Biol 16, e1007893 (2020). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pcbi.1007893&link_type=DOI) 28. 28. P. A. Eckhoff, A malaria transmission-directed model of mosquito life cycle and ecology. Malaria journal 10, 303 (2011). 29. 29. P. Eckhoff, P. falciparum infection durations and infectiousness are shaped by antigenic variation and innate and adaptive host immunity in a mathematical model. PloS one 7, e44950 (2012). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pone.0044950&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=23028698&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) 30. 30. P. Eckhoff, Mathematical models of within-host and transmission dynamics to determine effects of malaria interventions in a variety of transmission settings. The American journal of tropical medicine and hygiene 88, 817–827 (2013). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czo4OiI4OC81LzgxNyI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIxLzA3LzA4LzIwMjEuMDEuMjcuMjEyNTA0ODQuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 31. 31. J. T. Griffin et al., Reducing Plasmodium falciparum malaria transmission in Africa: a model-based evaluation of intervention strategies. PLoS Med 7, e1000324 (2010). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pmed.1000324&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20711482&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) 32. 32. J. T. Griffin, N. M. Ferguson, A. C. Ghani, Estimates of the changing age-burden of Plasmodium falciparum malaria disease in sub-Saharan Africa. Nature communications 5, 1–10 (2014). 33. 33. I. Fer et al., Linking big models to big data: efficient ecosystem model calibration through Bayesian model emulation. Biogeosciences (Online) 15, (2018). 34. 34. J. Mockus, in Bayesian Approach to Global Optimization. (Springer, 1989), pp. 125–156. 35. 35. M. C. Kennedy, A. O’Hagan, Bayesian calibration of computer models. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63, 425–464 (2001). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1111/1467-9868.00294&link_type=DOI) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000170353900001&link_type=ISI) 36. 36. A. Chong, K. Menberg, Guidelines for the Bayesian calibration of building energy models. Energy and Buildings 174, 527–547 (2018). 37. 37. R. B. Gramacy et al., Calibrating a large computer experiment simulating radiative shock hydrodynamics. Annals of Applied Statistics 9, 1141–1168 (2015). 38. 38. J. Snoek, H. Larochelle, R. P. Adams, Practical bayesian optimization of machine learning algorithms. Advances in neural information processing systems 25, 2951–2959 (2012). 39. 39. J. Snoek et al., in International conference on machine learning. (2015), pp. 2171–2180. 40. 40. D. H. Wolpert, Stacked generalization. Neural networks 5, 241–259 (1992). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S0893-6080(05)80023-1&link_type=DOI) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1992HL21400005&link_type=ISI) 41. 41. I. M. Sobol, Sensitivity analysis for non-linear mathematical models. Mathematical modelling and computational experiment 1, 407–414 (1993). 42. 42. D. Benkeser, C. Ju, S. Lendle, M. van der Laan, Online cross-validation-based ensemble learning. Statistics in medicine 37, 249–260 (2018). 43. 43. L. Breiman, Stacked regressions. Machine learning 24, 49–64 (1996). 44. 44. M. J. Van der Laan, E. C. Polley, A. E. Hubbard, Super learner. Statistical applications in genetics and molecular biology 6, (2007). 45. 45. J. Sill, G. Takács, L. Mackey, D. Lin, Feature-weighted linear stacking. arXiv preprint arXiv:0911.0460, (2009). 46. 46. N. Srinivas, A. Krause, S. M. Kakade, M. Seeger, Gaussian process optimization in the bandit setting: No regret and experimental design. arXiv preprint arXiv:0912.3995, (2009). 47. 47. E. Baker et al., Analyzing Stochastic Computer Models: A Review with Opportunities. arXiv e-prints, arXiv: 2002.01321 (2020). 48. 48. R. Moriconi, M. P. Deisenroth, K. S. Kumar, High-dimensional Bayesian optimization using low-dimensional feature spaces. Machine Learning 109, 1925–1943 (2020). 49. 49. D. Zhou, L. Li, Q. Gu, in International Conference on Machine Learning. (PMLR, 2020), pp. 11492–11502. 50. 50. R. T. Marler, J. S. Arora, The weighted sum method for multi-objective optimization: new insights. Structural and multidisciplinary optimization 41, 853–862 (2010). 51. 51. M. Binois, R. B. Gramacy, M. Ludkovski, Practical heteroscedastic gaussian process modeling for large simulation experiments. Journal of Computational and Graphical Statistics 27, 808–821 (2018). 52. 52. A. Hadji, B. Szábo, Can we trust Bayesian uncertainty quantification from Gaussian process priors with squared exponential covariance kernel? arXiv preprint arXiv:1904.01383, (2019). 53. 53. F. D. Foresee, M. T. Hagan, in Proceedings of International Conference on Neural Networks (ICNN’97). (IEEE, 1997), vol. 3, pp. 1930–1935. 54. 54. D. J. MacKay, Bayesian interpolation. Neural computation 4, 415–447 (1992). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1162/neco.1992.4.3.415&link_type=DOI) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1992HZ47600010&link_type=ISI) 55. 55. P. Rodriguez, D. Gianola, BRNN: Bayesian regularization for feed-forward neural networks. R package version 0.6, (2016). 56. 56. T. Hastie, R. Tibshirani, J. Friedman, The elements of statistical learning : data mining, inference, and prediction. Springer series in statistics (ed. Corrected printing 2002, 2009), vol. 2, pp. 533 S. 57. 57. L. Breiman, Random forests. Machine learning 45, 5–32 (2001). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1023/A:1010933404324&link_type=DOI) 58. 58. A. Liaw, M. Wiener, L. Breiman, A. Cutler. (2015). 59. 59. N. Chitnis et al., Theory of reactive interventions in the elimination and control of malaria. Malaria journal 18, 266 (2019). 60. 60. T. Reiker, N. Chitnis, T. Smith, Modelling reactive case detection strategies for interrupting transmission of Plasmodium falciparum malaria. Malaria journal 18, 259 (2019). 61. 61. M.-L. Cauwet et al., in International Conference on Machine Learning. (PMLR, 2020), pp. 1338–1348. 62. 62. S. Kucherenko, D. Albrecht, A. Saltelli, Exploring multi-dimensional spaces: A comparison of Latin hypercube and quasi Monte Carlo sampling techniques. arXiv preprint arXiv:1505.02350, (2015). 63. 63. P. Auer, Using confidence bounds for exploitation-exploration trade-offs. Journal of Machine Learning Research 3, 397–422 (2002). 64. 64. E. Brochu, V. M. Cora, N. De Freitas, A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning. arXiv preprint arXiv:1012.2599, (2010). 65. 65. M. Binois, R. Gramacy, hetGP: Heteroskedastic Gaussian process modeling and sequential design in R. (2019). 66. 66. B. Bischl et al., mlr: Machine Learning in R. The Journal of Machine Learning Research 17, 5938–5942 (2016). 67. 67. B. Ripley, W. Venables, M. B. Ripley, Package ‘nnet’. R package version 7, 3–12 (2016). 68. 68. T. Hastie, J. Qian, Glmnet vignette. Retrieved June 9, 1–30 (2014). 69. 69. B. Hofner, A. Mayr, N. Robinzonov, M. Schmid, Model-based boosting in R: a hands-on tutorial using the R package mboost. Computational statistics 29, 3–35 (2014). 70. 70. T. Hastie, R. Tibshirani, F. Leisch, K. Hornik, B. Ripley, mda: Mixture and flexible discriminant analysis. R package version 0.4-4, URL [http://cran.r-project.org/package=](http://cran.r-project.org/package=) mda, (2013). 71. 71. M. Kuhn, R. Quinlan, Cubist: Rule-And Instance-Based Regression. Modeling. R package version (2) 2, (2018). 72. 72. H. Ishwaran, U. B. Kogalur, M. U. B. Kogalur, Package ‘randomForestSRC’. (2020). 73. 73. M. N. Wright, A. Ziegler, ranger: A fast implementation of random forests for high dimensional data in C++ and R. arXiv preprint arXiv:1508.04409, (2015). 74. 74. N. Meinshausen, M. N. Meinshausen, Package ‘nodeHarvest’. (2015). 75. 75. M. J. Jansen, Analysis of variance designs for model output. Computer Physics Communications 117, 35–43 (1999). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S0010-4655(98)00154-4&link_type=DOI) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000079344200005&link_type=ISI) 76. 76. A. Saltelli et al., Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Computer physics communications 181, 259–270 (2010). 77. 77. B. Iooss et al., Package ‘sensitivity’. (2021). 1. 1. T. Smith et al., Ensemble modeling of the likely public health impact of a pre-erythrocytic malaria vaccine. PLoS Med 9, e1001157 (2012). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pmed.1001157&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22272189&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) 2. 2. T. Smith et al., Mathematical modeling of the impact of malaria vaccines on the clinical epidemiology and natural history of Plasmodium falciparum malaria: Overview. The American journal of tropical medicine and hygiene 75, 1–10 (2006). [FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiRlVMTCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czo2OiI3NS8xLzEiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMS8wNy8wOC8yMDIxLjAxLjI3LjIxMjUwNDg0LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 3. 3. N. Maire et al., A model for natural immunity to asexual blood stages of Plasmodium falciparum malaria in endemic areas. The American journal of tropical medicine and hygiene 75, 19–31 (2006). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czoxMzoiNzUvMl9zdXBwbC8xOSI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIxLzA3LzA4LzIwMjEuMDEuMjcuMjEyNTA0ODQuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 4. 4. W. E. Collins, G. M. Jeffery, A retrospective examination of sporozoite-and trophozoite-induced infections with Plasmodium falciparum: development of parasitologic and clinical immunity during primary infection. The American journal of tropical medicine and hygiene 61, 4–19 (1999). [Abstract](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czoxMjoiNjEvMV9zdXBwbC80IjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjEvMDcvMDgvMjAyMS4wMS4yNy4yMTI1MDQ4NC5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 5. 5. T. Smith et al., Relationship between the entomologic inoculation rate and the force of infection for Plasmodium falciparum malaria. The American journal of tropical medicine and hygiene 75, 11–18 (2006). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czoxMzoiNzUvMl9zdXBwbC8xMSI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIxLzA3LzA4LzIwMjEuMDEuMjcuMjEyNTA0ODQuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 6. 6. G. F. Killeen, A. Ross, T. Smith, Infectiousness of malaria-endemic human populations to vectors. The American journal of tropical medicine and hygiene 75, 38–45 (2006). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czoxMzoiNzUvMl9zdXBwbC8zOCI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIxLzA3LzA4LzIwMjEuMDEuMjcuMjEyNTA0ODQuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 7. 7. A. Ross, G. Killeen, T. Smith, Relationships between host infectivity to mosquitoes and asexual parasite density in Plasmodium falciparum. The American journal of tropical medicine and hygiene 75, 32–37 (2006). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czoxMzoiNzUvMl9zdXBwbC8zMiI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIxLzA3LzA4LzIwMjEuMDEuMjcuMjEyNTA0ODQuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 8. 8. I. A. Carneiro et al., Modeling the relationship between the population prevalence of Plasmodium falciparum malaria and anemia. The American journal of tropical medicine and hygiene 75, 82–89 (2006). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czoxMzoiNzUvMl9zdXBwbC84MiI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIxLzA3LzA4LzIwMjEuMDEuMjcuMjEyNTA0ODQuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 9. 9. A. Ross, T. Smith, The effect of malaria transmission intensity on neonatal mortality in endemic areas. The American journal of tropical medicine and hygiene 75, 74–81 (2006). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czoxMzoiNzUvMl9zdXBwbC83NCI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIxLzA3LzA4LzIwMjEuMDEuMjcuMjEyNTA0ODQuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 10. 10. T. Smith et al., An epidemiologic model of the incidence of acute illness in Plasmodium falciparum malaria. The American journal of tropical medicine and hygiene 75, 56–62 (2006). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czoxMzoiNzUvMl9zdXBwbC81NiI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIxLzA3LzA4LzIwMjEuMDEuMjcuMjEyNTA0ODQuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 11. 11. A. Ross, N. Maire, L. Molineaux, T. Smith, An epidemiologic model of severe morbidity and mortality caused by Plasmodium falciparum. The American journal of tropical medicine and hygiene 75, 63–73 (2006). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czoxMzoiNzUvMl9zdXBwbC82MyI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIxLzA3LzA4LzIwMjEuMDEuMjcuMjEyNTA0ODQuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 12. 12. G. Port, P. Boreham, J. H. Bryan, The relationship of host size to feeding by mosquitoes of the Anopheles gambiae Giles complex (Diptera: Culicidae). Bulletin of Entomological Research 70, 133–144 (1980). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1017/S0007485300009834&link_type=DOI) 13. 13.R. H. et al., “The epidemiology of severe malaria due to Plasmodium falciparum at different transmisison intensities in NE Tanzania,” LSTMH Malaria Centre R 2002-2003 (2004). 14. 14. L. Molineaux, G. Gramiccia, W. H. Organization, The Garki project: research on the epidemiology and control of malaria in the Sudan savanna of West Africa. (World Health Organization, 1980). 15. 15. S. Owusu-Agyei, T. Smith, H. P. Beck, L. Amenga-Etego, I. Felger, Molecular epidemiology of Plasmodium falciparum infections among asymptomatic inhabitants of a holoendemic malarious area in northern Ghana. Tropical Medicine & International Health 7, 421–428 (2002). 16. 16. T. Smith et al., Absence of seasonal variation in malaria parasitaemia in an area of intense seasonal transmission. Acta tropica 54, 55–72 (1993). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/0001-706X(93)90068-M&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=8103627&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1993LN14200006&link_type=ISI) 17. 17. A. Y. Kitua et al., Plasmodium falciparum malaria in the first year of life in an area of intense and perennial transmission. Tropical Medicine & International Health 1, 475–484 (1996). 18. 18. W. C. Earle, M. Perez, Enumeration of parasites in the blood of malarial patients. Journal of laboratory and clinical medicine 17, (1932). 19. 19. J.-F. Trape, C. Rogier, Combating malaria morbidity and mortality by reducing transmission. Parasitology today 12, 236–240 (1996). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/0169-4758(96)10015-6&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=15275204&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1996UM60800007&link_type=ISI) 20. 20. J.-F. Trape et al., The Dielmo project: a longitudinal study of natural malaria infection and the mechanisms of protective immunity in a community living in a holoendemic area of Senegal. The American journal of tropical medicine and hygiene 51, 123–137 (1994). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czo4OiI1MS8yLzEyMyI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIxLzA3LzA4LzIwMjEuMDEuMjcuMjEyNTA0ODQuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 21. 21. J. Charlwood et al., Incidence of Plasmodium falciparum infection in infants in relation to exposure to sporozoite-infected anophelines. The American journal of tropical medicine and hygiene 59, 243–251 (1998). [Abstract](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czo4OiI1OS8yLzI0MyI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIxLzA3LzA4LzIwMjEuMDEuMjcuMjEyNTA0ODQuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 22. 22. K. Marsh, R. Snow, Malaria transmission and morbidity. Parassitologia 41, 241 (1999). [PubMed](http://medrxiv.org/lookup/external-ref?access_num=10697862&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) 23. 23. R. W. Snow et al., Relation between severe malaria morbidity in children and level of Plasmodium falciparum transmission in Africa. The lancet 349, 1650–1654 (1997). 24. 24. E. L. Korenromp, B. G. Williams, E. Gouws, C. Dye, R. W. Snow, Measurement of trends in childhood malaria mortality in Africa: an assessment of progress toward targets based on verbal autopsy. The Lancet infectious diseases 3, 349–358 (2003). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S1473-3099(03)00657-1&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=12781507&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000183244300023&link_type=ISI) 25. 25. C. Rogier, D. Commenges, J.-F. Trape, Evidence for an age-dependent pyrogenic threshold of Plasmodium falciparum parasitemia in highly endemic populations. The American journal of tropical medicine and hygiene 54, 613–619 (1996). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czo4OiI1NC82LzYxMyI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIxLzA3LzA4LzIwMjEuMDEuMjcuMjEyNTA0ODQuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 26. 26. P. Vounatsou, T. Smith, A. Kitua, P. Alonso, M. Tanner, Apparent tolerance of Plasmodium falciparum in infants in a highly endemic area. Parasitology 120, 1–9 (2000). [PubMed](http://medrxiv.org/lookup/external-ref?access_num=10726260&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000085922200002&link_type=ISI) 27. 27. J. C. Beier, G. F. Killeen, J. I. Githure, Entomologic inoculation rates and Plasmodium falciparum malaria prevalence in Africa. The American journal of tropical medicine and hygiene 61, 109–113 (1999). [Abstract](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoidHJvcG1lZCI7czo1OiJyZXNpZCI7czo4OiI2MS8xLzEwOSI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIxLzA3LzA4LzIwMjEuMDEuMjcuMjEyNTA0ODQuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 28. 28. G. Barnish et al., Malaria in a rural area of Sierra Leone. I. Initial results. Annals of Tropical Medicine & Parasitology 87, 125–136 (1993). 29. 29. I. D. R. Centre, I. Network, Population and Health in Developing Countries: Population, health and survival at INDEPTH sites. (IDRC, 2002), vol. 1. 30. 30. H. C. Spencer et al., Impact on mortality and fertility of a community-based malaria control programme in Saradidi, Kenya. Annals of Tropical Medicine & Parasitology 81, 36–45 (1987). 31. 31. U. D’Alessandro et al., Mortality and morbidity from malaria in Gambian children after introduction of an impregnated bednet programme. The Lancet 345, 479–483 (1995). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S0140-6736(95)90582-0&link_type=DOI) 32. 32. P. Duboz, J. Vaugelade, M. Debouverie, “Mortalité dans l’enfance dans la région de Niangoloko,” (ORSTOM, Ouagadougou, Burkina Faso, 1989). 33. 33. J. A. Schellenberg et al., KINET: a social marketing programme of treated nets and net treatment for malaria control in Tanzania, with evaluation of child health and long-term survival. Transactions of the Royal Society of Tropical Medicine and Hygiene 93, 225–231 (1999). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S0035-9203(99)90001-9&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=10492745&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) 34. 34. Z. Premji et al., Community based studies on childhood mortality in a malaria holoendemic area on the Tanzanian coast. Acta tropica 63, 101–109 (1997). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S0001-706X(96)00605-5&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=9088423&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F07%2F08%2F2021.01.27.21250484.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1997WN68700003&link_type=ISI) 35. 35. J.-F. Trape et al., Impact of chloroquine resistance on malaria mortality. Comptes Rendus de l’Académie des Sciences-Series III-Sciences de la Vie 321, 689–697 (1998). [1]: F1/embed/inline-graphic-1.gif [2]: F1/embed/inline-graphic-2.gif [3]: F1/embed/inline-graphic-3.gif [4]: F1/embed/inline-graphic-4.gif [5]: F1/embed/inline-graphic-5.gif [6]: F1/embed/inline-graphic-6.gif [7]: /embed/inline-graphic-7.gif [8]: /embed/graphic-4.gif [9]: /embed/inline-graphic-8.gif [10]: /embed/inline-graphic-9.gif [11]: /embed/inline-graphic-10.gif [12]: /embed/inline-graphic-11.gif [13]: /embed/inline-graphic-12.gif [14]: /embed/inline-graphic-13.gif [15]: /embed/inline-graphic-14.gif [16]: /embed/inline-graphic-15.gif [17]: /embed/graphic-7.gif [18]: /embed/graphic-8.gif [19]: /embed/graphic-9.gif [20]: /embed/graphic-10.gif [21]: /embed/inline-graphic-16.gif [22]: /embed/graphic-11.gif [23]: /embed/graphic-12.gif [24]: /embed/graphic-13.gif [25]: /embed/inline-graphic-17.gif [26]: /embed/graphic-14.gif [27]: /embed/graphic-15.gif [28]: /embed/graphic-16.gif [29]: /embed/graphic-17.gif [30]: /embed/inline-graphic-18.gif [31]: /embed/graphic-18.gif [32]: /embed/inline-graphic-19.gif [33]: /embed/graphic-19.gif [34]: /embed/inline-graphic-20.gif [35]: /embed/inline-graphic-21.gif [36]: /embed/graphic-20.gif [37]: /embed/inline-graphic-22.gif [38]: /embed/inline-graphic-23.gif [39]: /embed/graphic-21.gif [40]: /embed/graphic-22.gif [41]: /embed/graphic-23.gif [42]: /embed/graphic-24.gif [43]: /embed/inline-graphic-24.gif [44]: /embed/graphic-25.gif [45]: /embed/inline-graphic-25.gif [46]: /embed/inline-graphic-26.gif [47]: /embed/graphic-26.gif [48]: /embed/graphic-27.gif [49]: /embed/inline-graphic-27.gif [50]: /embed/inline-graphic-28.gif [51]: /embed/inline-graphic-29.gif [52]: /embed/inline-graphic-30.gif [53]: /embed/graphic-28.gif [54]: /embed/inline-graphic-31.gif [55]: /embed/graphic-29.gif [56]: /embed/graphic-30.gif [57]: /embed/inline-graphic-32.gif [58]: /embed/graphic-31.gif [59]: /embed/graphic-32.gif [60]: /embed/graphic-33.gif [61]: /embed/inline-graphic-33.gif [62]: /embed/inline-graphic-34.gif [63]: /embed/inline-graphic-35.gif [64]: /embed/graphic-34.gif [65]: /embed/inline-graphic-36.gif [66]: /embed/graphic-35.gif [67]: /embed/inline-graphic-37.gif [68]: /embed/graphic-36.gif [69]: /embed/graphic-37.gif [70]: /embed/graphic-38.gif [71]: /embed/graphic-39.gif [72]: /embed/graphic-40.gif [73]: /embed/inline-graphic-38.gif [74]: /embed/inline-graphic-39.gif [75]: /embed/graphic-41.gif [76]: /embed/graphic-42.gif [77]: /embed/inline-graphic-40.gif [78]: /embed/inline-graphic-41.gif [79]: /embed/graphic-45.gif [80]: /embed/graphic-46.gif [81]: /embed/graphic-47.gif [82]: /embed/graphic-48.gif [83]: /embed/inline-graphic-42.gif [84]: /embed/graphic-49.gif [85]: /embed/inline-graphic-43.gif [86]: /embed/inline-graphic-44.gif [87]: /embed/graphic-50.gif [88]: /embed/graphic-51.gif [89]: /embed/inline-graphic-45.gif [90]: /embed/inline-graphic-46.gif [91]: /embed/graphic-52.gif [92]: /embed/graphic-53.gif [93]: /embed/graphic-54.gif [94]: /embed/inline-graphic-47.gif [95]: /embed/inline-graphic-48.gif [96]: /embed/graphic-55.gif [97]: /embed/graphic-56.gif [98]: /embed/inline-graphic-49.gif [99]: /embed/inline-graphic-50.gif [100]: /embed/graphic-57.gif [101]: /embed/graphic-58.gif [102]: /embed/inline-graphic-51.gif [103]: /embed/inline-graphic-52.gif [104]: /embed/graphic-59.gif [105]: /embed/graphic-60.gif [106]: /embed/inline-graphic-53.gif [107]: /embed/inline-graphic-54.gif [108]: /embed/graphic-61.gif [109]: /embed/inline-graphic-55.gif [110]: /embed/inline-graphic-56.gif [111]: /embed/graphic-62.gif [112]: /embed/inline-graphic-57.gif [113]: /embed/inline-graphic-58.gif [114]: /embed/graphic-63.gif [115]: /embed/inline-graphic-59.gif [116]: /embed/inline-graphic-60.gif [117]: /embed/graphic-64.gif [118]: /embed/inline-graphic-61.gif [119]: /embed/graphic-65.gif [120]: /embed/inline-graphic-62.gif [121]: /embed/graphic-66.gif [122]: /embed/inline-graphic-63.gif [123]: /embed/inline-graphic-64.gif [124]: /embed/graphic-67.gif [125]: /embed/inline-graphic-65.gif [126]: /embed/graphic-68.gif [127]: /embed/inline-graphic-66.gif [128]: /embed/inline-graphic-67.gif [129]: /embed/graphic-69.gif [130]: /embed/inline-graphic-68.gif