Agent-based modelling of complex factors impacting malaria prevalence

Open AccessThis article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Associated Data

Additional file 1. From in situ to continous model. This file contains more detailed explanation of the data and likelihood, model calibrations, Regression with simulated ABM outputs, the extension to continuous times and a detailed description of the ABM using the ODD protocol.

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Abstract

Background

Increasingly complex models have been developed to characterize the transmission dynamics of malaria. The multiplicity of malaria transmission factors calls for a realistic modelling approach that incorporates various complex factors such as the effect of control measures, behavioural impacts of the parasites to the vector, or socio-economic variables. Indeed, the crucial impact of household size in eliminating malaria has been emphasized in previous studies. However, increasing complexity also increases the difficulty of calibrating model parameters. Moreover, despite the availability of much field data, a common pitfall in malaria transmission modelling is to obtain data that could be directly used for model calibration.

Methods

In this work, an approach that provides a way to combine in situ field data with the parameters of malaria transmission models is presented. This is achieved by agent-based stochastic simulations, initially calibrated with hut-level experimental data. The simulation results provide synthetic data for regression analysis that enable the calibration of key parameters of classical models, such as biting rates and vector mortality. In lieu of developing complex dynamical models, the approach is demonstrated using most classical malaria models, but with the model parameters calibrated to account for such complex factors. The performance of the approach is tested against a wide range of field data for Entomological Inoculation Rate (EIR) values.

Results

The overall transmission characteristics can be estimated by including various features that impact EIR and malaria incidence, for instance by reducing the mosquito–human contact rates and increasing the mortality through control measures or socio-economic factors.

Conclusion

Complex phenomena such as the impact of the coverage of the population with long-lasting insecticidal nets (LLINs), changes in behaviour of the infected vector and the impact of socio-economic factors can be included in continuous level modelling. Though the present work should be interpreted as a proof of concept, based on one set of field data only, certain interesting conclusions can already be drawn. While the present work focuses on malaria, the computational approach is generic, and can be applied to other cases where suitable in situ data is available.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12936-021-03721-2.

Keywords: Computational biology, Socio-economic factors, Agent-based modelling, Prevention of reintroduction, Long-lasting insecticidal nets, Multiscale modelling

Background

Malaria is often regarded as a socio-economic disease associated with poverty and underdevelopment. The incidence of the disease tends to decline with economic development and associated improvement in domestic conditions, such as quality of housing and availability of medical aid [1, 2]. The elimination of malaria in, for instance, Finland was preconditioned on a drop in household size [2–4]. However, malaria imported by visitors and migrants carries the risk of re-introducing malaria in areas that have suitable vectors and climatic conditions [5]. For instance, the proportion of imported malaria cases due to migrants in Europe has recently increased from 14 to 83% [6–9]. It is therefore topical to reconsider various factors controlling the spread of malaria.

Classical compartmental models contain a limited account of the complex processes of malaria transmission dynamics, and more detailed models tend to get over-loaded with model parameters that are difficult to calibrate against real data [10]. Here, an approach to alleviate this dilemma is demonstrated by a combination of individual or agent-based modelling (ABM) strategy together with compartmental modelling. The ABM approach has become popular due to its enhanced realism, flexibility, explicitness and the advantages of spatial simulations with high resolution (see [11]). An agent-based modelling approach is employed in order to simulate the impact of factors such as intervention measures, household size, and the behavioural changes of the vector. The ABM results are then linked to basic dynamic transmission models in order to enable predictions on the level of public health [12–14].

The ABM modelling is done first for a single host in the hut, and then on a household level, with multiple individuals sleeping under the same roof. Subsequently, the household-level model is extended to community-level scenarios, enabling simulations of heterogeneity of mosquito-to-human contact rates due to partial coverage with nets or different household sizes. The crucial impact of socio-economic factors such as household size has been emphasized in [2–4]. The ABM simulations provide a ‘computational laboratory’ where data reflecting the impact of various complex factors can be produced. Upon repeated simulations, the ABM outputs can be used as synthetic data to produce regression models for the factors considered. Here, the focus is on household size, LLIN coverage, and alterations in mosquito behaviour induced by malaria parasite.

The agent-based model simulations are conducted over a ‘snapshot’ time period of one night. The results can be extended to continuous time by inserting the values fitted by the response surfaces as the key coefficients of classical compartmental models. Consequently, the impact of intervention measures or socioeconomic factors can be simulated over longer time periods, and to steady state. This allows for the estimation of the EIR [15] values in a wide variety of transmission scenarios.

The work-flow followed in the present study is summarized in the schematic illustration given in Fig. ​ Fig.1. 1 . The modelling process is iterative as there is back and forth movement from MCMC parameter identification to ABM of mosquito host-seeking behavior, such that the model fits the data well.

Schematic representation of transition from the ABM of mosquito host-seeking behaviour to the continuous modelling. The procedure is conducted separately for each of the mosquito species and the chemicals under study

Other studies have also estimated key parameter statistics from data on experimental hut trials and subsequently employed them in dynamic transmission models to enable public predictions. In the work by Churcher et al. [12], key parameters of the continuous model [16] were estimated using statistical models (such as binomial and mixed effect models) and calibrated with hut-level experiment data. Sherrard-Smith et al. [17] systematically assessed experimental hut data to characterize different indoor residual spray (IRS) product efficacies in terms of mosquito mortality, blood-feeding inhibition and deterrence against Anopheles mosquitoes, when fitted with statistical models. The impacts of IRS assessed from experimental hut trials are extrapolated for public health predictions in areas with different levels of coverage and pyrethroid resistance using the mathematical model of malaria transmission from [16]. Using a slightly different approach, Okumu et al. [13] directly inputted values of relevant parameters from experimental hut trials into their transmission model to make public predictions. This model additionally considered animal hosts (cattle) and predicted community-level impact on malaria transmission at high coverage (80%) using direct data from hut-level trials for various combinations of untreated nets or LLINs with IRS. Another related study [18] considered how coverage with ITNs (from 0 to 100%) influence the intensity of malaria transmission using an elaborated description of the classic feeding cycle model. The approach in this study differs from the above papers as it presents the model of mosquito host-seeking behaviour in a hut in terms of the mosquitoes’ attraction to the host, host-seeking orientation, biting and death rate. The agent-based simulation of mosquito behavior at hut-level in the presence of different insecticides is then calibrated with field data from [19]. The ABM approach enables modelling of behavioral changes typical for infected mosquitoes at the household level and subsequent extension to community-level simulations using households of different sizes. Thus, upon simulating the ABM, the key ODE model parameters are created, unlike in [18] which is model-based and parameter values are mainly assumed. Additionally, the LLIN coverage and household size are elaborately considered to range from 0 to 100%, and 2 to 10 respectively. This approach enables integration of socio-economic factors and the study of malaria prevalence in a population at varied protection levels, while in [12, 13] a certain level of coverage is assumed.

The rest of the paper is organized as follows. In "Methods" section, the basic agent-based modelling approach at the hut level with a single host and its extension are presented. Next, the ABM extension to household and and subsequently to community level are described. Then, the next subsection discusses the regression applied to the outputs of community-level simulations. The extension of the response surface results to continuous time is given in "Results" section. Finally, the discussion is presented in the last section.

Methods

Basic ABM host-seeking model with a single host

A previous work [20] presented an ABM simulation approach for mosquito host-seeking behaviour on hut-level in the presence of LLIN, calibrated for one case of the treatment data from [19]. Here, the model is extended in several ways to make it capable of reproducing the data of other insecticidal treatments, and to enable the extension to continuous time modelling done in "Results" section. The basic modelling approaches utilized in [20] is briefly recalled and the modifications made in the present work is pinpointed. The model developed in this study consists of four basic components, where each of the components features a number of associated attributes, see the summary in Table ​ Table1. 1 . Additionally, the properties that are assigned individually for each of the mosquito agents and updated within the simulation (see Table ​ Table2) 2 ) are listed. The model components and the property list of mosquitoes are described in detail in this subsection.

Table 1

Table 2

Property list of each agent and the relevant model component

Motion and host-seeking

The mosquito attraction model is based on the assumption that a mosquito estimates the direction of odour increase (the gradient) by the mechanism of klinotaxis [21]. During this plume-tracking behaviour, the mosquito samples the host odour at one location, changes location and then repeats the sampling, and uses its memory of the concentrations previously encountered to choose the next position [22, 23]. Imitating this process, the flight of mosquitoes is modelled as a discrete-time correlated random walk. Suppose that a mosquito agent is at position x n - 1 at time step n - 1 . A new position x n is selected by:

where the increment δ W added to x n - 1 is sampled in random direction, with a step size given by a normal distribution N ( x 0 , σ 2 I ) . In the experimental runs, the parameters x 0 , σ were matched to imitate the real flight speed of a mosquito [20]. Mosquitoes are able to identify the host by making use of the olfactory cues that are given off by the host. As a primary stimuli, they move in response to the carbon dioxide ( CO 2 ) exhaled by vertebrates. Additionally, at a short distance to the host, mosquitoes are able to discern by vision, olfaction and by using the heat sensors located around their mouthparts. In general, mosquitoes are unable to detect human prey from a distance greater than 80 m [22]. The concentration of attractive odour emitted from an individual host is modelled as a Gaussian kernel centered at a spatial position of the host x h :

where x denotes the position of the mosquito, and C stands for the concentration that enables a mosquito to sense the host at a distance d ( x , x h ) . Note that the impact of wind is omitted for simplicity. The standard deviation of the Gaussian σ a determines a maximal distance at which the mosquito is able to sense the host.

The mosquito flight is given by the above random walk in the absence of attraction effects towards the host. However, when there are attraction effects, the main features of the Metropolis algorithm is employed in order to simulate the random walk directionally biased by attraction [24]. Suppose that a mosquito takes a step from point x n - 1 to a candidate point x n with respective function values as p n - 1 and p n . Then a new point is accepted with probability:

where p ( x n ) / p ( x n - 1 ) is the ratio of the attraction potential function p ( x ) defined at each point x , which depends on the concentration and other attraction factors. In order to parsimoniously account for other short-distance attraction factors, the attraction potential function is defined as:

with a scaling factor σ acc that depends on the distance to the host. Outside the plume p ( x ) = 1 , so by Eq. 3 all steps are accepted, while closer to the host steps away from the host are increasingly rejected due to activation of the heat sensors. At a short distance to the host this is modelled by a linear scaling factor as:

σ acc ( x , x h ) = σ acc 1 + σ acc 2 d ( x , x h ) , d ( x , x h ) ≤ 80 σ acc max , d ( x , x h ) > 80 .

The above function increases from the minimum value of σ acc 1 with a slope given by the parameter σ acc 2 until it is replaced by a constant which suitably provides a purely random movement outside the concentration plume [20].

Death, poisoning and repellency

The LLINs are assumed to be equipped with repellent and poisoning effects. In the absence of chemical treatment, the total probability of death reduces to the natural mortality rate. The continuous-time mortality rate μ can be transformed into a probability of death per unit time Δ t by:

where Δ t = 2 s is used for all simulations, and a value for μ taken from the literature (see [20] for more details). This conforms with the 34-h natural mortality rates reported for Anopheles gambiae and Anopheles arabiensis as 10%, see [25] .

The poisoning effect is modelled with the assumption that at a time instance i, mosquito consumes a dosage of chemical D i spread on the treated net upon contact to the net surface. Thus, the total accumulated dosage C tot is computed as the number of contacts with the net:

where D i is non-zero in case of hitting the net surface (i.e., equal to the unit dosage), and zero otherwise.

The insecticidal-induced increase in mortality is then modelled as:

where the effective poisoning impact is obtained by a scaling coefficient μ p which depends on the given insecticide used for LLIN treatment.

So the total probability of death per unit change in time Δ t is modelled as the sum of natural and insecticide-induced mortality:

Repellency is modelled with the logistic curve multiplied with the repulsion intensity parameter r:

where d ( x , x h ) denotes the distance from the mosquito to the protected human and r ranges from 0 to 1. The parameters d 50 and s determine the range of coverage and the spread of the chemical. The logistic function is modified such that the rejection probability at the candidate position x amplifies as the mosquito approaches the source of repellent. Considering the properties of modern insecticidal treatments [26], the spatial range of the repellent s is taken to be small such that the impact is only within the vicinity of the net.

The repulsion by LLIN is computed in two stages. First, the accept/reject step is applied, where the probability of rejection is given by a logistic function describing the contact irritancy caused by the chemical, as given in Eq. 10. Next, the physical net barrier is taken into account, for which the probability of being blocked by the net is assigned as p net < 1 such that there is a non-zero chance for penetration.

Model extensions

Motion and host-seeking: excito-repellency

The aim is to keep the host-seeking model as minimalistic as possible, by including only the indispensable factors listed in Table ​ Table1. 1 . It turns out, however, that the impact of different chemicals could not be fitted by the basic formulation given above. For instance, the model has to reproduce cases of higher exit and lower contact rates along with more than twice higher mortality rate for An. gambiae than An. arabiensis, following the data reported in [19]. Three new features necessary to characterize the impact of different chemicals on mosquitoes: metabolic detoxification [27, 28], delayed impact [19] and excito-repellency (or insecticide-induced exiting) [19, 29], are introduced. In order to account for insecticide-induced exiting, a scaling factor which not only depends on distance but also on repellent effect is further obtained. The inclusion of both distance and repellent effect is essential in order to properly fit the exit rates, as it accounts for generally higher exit rate when confronted with the treated nets as compared to the control case with the untreated nets. Thus, an excito-repellency parameter [29], μ e is introduced, which depends on the mosquito species and the insecticide utilized in treating a given LLIN, parameterized as:

where C tot denotes the total dosage of chemical consumed by the mosquito (see Eq. 7).

The other two included features: metabolic detoxification (see Eq. 12) and delayed mortality are explained next.

Poisoning and death: detoxification and delayed death rate

Here, the scenarios in the datasets from [19], where An. arabiensis is revealed to have consistently higher (or equal) feeding rate than An. gambiae but considerably lower death rate, are accounted for. These scenarios are inconsistent with the mechanism of the model presented in [20]. The inconsistency is explained by the fact that it is not possible to have simultaneously high feeding rate and low mortality rate if both the probability of death and that of successful feeding is proportional only to the number of contacts with the net. A number of probable reasons can be offered to account for the conflicting situation. One explanation is that the rate of poisoning is different for the two species because it takes time for the poison to get from the salivary glands to the neural system of mosquito and this time delay is suspected to be different for the two mosquito species. However, a large dosage is equally lethal for both An. gambiae and An. arabiensis and mosquitoes do not acquire the lethal dosage upon a single contact with the net but rather a sub-lethal dosage [19]. So, the explanation of detoxification is followed such that the chemical concentration is exponentially decaying with a rate α which depends on the chemical and mosquito species [28, 30]. Hence, given the previous dosage of the chemical C tot ( n ) at the step n, the dosage at the next step n + 1 is calculated by modifying Eq. 7 as:

Additionally, the delayed mortality that is a result of the prolonged impact of poison in mosquitoes is considered. Since poisoning effect is primarily associated with contact with the treated surface, some time is needed for the chemicals to penetrate and reach their target, which in turn depends on the physiological characteristics of the mosquito, such as the sensitivity of target proteins and the thickness of the cuticle [27]. Also, due to enzymatic detoxification, the knock-down time is prolonged. Owing to the high exit rates reported in [19], it was concluded that the mortality induced by the insecticides occurred only after a delay. Although the mosquitoes respond differently with different chemicals, the detailed modelling is spared and the enhanced probability of death is simply taken into account only after a 24-h time period as given by Eq. 6, with Δ t = 24 · 1800 .

The improved model of the chemical-induced exiting and mortality introduced can be calibrated for all the different treatment kits data from [19] (see Additional file 1 for a summary of the datasets). The model is capable of reproducing, e.g., the experimentally recorded lower contact rates along with more than twice higher mortality rates for An. gambiae as compared to An. arabiensis [19]. The calibration is performed using Bayesian sampling methods (adaptive MCMC) in the same way as in [20], more the details are given in Additional file 1. The motive of the MCMC simulations is to find the posterior distibutions of model parameters, that is ‘all’ parameter combinations that reproduce the measured data, within the accuracy given by the estimated error bounds of the data. While most of the parameters are reasonably well identified, some of them are clearly correlated. For instance, as the chemically enhanced mortality rates are now explained by both detoxification and exito-repellency, the respective parameters are strongly correlated with μ p , the earlier introduced death rate coefficient.

Household-scale simulations: household size effect and behavioral alterations

Here, the description of the household and community level modelling is presented, adding more details to the preliminary demonstration given in [20]. First, the ABM of mosquito host-seeking behaviour is extended to the household level with multiple individuals sleeping under the same roof. Next, the modelling is extended to community-level scenarios with several households located in the landscape of interest. See the illustration of the workflow in Fig. ​ Fig.1 1 .

A significant correlation between malaria reduction and the decline in typical household size in malaria-endemic countries is discussed in [2–4]. It was concluded that the larger the number of people sleeping together in non-segregated quarters, the higher the probability of transmitting the infection to new uninfected humans [2]. In Finland, for instance, the probability of malaria disappearance increased when the average number of individuals in one household declined below the threshold of four people, even when no specific control measures were applied [2, 4].

Naturally, there are several other household-related factors that can influence the rate of transmission apart from the household size. Such factors include, e.g., household practices like livestock/poultry rearing, as well as the rate of hygiene maintenance in a given household [31]. For simplicity, these factors are omitted here. The interest of the ABM simulations is in the impact of LLINs. The mosquito density m which can be impacted by these omitted factors, is taken into account in the ODE model. Also, the situation is restricted to a given number of persons sleeping together in the same room, while the approach can be extended also to cases of many people sleeping in separated quarters. The aim here is to demonstrate how household-level factors can be included in ABM simulations, and how even most rudimentary considerations impact the modelling outcomes.

On entering a household, the mosquito is assumed to choose one of the hosts randomly. After this, the modelling reduces to the previous case of a single host in the hut. A few changes are needed, however. The tendency of mosquitoes to switch to neighbouring individuals after spending a certain time in unsuccessful attempts to feed on a protected host, should be considered. Thus, an additional parameter, t max host , the maximal time spent while attempting to feed on a protected host, is introduced [32, 33]. In the absence of more specific knowledge, the parameter is set to 10 min. In addition, same as in the hut-level experiment, mosquitoes are restricted to a maximum host-seeking time, t max inside the household after which they switch to a random walk with no influence of the human bait. Another difference is easier exit from a usual household compared to that from the special design of experimental huts. A typical human dwelling [34] is mimicked by setting the probability of exit to constant value that produces about 90% exit rates per night in the absence of chemical treatment.

Infection with malaria parasites has been shown to alter the behaviour of mosquitoes, with varying effects that are based on the life stage of the parasite [35]. The underlying mechanisms that engender these behavioural alterations are not fully explored but mostly result from at least two manipulation processes. Firstly, the parasite increases the mosquito’s motivation to continue a meal after interruption, thus increasing its probability of taking several bites. Secondly, the parasite impairs the vector’s ability to obtain a full blood meal upon a single bite, inducing the vector to bite several times before it is fully engorged [36, 37]. These behavioural changes associated with infection seem likely to be an evolutionary mechanism that has been developed by malaria parasites, which enhances the spread of infection [38, 39]. A more profound understanding of the behavioural tendencies of parasite-infected mosquitoes alongside the stage-specific changes in their host-seeking behaviour could provide a potential target for genetic manipulation of mosquitoes, as a preventive measure for the elimination of malaria infection [40].

In the simulation, the impact of multiple biting typical for infected mosquitoes is accounted for. Both infected and uninfected mosquitoes are assumed to have the tendency of feeding on multiple hosts [41]. However the tendency of multiple feeding is higher for infected mosquitoes. Thus, the statistics from [36] is employed, which indicate that 10% of uninfected and 22% of infected mosquitoes obtain a blood meal on at least two hosts, while assuming that the maximal number of successful feeding attempts can be up to 5 for both, depending on the accessibility of the hosts. The dosage of blood sufficient for ovipositing is assumed to be achieved after the maximal number of successful feeding attempts is reached. Note that the hut-level data, with one person in the hut, does not contain information on the alterations in behaviour during the host-seeking, so at this point, the literature is relied on. On the other hand, the conjecture that humans infected by the parasite attract more mosquitoes [42] is not included in the simulations, since the hypothesized enhanced attractiveness has demonstrated insignificant impact on the outcome of the simulations (see [43]).

Note that in the simulations, the model parameter values are also re-sampled from the estimated parameter posteriors at each successive iteration of the algorithm to account for parameter uncertainty (see Additional file 1). The main model parameters are summarized in Table ​ Table3 3 .

Table 3

Summary of the basic agent-based model parameters, [20]