Spatial Dynamics of Malaria Transmission

The Ross-Macdonald model has exerted enormous influence over the study of malaria transmission dynamics and control, but it lacked features to describe parasite dispersal, travel, and other important aspects of heterogeneous transmission. Here, we present a patch-based differential equation modeling framework that extends the Ross-Macdonald model with sufficient skill and complexity to support planning, monitoring and evaluation for Plasmodium falciparum malaria control. We designed a generic interface for building structured, spatial models of malaria transmission based on a new algorithm for mosquito blood feeding. We developed new algorithms to simulate adult mosquito demography, dispersal, and egg laying in response to resource availability. The core dynamical components describing mosquito ecology and malaria transmission were decomposed, redesigned and reassembled into a modular framework. Structural elements in the framework – human population strata, patches, and aquatic habitats – interact through a flexible design that facilitates construction of ensembles of models with scalable complexity to support robust analytics for malaria policy and adaptive malaria control. We propose updated definitions for the human biting rate and entomological inoculation rates. We present new formulas to describe parasite dispersal and spatial dynamics under steady state conditions, including the human biting rates, parasite dispersal, the “vectorial capacity matrix,” a human transmitting capacity distribution matrix, and threshold conditions. An R package that implements the framework, solves the differential equations, and computes spatial metrics for models developed in this framework has been developed. Development of the model and metrics have focused on malaria, but since the framework is modular, the same ideas and software can be applied to other mosquito-borne pathogen systems. 4 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint

1 Plasmodium falciparum transmission dynamics are complex: they involve 2 multiple-agents, non-linear dynamics, localized spatial interactions, spatial, temporal 3 and behavioral heterogeneity, stochasticity, and exogenous forcing by weather, 4 hydrology, and malaria control. Over time, these processes can be modified by economic 5 development; by changing socioeconomic status, human incentives and social norms; 6 and by the evolution of resistance. Every one of these features of malaria transmission 7 dynamics and control presents its own set of challenges to the quantitative study of 8 malaria for scientific research and for analytics to support policy. An important 9 practical problem is how to quantify and synthesize all of the factors affecting 10 transmission at some particular place and time to support malaria control programs in 11 various ways, including monitoring and evaluation of malaria control. The study of 12 complex spatial processes are best addressed using some sort of mathematical model. 13 Here, to fill a need to give robust policy advice, we have developed a modular 14 framework with accompanying software to build and analyze suites of models with 15 scalable complexity for malaria spatial transmission dynamics and control. 16 A starting point for the quantitative study of malaria transmission dynamics has 17 been the Ross-Macdonald model, which played a central role in developing basic theory 18 and metrics for malaria [1,2]. That model is simple, general, and conceptually useful, 19 but it is not realistic enough to describe many important features of transmission [3]. 20 The model’s lack of realism has also limited its applicability: simple models support 21 generic policy advice, but specific advice – tailored to context – must be based on 22 models that can quantify and weigh the effects of locally relevant details [4]. A basic 23

of the Ross-Macdonald model was that it lacked features required to describe 24 spatial transmission dynamics and control. Mathematical models for spatial dynamics 25 of mosquito-borne pathogens have been developed [5–18], but there is a need for a 26 generalized synthetic framework to develop and use spatial dynamic models, to extend 27 the Ross-Macdonald model to define and analyze parasite dispersal, to define and 28 measure malaria connectivity [19], and to link spatial dynamics to spatial data. The 29 Ross-Macdonald model is also missing other features that are relevant for malaria 30 dynamics and control, which can be identified from a survey of studies that have 31 modeled mosquito-borne diseases (see Box #1) [2]. Modeling and analyzing real 32 systems can become overwhelming because of computational, parametric, or conceptual 33 challenges that arise from combining all the factors, dimensions, interactions, features, 34 and processes. Individual-based models (IBMs) have been developed around algorithms 35 that make it possible to deal with the complexity by simulating individual states and 36 transitions in silico [20], but these high-dimensional computational approaches have 37 some limitations that limit their use and applicability. IBMs require intensive 38 computation, are challenging to parameterize, are difficult to critically evaluate, and 39 their output that is often as difficult to analyze and understand as malaria itself. Using 40 a modular framework, we present an alternative way of dealing with the complexity that 41 is analytically tractable, including some new algorithms to understand mosquito ecology, 42 parasite transmission by mosquitoes, and parasite dispersal on spatial landscapes. 43 In most places, malaria transmission has been modified by control. The extent of 44 effect modification by malaria control is occasionally revealed when health systems are 45 5 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint disrupted (e.g., [21,22]), when malaria control is relaxed or abandoned [23], or when 46 resistance evolves to drugs or insecticides (e.g., [24, 25]). Programs must weigh evidence 47 and make decisions through analysis of counterfactuals, rather than through direct 48 estimation of control effect sizes, since there would be drastic consequences to 49 experimentally disrupting control. A predominant need in most contexts is thus a set of 50

to quantify transmission in its local context as a baseline that has been 51 modified by control. A challenge to achieving this has been that the responses to 52 control efforts are context dependent and have been highly variable across settings. 53 Relevant factors affecting responses to control include details about blood feeding, 54 mosquito ecology, and mosquito behaviors that affect contact with interventions ( e.g., 55 resting indoors and IRS). To reconstruct the counterfactual baseline, transmission must 56 be understood in terms of innate mosquito behaviors responding to local resources, 57 vector control, and other contextual factors that have been modified by control. All 58 these have been characterized as being notoriously context dependent and 59 heterogeneous [26–28]. What are the local factors that determine baseline malaria 60 transmission, effect modification, and differences in effect modification at some 61 particular place and time? Basic concerns about the heterogeneous impacts of vector 62 control raise a larger set of questions about how to study and quantify transmission in a 63 way that is relevant for planning malaria control. 64 This new framework is thus an attempt to bridge two well-established but somewhat 65 contradictory views of malaria. One view is that human malaria transmission dynamics 66 and control are so moulded by local ecology and other conditions that the factors 67 driving transmission or responses to control at one time and place are unlikely to hold 68 elsewhere [27]. Another view – encouraged by the rigorous analysis of the 69 Ross-Macdonald model and extensions of it – is that malaria transmission intensity can 70 be quantified using a small set of bionomic parameters to compute basic reproductive 71 numbers, which also provide a basis for computing threshold conditions for endemic 72 malaria. To build a bridge, the contextual factors affecting basic bionomic parameters 73 must be identified and integrated with new theory describing spatial extensions of the 74 basic metrics, including rigorous, quantitative description of parasite dispersal, and 75 some estimates of the appropriate spatial scales to measure malaria transmission [3]. 76 Context-dependency is an uncomfortable but unavoidable fact of malaria ecology. 77 The heterogeneous nature of transmission and the causes and consequences of variable 78 responses to control have been a difficult and sometimes contentious problem for 79 scientists studying malaria, for national malaria programs and funding agencies making 80 malaria policy, and for malaria advocates. Historical trends in malaria and the outcomes 81 of malaria control have been so variable that case studies can be found to support rosy 82 projections, alarmist warnings, or contradictory claims about the underlying causes of 83 trends or patterns. To be useful, studies of malaria and programmatic evaluations must 84 acknowledge the important role of context, the multi-factorial nature of causation in 85 these complex systems, non-linear responses to control, the difficulty of measuring 86 heterogeneous systems, and the resulting uncertainty. A consequence of context 87 dependency is the difficulty in drawing conclusions that generalize across systems. 88 The framework is designed to support development of robust malaria policy advice 89 and to find practical ways of dealing with uncertainty. While scientific research and 90 policy analytics grapple with the same issues and use similar methods, they often put 91 6 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint very different weights on uncertainty. Uncertainty affects the ability to do effective 92 inference for scientific research versus policy analytics – questions about what is known 93 versus what should be done. To address these concerns and give policy advice despite 94 uncertainty, an integrated inferential framework is needed to weigh evidence, integrate 95 the effects of multiple exogenous factors (often involving experts from distinct 96 specialties), estimate their effect sizes, quantify uncertainty, and identify critical gaps. 97 Statistical theory and inferential methods have been developed around the principle of 98 parsimony for scientific inference, but substantially less attention has been given to 99 appropriate designs for analyses that can give advice that is robust to uncertainty. Are 100 the conclusions of an analysis robust to reasonable alternative formulations of a model, 101 and how well are policy recommendations really supported by the evidence? Concerns 102 about robustness could lead to study designs that make different tradeoffs between 103 realism and abstraction. For example, compared with parsimonious models, models 104 with a high degree of realism might be more useful for identifying critical missing data 105 and prioritizing studies to collect it. Robust analytics requires having a modeling 106 framework to build suites of models that are realistic enough to weigh the importance of 107 the major drivers of transmission despite major knowledge gaps. 108 To address these needs, we have developed a new, modular framework designed to 109 support development of models for robust, simulation-based analytics and adaptive 110 malaria control with scalable complexity. With scalable complexity in model building, 111 members of a model ensemble could range from very simple to very complex, and that 112 models along that spectrum are related to one another through a logical sequence of 113 structural or parametric changes. At one extreme, this framework includes the 114 Ross-Macdonald model, a simple system of differential equations describing the parasite 115 life-cycle in mosquito and vertebrate host populations linked by transmission during 116 blood feeding [1,29,30]. By extending the Ross-Macdonald model, simple models can be 117 extended step by step to add complexity or heterogeneity that could be important – 118 based on a priori considerations – yet difficult to quantify or poorly informed by 119 existing data (Box #1). With modularity, it is possible to develop new dynamical 120 systems models describing some parts of the system, add or modify components, or add 121 a set of exogenous factors that force a system. It is also relatively straightforward to 122 modify functional responses, or to modify some basic parameters affecting the outcome. 123 Swarms of models can thus be developed to analyze data and to test the robustness of 124 any conclusions. To demonstrate scalable complexity, we here present a complicated 125 family of models that has terms and variables anticipating modification by weather or 126 malaria control. For practical reasons, the model family we present here was scaled back 127 to include a limited set of elements describing transmission, but leaving in place the 128 elements that facilitate modeling control (Box #1). The resulting extensible framework 129 that is capable of describing and analyzing malaria spatial transmission dynamics and 130 control with a high degree of realism in any particular setting. An R package which 131 implements the modular differential equations and spatial metrics presented in the 132 article is available with documentation (Supplement 1 - 133 https://dd-harp.github.io/exDE/). Despite being programmed in R, the 134 implementation of the mathematical framework into code should be easily adapted to 135 any high-level programming language. 136 In F ramework, we first present the modular concepts and structural elements, 137 7 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint including a new blood feeding model. Next we present one exemplar model family for 138 each dynamical component. In Spatial Metrics, we develop a set of metrics that 139 describe various aspects of parasite spatial dynamics, including metrics for parasite 140 dispersal, connectivity, and the parasite’s reproductive success. Finally, in 141 Quantifying T ransmission in a Place , we discuss the application of these models to 142 the investigation of malaria transmission dynamics and control in a particular place. 143 8 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Box 1: F eatures This generalized, modular framework presents equations integrating multiple agents and interacting processes. Many of these innovations appeared first elsewhere, but here they are integrated into a single framework: • Immature mosquito population dynamics structured in distinct aquatic habitats linked to adult populations through egg laying and emergence [31,32]; • Spatially heterogeneous blood feeding and parasite mixing on vertebrate populations (i.e., blood hosts) with dynamically changing availability, such that feeding rates and the human fraction change adjust to changing conditions [33–36]; • Heterogeneous adult mosquito behaviors, including dispersal, survival, blood feeding, egg laying, mating, and sugar feeding on landscapes in response to spatially heterogeneous resource availability ( e.g., mating sites, sugar sources, blood hosts, aquatic habitats) [37–39]; • Multiple vector species or types with different host preferences, daily activity patterns, habitats, etc. [40], and potentially with inter-specific resource-based competition in habitats; • Human mobility based on a concept of time at risk, which combines time spent by humans in places where they are at risk with mosquito blood feeding activity, preferences and other factors [9,18]; • The capability to model indoor and outdoor spaces for blood feeding, exposure, and vector control; • A non-linear relationship between the daily entomological inoculation rate (EIR) and the daily force of infection (FoI) due to heterogeneous exposure [41]. • Malaria importation through multiple routes [42]; • An exogenously forced, time-varying extrinsic incubation period (EIP) to model effects of temperature on parasite development; The model has flexible structural elements to stratify an area into patches, to model any distribution of aquatic habitats, and to stratify a human population into sub-populations by age, immunity, or any heterogeneous, epidemiologically relevant trait. The software also includes time-dependent terms and structures to model exogenous forcing by weather, modification of exposure or transmission by vector control in relation to coverage, including effects of spatial repellents and mosquito behaviors that result in heterogeneous local contact patterns with vector-based interventions. 144 9 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Framework 145 To describe malaria spatial dynamics with scalable complexity, we designed a modular 146 framework for model building around four core dynamical components, each one a 147 (potentially non-linear) state-space model. An interface rigidly defines interactions 148 among those components, based on passing terms we call dynamical quantities. All state 149 variables are vectors of arbitrary length, to accommodate models with different 150 structure or spatial granularity. 151 To model mosquito ecology, we consider immature mosquitoes in a set of aquatic 152 habitats, and adult mosquitoes in a set of patches. A state space model describes 153 aquatic immature mosquito populations (L) with dynamics dL/dt requiring an input 154 term from adult mosquito populations: the daily rate eggs are laid in each habitat (η ). 155 A coupled state space model describes mature adult female mosquito populations ( M) 156 with dynamics dM/dt requiring an input term from the aquatic mosquito populations: 157 the rate adults emerge from all the habitats in each patch (Λ). A state space model for 158 parasite infection dynamics in mosquitoes (Y, which extendsM) with dynamics dY/dt, 159 requires an input term from human malaria epidemiology: the net infectiousness of 160 humans (NI), the probability a mosquito becomes infected after blood feeding on a 161 human (denoted κ). A state space model describing parasite infection dynamics in 162 humans, immunity, and disease (X ) with dynamics dX /dt, requires an input term from 163 adult mosquito infection dynamics: the daily EIR ( E). The inputs to one component 164 can be passed as trace functions or as the outputs of another coupled component, which 165 is called the interface of each dynamical component; a generic interface is coded for 166 each term and if needed specialized methods can be written for particular models. 167 Models in the framework have the following form: 168 dL/dt = FL (η,L) dM/dt = FM (Λ,M) dY/dt = FY (κ,M,Y) dX /dt = FX (E,X ) . (1) The interactions among these dynamical components are thus defined by four input 169 terms (η, Λ, κ, and E), which may be computed as outputs of another component or 170 provided as an external forcing term (Fig. 1). Because these terms can be computed 171 from the state of the model and are used to couple different model components together, 172 we call these dynamical quantities. These terms are rates which determine how 173 components interact (e.g., flows between components). Because construction of these 174 dynamical quantities can be done in a generic way, computation of these quantities in 175 code can be done for any model which fulfills the interface of its dynamical component. 176 The dynamical quantities responsible for transfer of pathogens between hosts and 177 vectors are E and κ, the EIR and NI of humans, respectively. These quantities couple 178 the dynamics between the human X and mosquitoY dynamical components. To allow 179 computation of E and κ to be highly generic across various types of models of human 180 and mosquito infection, we developed a new model of blood feeding which produces β, 181 the biting distribution matrix describing how bites arising from mosquitoes at patches 182 are taken on human population strata. 183 10 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Similarly, the adultM and aquaticL mosquito components are coupled via egg 184 laying from adults in aquatic habitats, and emergence of new adults from those aquatic 185 habitats. Because the patches where adult mosquitoes are found may contain many (or 186 no) aquatic habitats, another matrix translates the rate of egg production from adults 187 into egg deposition in each aquatic habitat η. Likewise, each aquatic habitat produces 188 newly emerging adult mosquitoes at some rate α, which in general depends on the 189 current aquatic population, and therefore on lagged adult densities. Another matrix 190 maps this into the rate at which new adults are added to each mosquito population, Λ. 191 In addition to reformulating blood feeding and egg laying, the framework includes 192 mathematical descriptions of survival, search for blood hosts or habitats, and dispersal. 193 These new models of adult mosquito behaviors have all been reformulated around the 194 concept of heterogeneous resource availability and functional responses to available 195 resources. 196 Fig 1. Models for malaria transmission dynamics are naturally modular (see Eq. 1). The dynamic modules describe a stratified human population (purple) that interacts through blood feeding (red) with adult mosquito populations in a discrete spatial domain; each patch could contain a set of aquatic habitats. Two components, L andM, describe mosquito ecology: dynamics of immature mosquitoes (blue) in aquatic habitats are described by a system of equations dL/dt; and dynamics of adult mosquitoes (green) are described by dM/dt. Habitat locations within patches are described by a membership matrix,N . Eggs hatch into larval mosquitoes, that develop, pupate, and later emerge from habitats as mature adults ( α) and added to the adult populations in each patch (Λ). Adults lay eggs ( ν), which are distributed spatially according to which patch habitats belong (N ). Egg deposition rates at the habitats are ( η). Two additional components,Y andX , describe parasite infection dynamics and transmission: that for mosquitoes, described by dY/dt and in humans, described by dX /dt, are linked through parasite transmission. A new model for blood feeding describes how blood meals are allocated among humans ( β) and associated parasite transmission rates: the density of infectious humans by strata ( X) is used to compute net infectiousness (NI) of humans to mosquitoes in patches ( κ); and the density of infectious blood feeding mosquitoes ( Z) is used to compute the entomological inoculation rate (EIR) on each strata ( E). The modular framework was implemented as a software package in R [43] 197 (Supplement 1 is the website https://dd-harp.github.io/exDE/). The software 198 builds dynamical models of malaria in a modular way using method dispatch to define 199 generic code which implements the framework described here. The dynamical models 200 are functions which return arrays of derivatives of state variables, and can be solved 201 using the integrators available in deSolve, or other tools in R [43,44]. The software also 202 11 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint includes routines that compute steady state conditions and spatial metrics (see Spatial 203 Metrics, below). Because each component has an interface – the generic functions that 204 compute and pass of dynamical quantities between components – any new model can be 205 implemented which fulfills a specific interface, independent of the rest of the framework. 206 In this way, building and testing new models of particular components is straightforward, 207 and the framework is flexible and extensible. As new models are required, they will be 208 added to the package, increasing its applicability and scope over time. 209 We have developed a glossary of terms (Supplement 2). In the equations that follow, 210 for each dynamical component, we describe the model structure in detail, and we 211 present one family of models describing transmission dynamics in a single vector species. 212 In Supplement 3, we formulate a model using both conventional notation and the 213 modular notation of this framework. In Supplement 1, we have implemented a 214 previously published model of malaria transmission on Bioko Island [45]. In Supplement 215 4, we extend the discussion of vector dynamics, including a discussion of models with 216 multiple vector species. All the terms and parameters may be time dependent to 217 accommodate seasonality or modification by exogenous factors: seasonal travel, 218 exogenous forcing by weather, and parameter modification by vector control. Analysis 219 of temporal heterogeneity in this same framework is outside the scope of this study, it 220 but would be straightforward extension following approaches analogous to those shown 221 in the supplements. 222 Box 2: Notation Equations describing spatial processes include terms describing scalar quantities, vectors of scalars, vectors of functions, and matrices. We have avoided using any notation to designate a vector or indicate it could be time-dependent, in part, because it would be ubiquitous; most parameters could vary by space and time. The most general form of a term or parameter is usually described when it is first presented, but most terms describing a vector or matrix should be assumed to be modifiable. In writing out the equations, we consistently use the center dot, “ ·”, in equations to denote the dot product of two matrices, or a matrix and a vector. The juxtaposition of two vectors denotes element-wise product, and 1/∗ denotes the vector of the inverses of each element. The symbol ⊙ denotes the Hadamard product (i.e., element-wise multiplication) of two matrices. When x is a vector, diag(x) is a matrix with the elements of x on the main diagonal. The identity matrix is denoted I, and 1 denotes a row or column vector with each element equal to 1. When F is a functional response, we assume it accepts vector arguments and returns a vector of the same length, i.e.,|F (X)| =|X|. The glossary (Supplement 1) discusses the dimensions of each term. 223 224 12 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Model Structure 225 The following describes, in detail, the structural elements and the algorithms that 226 connect them. Adult mosquito and human population strata are connected through 227 blood feeding and transmission, and adult and aquatic mosquito populations are 228 connected through egg laying and emergence. 229 Structural Elements The framework has been designed to build model ensembles 230 with the goal of studying the spatial transmission dynamics of malaria in a defined 231 geographical area, called the spatial domain. An important part of this framework is 232 having flexibility in defining the model structure to describe spatial and population 233 heterogeneity at the appropriate level of detail, depending on the needs of a study and 234 the available data. The structural elements – the patches, the aquatic habitats, and the 235 population strata – were designed to handle arbitrary patch definitions, arbitrary 236 human population residency patterns and stratification, and arbitrary numbers and 237 locations of aquatic habitats. 238 To deal with spatial heterogeneity in transmission, we subdivide the spatial domain 239 and identify a set of p patches that includes all locations relevant for studying and 240 quantifying mosquito ecology or transmission: places where people live; places where 241 mosquitoes blood feed; or places with aquatic habitats where mosquitoes lay eggs. We 242 assume that there are l aquatic habitats with actual physical locations that are nested 243 within the patches. To deal with heterogeneity in the human population, the model 244 accommodates stratification. The human population is sub-divided into a set of n 245 population strata by residency, immunity, behaviors affecting risk, or any other 246 epidemiologically relevant factors (Supplement 5). Human populations are assigned a 247 single residency patch, where they live and spend most of their nights. Other 248 subdivisions of the human population could take into account age, sex, travel patterns, 249 ITN usage, or any trait that is heterogeneous and epidemiologically relevant. The total 250 census population size, the number of people who reside in each patch in the spatial 251 domain, is given by a vector denoted P (of length p). The number of people in each 252 stratum is given by a vector H (of length n). In this model, it is not necessary for every 253 patch to have some residents. 254 To manage terms for interactions among structural elements, we create two 255 mathematical objects called membership matrices that aggregate quantities to patches 256 (Supplement 3). Since the l aquatic habitats are scattered among the patches, we define 257 the habitat membership matrix N , a p× l matrix, that aggregates quantities from the l 258 aquatic habitats to p patches where they are found. Similarly, we define the strata 259 membership matrixJ , a p× n matrix, that aggregates the n human population strata 260 to the p patches where they reside. The census population size, for example, is 261 P =J· H. If a human population were stratified by other traits, such as frequent travel 262 or age, a membership matrix could be created to aggregate model output by trait. 263 The framework has also been designed to accommodate models with multiple 264 mosquito vector species or types (see Supplement 4). Most of the following discussion 265 assumes there is just one vector species, but we point out where the framework has can 266 generalize to multiple vector species. 267 13 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Human Mobility After defining the model structure (i.e., the patches and 268 population strata), the next challenge is to construct the algorithms describing local 269 human mobility and travel. Local mobility determines where and when humans are 270 available and exposed to blood feeding mosquitoes within the spatial domain. We define 271 travel in this model by time spent outside the spatial domain; travel and mobility are 272 thus different modalities and handled with different constructs. 273 To model local human mobility patterns within the patches, we develop a model 274 describing the fraction of time spent by humans in each stratum among the 275 patches [9,18]. The information is summarized in a time-dependent p× n matrix Θ(t), 276 called the Time Spent (TiSp) matrix (Supplement 5). Each column in a TiSp matrix 277 describes the fraction of time spent in each patch by an individual from a single 278 stratum. In formulating the TiSp matrix, we account for time spent by time of day in 279 the patches where mosquitoes are blood feeding. Total time spent should subtract time 280 spent traveling and and time spent in the spatial domain in places where there is no risk 281 (e.g., in office buildings). 282 Blood feeding combines human and mosquito behaviors. Since mosquito blood 283 feeding has a daily rhythm [46], time at risk modifies time spent to account for 284 differences in mosquito daily blood feeding activity rates. We let ξ(t) denote a 285 species-specific circadian weighting function for blood feeding rates, constrained such 286 that ∫1 0 ξ(t)dt = 1, which appropriately assigns a weight to time spent by time of day 287 (Supplement 5). Using ξ, we compute the Time At Risk (TaR) matrix as time spent 288 weighted by mosquito activity: Ψ( t) = diag (ξ(t))· Θ(t). 289 This distinction between TiSp and TaR matrices makes it possible to study human 290 mosquito contact in detail, to quantify differential transmission by multiple vectors with 291 the same human mobility patterns, and to quantify other aspects of mosquito-human 292 contact [47,48]. A model could have two or more vector species, each with different 293 blood feeding patterns ( ξ1 and ξ2), so that one TiSp matrix would be transformed into 294 two different TaR matrices (Ψ1 = ξ1Θ and Ψ 2 = ξ2Θ). 295 Denominators and Availability After defining host population movement, it is 296 necessary to compute appropriate denominators to model blood feeding, based on the 297 models for time spent and time at risk. Because of mobility, mosquito preferences, and 298 human behaviors, the denominators for blood feeding are different from the resident 299 population size – the number that would be used by most studies (Fig. 2). 300 An important intermediate quantity is ambient population density, which describes 301 the population present in patches at a point in time. In a mobile population, the 302 ambient population density will tend to be different from resident population density. 303 From the time spent matrix, the ambient population density is a vector of length p 304 given by: 305 A(t) = Θ(t)· H. (2) Similarly, ambient population density at risk is given by: Ψ( t)· H. One way to 306 understand what the TiSp matrix means is by taking ratios of ambient to resident 307 populations. The ambient density of residents is Ar = (J⊙ Θ)· H, where⊙ denotes the 308 Hadamard (element-wise) product. The non-resident, non-visitor, ambient population is 309 A− Ar. The ratios of various census and ambient population densities (e.g., the ratio of 310 14 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Fig 2. Denominators and Mixing A schematic diagram relating various concepts of population density under a model of human mobility, resulting in a biting distribution matrix, β. Here, and and in Figures 3-6, rounded rectangles denote endogenous state variables, sharp rectangles denote endogenous dynamical quantities, and parallelograms represent exogenous data or factors. Purple indicates the element is related to human populations, green for mosquitoes, and red for biting and transmission. Population strata (H ) describe how persons are allocated across demographic characteristics. The matrixJ distributes these strata across space (patch), according to place of residency. By combining information on how people spend their time across space (Θ( t)) and mosquito activity (ξ(t)) a time at risk (TaR) matrix Ψ is generated describing how person-time at risk is distributed across space. Because blood feeding can be modified by human and mosquito factors ( e.g., net use and biting preferences), search weights (wf(t)) may further weight person-time at risk. The final result is a biting distribution matrix β, which is the fraction of each bite in each patch that would arise on an individual in each stratum, so diag(H )· β = 1. residents to ambient population P/A, defined wherever A > 0), can be used to 311 understand and diagnose unrealistic terms in a TiSp or TaR matrix. The ambient 312 population thus provides one easy statistic to understand TiSp or TaR matrices. 313 To model the denominators for blood feeding, we also consider other factors – 314 mosquito preferences or human behaviors or traits such as ITN usage – that affect host 315 availability to mosquitoes and relative biting rates on the strata [33]. We assign biting 316 weights, wf, to each strata, where we think of wf = 1 as the value that would be 317 assigned to an average person under baseline conditions ( e.g., without a net). These 318 weights affect both the total biting rates and the relative biting rates on the ambient 319 population. We define the availability of the host populations to mosquitoes for blood 320 feeding as: 321 W = Ψ· wf H. (3) Availability is thus defined in units of weighted person-days at risk, and W is a vector 322 of length p describing total human availability in each patch. 323 We also consider the presence of a population of visitors, a non-resident population 324 spending time in the spatial domain (Supplement 5). We assume that some visitors 325 could be present, and that some of them could be infectious. We can let Aδ denote the 326 ambient density of visitors, but we let Wδ denote their availability by patch. The 327 resident fraction or fraction of human blood meals taken on a resident in each patch, a 328 vector of length p denoted υ, is: 329 υ = W W + Wδ . (4) 15 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint The total availability of humans for blood feeding, in each patch, is thus W + Wδ. 330 Fig 3. Blood F eeding and Human Biting Rates The daily human biting rates (HBR) for the resident population strata are defined as the expected number of bites by vectors, per person, per day. To compute the HBR, we count up exposure over all the patches where residents spend time. We also consider the presence of visitors and other bloodhosts (yellow input), which increases the total available hosts. Blood F eeding With a well-defined population denominator, we can compute the 331 frequency of blood feeding rates and the human fraction (i.e., the fraction of human 332 blood meals among all blood meals) in each patch in response to the availability of 333 humans and other available vertebrate hosts. To do so, we use functional responses to 334 model blood feeding rates and habits [33–36]. 335 Human availability, W , is often highly variable among patches and over time, which 336 could affect the rate mosquitoes blood feed (Fig. 3). Mosquitoes could also feed on other 337 vertebrate hosts. To model blood feeding, we supply a vector of functions describing the 338 availability of non-human vertebrate hosts in each patch over time, denoted O(t). We 339 assume that mosquito preferences could scale with host densities, so we assign a shape 340 parameter, ζ, that modifies how preferences scale with host densities. Total availability 341 of all vertebrate hosts for blood feeding is B = W + Wδ + Oζ (Supplement 5). 342 Let f(t) denote the blood feeding rate, the number of blood meals, per mosquito, 343 per day. To guarantee mathematical consistency in computing blood feeding rates (e.g., 344 if B = 0, then it should be true that f = 0), we can model time-dependent blood 345 feeding rates, where f(t) is a vector of length p, as: 346 f(t) = Ff(B) = fx sf B 1 + sf B . (5) Depending on a shape parameter(s), sf, blood feeding rates increase with host 347 availability up to a maximum (or maxima) fx, which is limited by the time it takes to 348 search, process the blood meal, lay eggs, and perhaps to sugar feed. The fraction of 349 blood meals taken on humans at a point in time, a vector of length p denoted q(t), is 350 called the human blood feeding fraction or human fraction: 351 q(t) = W + Wδ B . (6) The local human fraction, the fraction feeding on resident humans, is thus υq = W/B. 352 The functional forms guarantee that when no humans are present, it must be true that 353 f q = 0; and when only humans are available, it must be true that q = 1. 354 16 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Mixing and Parasite T ransmission The model for mixing is an answer to the 355 question: How are blood meals in a patch allocated among humans in the strata? The 356 time at risk matrix and the factors affecting blood feeding rates and habits in each 357 patch must be consistent with the algorithm that computes the distribution of biting 358 and parasite mixing. 359 To allocate mosquito bites in patches among the resident strata, we let β denote a 360 n× p biting distribution matrix: 361 β(t) = diag(wf)· ΨT· diag ( 1 W (t) ) . (7) Each column of β describes the fraction of a bite in a patch that lands on an individual 362 in each strata, so the matrix diag(H )· β gives the fraction of bites that land on each 363 stratum, and its columns sum to unity. 364 In the models for mosquito ecology and infection dynamics, we define variables 365 (vectors of length p) for the density of mosquitoes ( M) and infectious mosquitoes (Z ). 366 From these, we derive an expression for the daily human biting rate (HBR) and 367 entomological inoculation rate (EIR) for all the strata. The sporozoite rate (SR) in each 368 patch is given by: 369 z = Z M . (8) The net per-capita human blood feeding rates in each patch, or f qM/W, are hereafter 370 called the patch HBR (pHBR), and f qZ/W is hereafter called the patch EIR (pEIR) 371 for infectious mosquitoes. By way of contrast, exposure risk for the strata – the HBR 372 and EIR – are defined as the number of bites / infectious bites by vectors, per person, 373 per day. The HBR is β· f qυM, and the EIR is the product of the HBR and the SR, or 374 E = β· f qυZ. (9) To draw a sharp contrast between the terms, the pHBR and pEIR describe the number 375 of bites / infectious bites, per person, in patches. They are stratified by location, so 376 they are vectors of length p. The HBR and the EIR are stratified quantities that sum 377 exposure over all locations for the strata, so they are vectors of length n. 378 Each model for parasite infection dynamics in humans defines a quantity, x, the 379 probability a mosquito becomes infected after biting a human in each stratum. The 380 quantity X = xH, a vector of length n, is herein called the infective density of 381 infectious human residents. We can also specify the infective density of visitors, Xδ 382 where Xδ = xδWδ is intrinsically using the availability of visitors. The net 383 infectiousness (NI) for the mosquito populations in all the patches, denoted κ, is: 384 κ = υβT· X + (1− υ)Xδ (10) The force of infection for the mosquito population is thus f qκ. 385 Egg Laying To compute quantities affecting mosquito ecology and population 386 dynamics, we need to formulate algorithms to compute egg laying rates and egg laying 387 distributions: how many eggs are laid by adult mosquitoes in a patch, and how are they 388 17 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Fig 4. Egg Laying and Egg Deposition The availability of aquatic habitats (Q) the patch sum of habitat search weights ( Q =N· wν), and the egg distribution matrix (U) describes the locally normalized search weights. Available habitat determines per-capita oviposition rates ( ν) by the population of gravid mosquitoes ( G) in a patch through a functional response to availability, Fν(Q). The net egg laying rate, per-patch, is Γ = χνG. The eggs are distributed among the aquatic habitats (U ) so that the egg deposition rates in habitats is η =U· Γ. distributed among the aquatic habitats in that patch? To do so, we develop the concept 389 of habitat availability. We assign a search weight to each aquatic habitat, wν. Using the 390 patch membership matrix,N , we define aquatic habitat availability as: 391 Q(t) =N· wν(t) (11) For each patch, total habitat availability is the sum of the search weights for habitats in 392 that patch. 393 Daily, per-capita oviposition rates of gravid mosquitoes are computed using a 394 functional response to habitat availability, such as: 395 ν = Fν(Q) = νx sνQ 1 + sνQ . (12) where νx is the highest possible egg-laying rate for a gravid female, and sν is a shape 396 parameter. We note that if Q = 0, then ν = Fν(0) = 0. We let G = FG(M) denote the 397 density of gravid mosquitoes, and we let χ denote the number of eggs laid, per batch. 398 The net egg laying rate, per patch, per day, is: 399 Γ = χνG (13) To model egg distribution among habitats, we formulate an egg distribution matrix ( U) 400 that allocates eggs to habitats in proportion to local habitat availability. To compute U, 401 for computational reasons we first create Q∗ by setting any zero entries to an arbitrary 402 positive value (if Q = 0, then ν = 0, so associated products will later be multiplied by 403 zero), and the egg deposition rate, η, is computed by: 404 U (N , wν) = diag(wν)·N T· diag ( 1 Q∗ ) . (14) Finally, we can compute egg deposition rates in the habitats: 405 η =U· Γ (15) 18 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint While Γ (a vector of length p) describes the net egg-laying rate of the adult mosquito 406 population in each patch, per day η (a vector of length l) describes the number of eggs 407 laid, in each habitat, per day. 408 Core Dynamical Components 409 The dynamical quantities whose computation was described above, are configurable 410 elements that connect the four dynamical components: aquatic mosquito ecology; adult 411 mosquito ecology and infection dynamics; and infection and immunity, including human 412 demography. In the following, we describe one model family for each component, 413 including functions that compute terms required for the dynamical quantities; in code 414 these are the generic interface of each dynamical component. In Supplement 1 415 (https://dd-harp.github.io/exDE/), we have formulated alternative model families for 416 some of the components. 417 Aquatic Mosquito Ecology The first core dynamical component describes aquatic 418 mosquito population dynamics; the algorithm computes mosquito survival and 419 development from eggs laid through adults emerging. For aquatic population dynamics, 420 we here adapt a previously published model [31,32]. 421 Let L(t) denote the total density of immature mosquitoes. We let ψ(t) denote 422 maturation rates, φ(t) the density independent mortality rate, and θ(t)L(t) describes 423 increased per-capita mortality due to mean crowding. The aquatic dynamics are thus: 424 dL dt = η− (ψ + φ + θL)L (16) The total emergence rate of female mosquitoes in this model, per aquatic habitat, is: 425 α(t) = Fα (L (t)) = ψ(t)L(t) 2 . (17) These are recruited into the adult population in the patch, so that the net emergence 426 rate per patch is: 427 Λ(t) =N· α (18) While α is a vector of length l, Λ is a vector of length p. This is passed as input to the 428 equations describing adult populations (below). 429 Given uncertainty about the factors affecting immature mosquito populations, we 430 assume studies might choose to formulate and analyze alternative dynamics. Other 431 dynamical systems models for aquatic ecology in the framework are defined by state 432 variables,L, with dynamics defined by a system of equations dL/dt = η− FL(L), and a 433 function such that α = FΛ(L), such that Λ =N· α (Supplement 4). 434 Adult Mosquito Ecology The second core dynamical component describes adult 435 mosquito ecology. Given all the functions, terms and parameters above, we have 436 formulated a set of algorithms describing adult mosquito mortality and dispersal that 437 are internally consistent. All this is embodied in the mosquito demographic matrix, 438 called Ω(t). 439 19 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Fig 5. Adult mosquito demography is defined by survival and dispersal. Mobility rates and dispersal are determined by the available of resources: aquatic habitats (Q), available humans (W + Wδ) and other blood hosts ( Oζ), and sugar ( S). The emigration rate is a functional response ( Fσ) that increases if any one of the resources is missing. Resource availability and distance also play a role in computing the dispersal kernel, K, that determines where mosquitoes land if they leave a patch. When combined with mortality, a matrix Ω is produced which describes the behavior of adult mosquitoes after emergence. We assume mosquito mobility is driven by a search for resources. We have already 440 defined total blood host availability B, and aquatic habitat availability Q. We also 441 consider sugar availability, S(t), which is passed to the model as a function vector of 442 length p. We assume mosquitoes leave a patch while searching for resources, and that 443 they leave a patch more frequently if the resources are less available. Patch-specific 444 emigration rates, σ(t), are a functional response to resource availability: 445 σ = Fσ(B, Q, S) = σx ( σB 1 + sBB + σQ 1 + sQQ + σS 1 + sSS ) (19) The parameters σB, σQ, and σS determine the rate that mosquitoes leave a patch if no 446 resources are available, and the shape parameters sB, sQ, and sS determine how the 447 rate of patch leaving is reduced by the availability of resources. The shape parameter σx 448 is a scaling parameter that can be used to adjust models with differing patch sizes. 449 Similarly, we formulate a mosquito dispersal matrix, K(t) that describes where 450 mosquitoes land after they leave each patch (the diagonal elements of K are constrained 451 to be equal to zero, Supplement 4). 452 We let g(t) denote the local per-capita mortality rate of mosquitoes in each patch. 453 The matrix Ω(t) describes adult mosquito survival and dispersal: 454 Ω = diag(g) + (I−K )· diag(σ) (20) where I is the identity matrix. 455 We let Λ(t) be the net emergence rate of mosquitoes into the patches from aquatic 456 habitats (see Eq. 18, above). The dynamics of adult mosquitoes are described by the 457 equation: 458 dM dt = Λ− Ω· M (21) 20 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Under the assumptions of this model, the density of gravid mosquitoes, G, is: 459 dG dt = f(M− G)− νG− Ω· G (22) This model thus assumes that only gravid mosquitoes can lay eggs (Eq. 13), but that all 460 mosquitoes (including gravid mosquitoes) can blood feed. 461 Other models for adult mosquito ecology, denoted dM/dt, could be formulated that 462 describe separate functions for mosquito survival and dispersal, depending on their 463 behavioral states (possibly including sugar feeding, mating and maturation), or that 464 describe a mosquito’s reproductive states, or its chronological age or reproductive age. 465 All models developed in this framework must accept the adult emergence rates, Λ, and 466 they must be formulated in enough detail to specify a population of egg-laying 467 mosquitoes, G, to compute ν (see Eq. 13). 468 Parasite Infection Dynamics in Mosquitoes The third core dynamical 469 component describes parasite infection dynamics in adult mosquito populations. Here, 470 we extend a previously published delay differential equation for the density of infectious 471 mosquitoes to include space and a time-varying extrinsic incubation period (EIP) [49]. 472 Let Y (t) denote the density of infected mosquitoes. Using κ from Eq. 10, the 473 dynamics of infection in mosquitoes are described by: 474 dY dt = f qκ(M− Y )− Ω· Y (23) We include a time-dependent EIP so that parasite development can be modulated by 475 temperature or other factors exogenous to the system: let τ(t) denote the EIP for a 476 mosquito that becomes infected at time t (i.e., it becomes infectious at time t + τ(t), 477 Supplement 4). We must also define the inverse τ−1(t), the delay for a mosquito that 478 became infectious at time, t. Let Υτ(t) denote a matrix describing survival and dispersal 479 of a cohort from time t− τ−1(t) through the EIP to become infectious at time t: 480 − ln Υτ(t) = ∫ t t−τ −1(t) Ω(s)ds. (24) When Ω and τ are constant, survival and dispersal through the EIP is Υ τ = e−Ωτ. 481 Otherwise, let the τ-subscript denote the value of a variable or parameter at time 482 t− τ−1(t). 483 To model the density of infectious mosquitoes, let Z(t) denote the density of 484 infectious mosquitoes. The dynamics of infectious mosquitoes are: 485 dZ dt = Υτ· fτ qτ κτ(Mτ− Yτ)− Ω· Z (25) The number of human blood meals per patch, called the net infectious biting rate, is 486 f qZ. 487 Models for infection dynamics, generically denoted dY/dt are nested within the 488 model for adult mosquito population dynamics dM/dt (for example, see [50]). These 489 models accept the net infectiousness (κ), they must define a variable describing the 490 21 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint density of infectious, blood feeding mosquitoes, Z, in order to compute the EIR (see 491 Eq. 9). 492 Epidemiology The fourth core dynamical component describes parasite infection 493 dynamics in human populations. Models for malaria infection, immunity, disease, and 494 infectiousness in humans, denoted dX /dt, can become quite complicated, depending on 495 the needs of a study. Studies of malaria epidemiology could consider the complex time 496 course of infections, superinfection, disease, detection, infectiousness, and immunity. 497 The state space describing malaria infection and immunity X can be modified to suit 498 the needs of a study, and the framework also has enormous flexibility to model 499 heterogeneity in populations through stratification. The following is one model family 500 that is complex enough to illustrate the generic features of the framework. 501 Let h = fh(E) denote the local daily force of infection (FoI) and δ(t) the FoI during 502 travel. In general, fh(E) could be modified to include heterogeneous biting [41], but in 503 this model, we assume h = bE. Both terms are defined for each sub-population. In 504 these models, we stratify on variables relevant for the epidemiology, including immunity, 505 and we model the effects by assigning different parameter values to each stratum. 506 To model infection dynamics, we modify a hybrid model for the multiplicity of 507 infection (MoI). The dynamics are based on a queuing model, in which new infections 508 occur at the rate h, and each parasite clears at the rate r, where we track apparent and 509 actual clearance as linked but distinct processes. The variables m1 and m2 track the 510 mean MoI for present and detectable parasites in each strata, which fully describe the 511 epidemiological state space in a simple model with superinfection [51]. We assume that 512 parasites clear at the per-capita rate, r1, so that: 513 dm1 dt = h + δ− r1m1 (26) In this model, the true prevalence is: 514 x1 = 1− e−m1 (27) We also formulate a model for the MoI of apparent infections. We assume parasite 515 infections are detectable for a shorter time so they appear to clear at a higher rate, r2, 516 and 517 dm2 dt = h + δ− r2m2 (28) Similarly, we let x2 denote the apparent prevalence 518 x2 = 1− e−m2 (29) We assume that if the infection is patent, a bite infects a mosquito with a higher 519 probability, c2, and c1 if it is not. A bite on a person in each stratum infects a mosquito 520 with probability: 521 x = c2x2 + c1(x1− x2) (30) To compute κ, the infective density of infectious resident hosts by strata is X = xH. 522 The vector X is passed to Eq. 10 to compute a vector of patch-specific net 523 22 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint infectiousness, κ. 524 To compute some of the spatial transmission metrics, including basic reproductive 525 numbers (see below), a model must compute the human transmitting capacity 526 (HTC) [52]. In this model, the number of days infecting mosquitoes at the higher 527 probability, c2 is 1/r2. The remaining days, are spent infecting mosquitoes at the lower 528 probability. Expressed as the equivalent number of perfectly infectious days, the HTC is: 529 D = c2 r2 + c1 ( 1 r1 − 1 r2 ) (31) This framework can accommodate other systems of equations describing parasite 530 infection and immune dynamics in humans. This particular model was designed to 531 illustrate some basic features of the modular design. These particular equations were 532 designed to incorporate the effects of immunity on transmission through stratification, 533 allowing parameters describing the duration of infections or detection and the 534 infectiousness to vary among strata ( e.g., r1, r2, c1 and c2). New models for human 535 epidemiology can use any epidemiological state space, X , and any system of equations, 536 dX /dt, including models with dynamical changes in the host population size. While the 537 travel FoI is recommended, it is not required. The modules should accept the EIR, and 538 to interact with other components, they must provide a function to compute the 539 infective density of infectious hosts, X. 540 Spatial Metrics 541 The Ross-Macdonald model defined a set of concepts and metrics that have formed a 542 basis for measuring and understanding malaria transmission, including vectorial capacity 543 and the basic reproductive number R0, but that model and associated metrics did not 544 include metrics for spatial dynamics, parasite dispersal, or malaria importation [3]. 545 Here, we define parasite dispersal by the set of locations ( i.e. patches) where 546 infecting bites occurred in continuous chains of transmission stretching back in time. 547 Dispersal for any parasite transmission chain is thus defined by locations of the bites 548 that caused each infection, and dispersal alternating between moving humans and 549 mosquitoes between bites. We acknowledge that, due to an observational process, there 550 is an important difference between where an infection occurred and where an infectious 551 person or mosquito is found. There is also an important difference between the formulas 552 defining dispersal and those used to compute reproductive numbers, which count from 553 after a host becomes infectious. Using this definition of parasite dispersal in the context 554 of a model, we have developed formulas and metrics to compute and study parasite 555 dispersal and reproductive success. 556 To develop these metrics, we assume steady state conditions. This is done for 557 convenience to avoid discussing the complications of understanding spatial dispersal 558 under dynamically changing conditions, and it is a necessary first step to understanding 559 such models. Analysis of malaria transmission dynamics under temporally varying 560 conditions are being developed in a subsequent manuscript. 561 The formulation of this static model helps to clarify the role of some of the 562 intermediate terms – if all parameters in a model were constant, the transmission model 563 23 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint could be fully defined by a much smaller set of parameters, but it may not be clear why 564 the parameters take on those values. Some of the terms that appear in the static 565 analysis correspond to parameters or variables in some Ross-Macdonald models, while 566 others are new: net emergence rates (Λ) or adult mosquito density (M ), scaled to the 567 appropriate human population density denominator of host availability ( W ), mosquito 568 bionomics (f, q, and Ω), and epidemiological parameters ( r1, r2, c1, and c2). New terms 569 describe the spatial biting distribution matrix ( β) and parameters describing malaria 570 importation (δ, υ, and xδ). 571 In models where the context is changing dynamically – due possibly to weather, land 572 use changes, or vector control – exogenous forcing functions can be passed to the model 573 that change resource availability or that perturb the dynamics; the functional forms and 574 intermediate terms (e.g. availability) are used to describe changes in the local 575 parameter values and guarantee mathematical consistency. In these static models, the 576 functions and terms are used to set up the model, but after setting parameter values, 577 they need not be called again. 578 Fig 6. To model malaria importation, we define a travel FoI for each stratum, δ(t), and two set of terms to model the role of visitors in mosquito blood feeding and parasite transmission: the available visitor population Wδ and the NI for the visitor population, by patch xδ. To model blood feeding and transmission, we compute a patch-specific resident fraction for blood feeding, υ, the fraction of all biting that occurs on a resident of the spatial domain. From this, we can compute the visitor reservoir fraction, γ, the travel fraction for incidence, and other measures of malaria importation. Net Malaria Importation and Travel Fractions 579 Terms describing the travel FoI (δ) and visitor populations were defined above and 580 integrated into the models for blood feeding and human epidemiology. We define an 581 imported malaria case as a human infection that traces back to a location outside of the 582 spatial domain in the parasite’s previous generation, i.e., the mosquito and human host 583 preceding this one in a chain of infections [42]. Net malaria importation rates describe 584 the number of imported malaria cases, per day. 585 The fraction of all cases that were imported called the travel fraction can be defined 586 as either: 1) the fraction of incident infections that were imported; or 2) the fraction of 587 prevalent infections that were imported [45,53]. To compute these travel fractions, we 588 24 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint let γ = (1− υ)xδ/κ denote the visitor fraction, the fraction of infectious mosquitoes that 589 were infected by visitors. We let h denote the FoI. The travel fraction for incidence is: 590 hγ + δ h + δ (32) The travel fraction for true prevalence is: 591 1− e−(δ+hγ)/r1 1− e−(δ+h)/r1 (33) We note that these are per-capita terms defined for the strata. The net malaria 592 importation rate, the number of imported malaria incidence per day for each patch is: 593 J· (hγ + δ) H (34) so the travel fraction for incidence for the patches would be: 594 J· (hγ + δ) H J· (h + δ) H (35) Formulas for the travel fraction for prevalence are formulated in the same way. 595 Parasite Dispersal 596 To compute quantities related to parasite dispersal, from bite to bite, we focus on local 597 transmission, and we need some formulas that describe how mosquitoes move around in 598 humans and in mosquitoes. 599 Mosquito Dispersal and Steady States In these models, we can compute steady 600 state mosquito population density, assuming Λ is constant over time. At the steady 601 state of Eq. 21, 602 M = Ω−1· Λ (36) Here, the inverse Ω−1 can be understood as a measure of time spent alive in each patch 603 by mosquitoes emerging habitats in each patch. In other Markov chain models with 604 finite state space, it has also been shown that the elements of the matrix inverse can be 605 interpreted as residence times [54,55]. In the simpler Ross-Macdonald model, the 606 inverse of a mortality rate, g, is a measure of time spent alive or the average mosquito 607 lifespan [56,57]. The time spent alive interpretation of Ω−1 is more apparent if there is 608 no movement: if we set σ = 0, then Ω−1 = diag(1/g). 609 In spatial models, the matrix Ω accounts for both survival and movement. To 610 illustrate – and to demonstrate that if Ω is a sensible description of mosquito 611 demography, then the matrix inverse must exist – we construct a tracking matrix. Let 612 Ξ(t) denote a matrix that tracks cohorts of mosquitoes: 613 Ξ(t, M0) = e−Ωt· diag (M0) (37) It describes the density of mosquitoes left from an initial cohort in each patch M0 that 614 is found in each location at each point in time. There is a duality between the 615 25 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint equilibrium population density from Eq. 21 and time spent alive by a cohort, computed 616 by integrating Eq. 37 (i.e. orbits of the related equation dM/dt =−Ω· M). Just as we 617 can compute g−1 = ∫∞ 0 e−gtdt, we can compute: 618 M = Ω−1· Λ = ∫ ∞ 0 e−Ωtdt· Λ (38) so that the steady state can be found by simply adding up the time spent alive in each 619 patch by a cohort emerging from every other patch. Under generalized static conditions 620 (i.e. σ > 0), Ω−1 can thus be interpreted as the average time spent alive in every patch 621 by cohorts of mosquitoes initially found in each patch. 622 Parasite Dispersal in Mosquitoes Using mosquito tracking matrices, we can also 623 track parasite dispersal in mosquitoes to derive a matrix that has the same 624 interpretation as the formula for vectorial capacity [57,58]. 625 To transmit, mosquitoes must blood feed on a human to become infected: the net 626 infection rate in each patch, per available human, is f qκM/W. After becoming infected, 627 a mosquito must survive while dispersing through the EIP (Υ = e−Ωτ). After becoming 628 infectious, a mosquito must blood feed to transmit parasites, so we use the matrix 629 inverse Ω−1 which describes where the mosquitoes are for each infectious human blood 630 meal as long as they remain alive; after becoming infectious, the distribution of 631 infectious bites is given by f qΩ−1. We can describe parasite transmission by mosquitoes 632 by following the story of infection in mosquitoes: after emerging ( diag(Λ)), a mosquito 633 must blood feed on a human to become infected ( f qΩ−1/W ); then survive the EIP 634 (e−Ωτ); and then blood feed to transmit ( f qΩ−1). 635 In the Ross-Macdonald model, the formula for vectorial capacity can be derived from 636 the formula for the daily EIR as a limit [57]. In spatial models, a vectorial capacity 637 matrix can be derived as the limit of a tracking matrix describing the number of 638 infectious bites arising, per available person (i.e., the denominator is W ), per day at the 639 steady state (Supplement 4): 640 V = lim κ→0 f qZ W = f qΩ−1· e−Ωτ· diag (f qM W ) (39) Elements in the matrix V are the expected number of infectious bites eventually arising 641 in every patch from all the mosquitoes in a single patch blood feeding on a single human 642 on a single day, computed as if each human were perfectly infectious. The derivation 643 assumes that no mosquitoes are already infected, and the assumption that humans are 644 perfectly infectious ( κ = 1) is made so that the formula deals only with phenomena 645 related to mosquitoes. In models with multiple vector species, the notion of what it 646 means to be “perfectly infectious” is not as simple because of differences among vector 647 species in their capacity to be a host for the parasites, or vector competence 648 (Supplement 4). 649 Parasite Dispersal by Humans To quantify parasite dispersal by humans, we 650 compute the human transmitting capacity distribution (HTCD) matrix. We let human 651 transmitting capacity (HTC) describe the net number of perfectly infectious days for 652 26 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint each stratum: since infectiousness varies over the time-course of infections, we sum 653 partially infectious days and interpret the HTC as an equivalent number of days spent 654 perfectly infectious [52]. For the population strata in this model, the HTC ( D) is 655 defined by Eq. 31. Since transmission requires two bites, we use the TaR matrix to 656 determine both where a human becomes infected and where it infects a mosquito. Using 657 the transposed TaR matrix, we can describe where infectious days at risk are spent, 658 ΨT· D. Parasite dispersion by mosquitoes for the sub-populations also accounts for 659 where a mosquito becomes infected, or bΨ. 660 The HTCD matrix uses the biting distribution matrix, β, to count from the 661 infectious bite and weight biting appropriately for subsequent blood feeding by all the 662 population strata. The HTCD, a p× p matrix (D), is: 663 D = diag (W)· βT· diag (bDH)· β. (40) We note thatD in spatial models is analogous to bD in models with a single patch. 664 (The equivalency ofD and bD is most apparent if no humans move, and if there is one 665 stratum per patch, and if all search weights are 1, in which case H = W and 666 β = diag(1/H).) Like bD,D describes days spent infectious by an individual human, 667 but inD, describes both where a human got infected and where the mosquitoes were 668 subsequently infected. 669 The definition ofD as a time-dependent matrix is substantially more complicated if 670 local human mobility patterns change dynamically. 671 Parasite Dispersal through one Parasite Generation Parasite dispersal is 672 defined by the locations where infecting bites occurred, alternatively moving in infected 673 mosquitoes and humans. The equations for D andV describe the expected movement 674 for a parasite among patches in humans or mosquitoes, respectively, counting from bite 675 to bite. Notably, the formulas are defined for a parasite in either a mosquito or a 676 human. We can also define parasite dispersal through one parasite generation ( i.e., from 677 human to human, or from mosquito to mosquito) but the formula depends on where we 678 start counting. If we started from all the mosquitoes blood feeding on a single human 679 (averaged appropriately) on a single day in every patch, then we would get a matrix 680 describing dispersal from every patch to every patch: 681 D·V . (41) If we started counting from a typical human infected in a patch on a single day, we 682 would get a different dispersal matrix: 683 V·D . (42) Importantly, these formulas follow the same process in the same order, and thus closely 684 resemble the reproductive numbers for malaria (described below), which measure 685 reproductive success for a single parasite. These formulas are two among many that 686 could be developed to count events through a parasite’s life-cycle starting at different 687 points. 688 Formulas that describe the parasite’s per-capita reproductive success, such as 689 27 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Eqs. 41-42, counting events arising from a single host. In some cases, we might wish to 690 count the total number of events arising from a patch. To measure the contribution of a 691 patch to overall transmission, we must have a measure of connectivity, or total parasite 692 flows. A tracking matrix describing all of the infections arising from each patch on a 693 day, is: 694 diag (W)·D·V (43) If we started counting infections occurring on humans in a patch, we would get an 695 alternative patch-based tracking matrix. The number of infections arising from a patch 696 is thus tracked by: 697 diag (W)·V·D (44) These measures emphasize the role of places with larger available populations. 698 The same sort of formulas can be devised to describe transmission from human 699 strata to human strata, but the resulting formulas are only spatial insofar as the human 700 strata are anchored to a residency. If we focused instead on parasite reproductive 701 success starting with an infection in humans, regardless of location, we would get 702 R = bβ·V· diag (W)· βT· diag (DH) . (45) or we could also count bulk transmission from humans as diag(H)·R . Notably, Eq. 45 703 is a stratum-based measure. To make it quasi-spatial, we would need to assign events to 704 patches by stratum residency using the membership mapping operator J·R·J T . 705 Distances Dispersed To get a measure of the distribution of distances travelled by 706 parasites, we match a measure of transmission intensity with the corresponding element 707 in a patch distance matrix describing the distance. We take the couplet (distance and 708 intensity) and sort by distance, then compute the cumulative distribution function 709 (CDF). From the CDF, we derive a probability mass function [39]. These dispersal 710 kernels provide a simple way of visualizing distances dispersed by mosquitoes, humans, 711 or parasites. 712 These formulas and algorithms draw attention to the differences in metrics 713 describing parasite transmission dynamics and dispersal. Because of spatial 714 heterogeneity in mosquito and human population densities, there are many sensible 715 formulas for counting dispersal, some of which correspond to describing rates, ratios, 716 proportions, and numbers. Careful thought should be given to choosing or developing a 717 metric that fits the analysis. 718 Reproductive Numbers 719 Reproductive numbers are a measure of the parasite’s average reproductive success. 720 When transmission is spatially heterogeneous, reproductive success will vary for 721 parasites, depending on where they are. As parasites spread over several generations, 722 the expected success of its progeny will change. To calculate threshold criteria for 723 persistence (in the absence of malaria importation), we want a reproductive number to 724 be a measure of average success taken over the whole system, but we want to use an 725 average that does not change across generations. Doing so requires that we compute the 726 28 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint spectral average, which is computed as the dominant eigenvalue of the parasite’s next 727 generation matrix. 728 For many reasons, it is useful to formulate local reproductive numbers that describe 729 a parasite’s average reproductive success at a particular place and time – an arithmetic 730 average. These local reproductive numbers could ignore differences across generations, 731 so they would not serve as thresholds for parasite persistence. In this section, we define 732 local reproductive numbers at the steady state, but the formulas could also serve as 733 point estimates. 734 Reproductive numbers describe malaria transmission under a range of different 735 conditions that are relevant for understanding malaria transmission dynamics and 736 control or for national strategic planning. Baseline conditions are described by the basic 737 reproductive number, R0, which is defined for a population with no acquired immunity 738 and no malaria control. The adjusted reproductive number, RC, describes a family of 739 numbers defined for a population with no acquired immunity adjusted by malaria 740 control, at a fixed level of control denoted C. In other words, R0 is defined as a special 741 case of RC, but in the absence of control. The total effect size of malaria control on 742 transmission is R0/RC. Here, we also describe the endemic reproductive number, RE, 743 which describes potential transmission modified by immunity. The total effect size of 744 immunity on transmission is RC/RE. In computing RE, as with R0 and RC, we ignore 745 the fact that some hosts are already infected. In this way, RE is defined differently than 746 the effective reproductive number, denoted Re, which is lower than RE because it does 747 not count infections occurring in someone who is already infected. We note that, by 748 definition, at an endemic steady state Re = 1. By way of contrast, RE counts the 749 number of infections that would occur after one generation, which is useful for planning 750 because it helps to clarify how success in malaria control can be assisted by immunity 751 that will eventually wane. 752 Both R0 and RC are computed as if there were no acquired immunity. In this model, 753 the effects of acquired immunity on transmission are quantified through the stratified 754 values of b, r1, r2, c1 and c2. These parameters determine the HTC for all the strata ( D, 755 see Eq. 31). If D were computed using values that have been tuned to a stratum with 756 some level of immunity, we would be computing RE. To compute RC, we would need to 757 replace D with values set to a non-immune baseline ( i.e., D0), and then recompute the 758 next-generation matrix. Next generation matrices computed with values of D that 759 include the effects of acquired immunity are thus describing an endemic reproductive 760 number. Depending on how D is computed, and whether the bionomic parameters 761 incorporate effects of vector control, we may thus be computing R0, RC or RE. 762 Local Reproductive Numbers One way to define local reproductive numbers is to 763 modify Macdonald’s formula using the local values of parameters, as if there was no 764 movement of mosquitoes or humans. To write the formula using some models in this 765 framework, we may need to modify HTC (which is defined for the strata, of length n) to 766 take a patch average. To compute a patch average HTC, ˘D (a vector of length p), we 767 take the population weighted average, 768 ˘D = ΨT· wf DH ΨT· wf H (46) 29 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint We can then describe a local reproductive number, ˘RC (or possibly ˘RE, depending on 769 how the interpret parameters are defined in D): 770 ˘RC = Λ W f2q2 g2 e−gτ ˘D (47) This local measure is similar to Macdonald’s formula [59]. While useful in some 771 contexts, the formula should be applied with caution. 772 An alternative way to compute local reproductive numbers uses V andD (perhaps 773 modified to remove the effects of immunity on transmission). Since the matrices count 774 infections arising from each patch, and we add all infections arising to the patch where 775 the bite originates. We let 1 be a row vector of ones of length p, and we can count 776 infections arising starting from all the humans infected in a patch on a single day: 777 ˆRC = 1·V·D (48) that counts infections occurring on humans, or we can start from all the mosquitoes 778 blood feeding on humans on a single day, and: 779 ˜RC = 1·D·V (49) that counts infected mosquitoes. These patch reproductive numbers could provide 780 valuable information about whether to target the mosquitoes or humans in some patch 781 for enhanced interventions. We could also consider the equivalent formulas for total 782 patch outputs: 783 WT·V·D or WT·D·V (50) where WT is a row vector. Alternatively, we can also weigh transmission from strata 784 using Eq. 45: 785 1·R (51) or the equivalent scaled by stratum size: 786 HT·R, (52) where HT is a row vector, which gives us valuable information about infections arising 787 from every stratum on every strata, a way of identifying the relative importance of 788 various population strata. 789 Next Generation Matrix In the Ross-Macdonald model, a parasite’s reproductive 790 success in the next generation is described by a single number. It is computed by 791 counting forward from the moment a mosquito or human becomes infectious. Since 792 parasites move in infected mosquitoes and humans, parasite reproductive success – 793 measured as the number of infections in the next generation – varies across generations 794 as the parasite distributions evolve across generations among strata and among patches. 795 The matricesV andD describe parasite transmission and dispersal in mosquitoes and 796 humans, respectively. While the product of these formulas does describe net 797 reproductive success, the computation of threshold conditions has been developed 798 30 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint around the concept of a next generation matrix [60,61], which traces the same process 799 in the same sequence but that start counting at a different point in the parasite’s life 800 cycle (Fig. 7). A threshold condition is found by taking the spectral average of the next 801 generation matrix. 802 Fig 7. A Spatial Life-Cycle Model . A diagram that illustrates how the parameters describing each stage in the parasite’s life-cycle translate into a parasite’s reproductive success spatially, when mosquitoes and hosts move. The right half of the circle represents mosquitoes and the left half humans. The flow of events is clockwise. Mosquitoes must blood feed to become infected ( f qM), and then survive and disperse through the EIP ( e−Ωτ). infectious bites are distributed as long as a mosquito survives, while it blood feeds and disperses ( f qΩ−1). The bites are distributed among humans ( β) and some of them cause an infection (b). Parasites are transmitted for as long as humans remain infectious, measured in terms of the human transmitting capacity (HTC, or D days). Infectious humans are distributed wherever humans spend time at risk (affecting β). These processes are summarized differently to model parasite dispersal and parasite reproductive success. Dispersal counts from bite to bite using the VC matrix (V) and the HTC matrix (D). Threshold computations count from when a host becomes infectious to measure a parasite’s reproductive success in infectious mosquitoes (RZ); in infectious humans (RX); from human to humans among strata after a human becomes infectious (R); and from mosquito to mosquitoes (Z). R0 is the lead eigenvalue ofR orZ. Under endemic conditions, we can also consider how frequently parasites are actually transmitted by including the probability a mosquito gets infected κ, and the probability a mosquito is infectious, given by the sporozoite rate z. In computing next generation matrices, we focus on transmission within a defined 803 spatial domain. For mathematical convenience here, we thus set υ = 1, though we could 804 easily develop matrices leaving υ undetermined to discount exported malaria cases. 805 We first compute offspring transmitted from a single infectious mosquito to humans 806 or from a single infectious human to mosquitoes, each of which defines a stage in the 807 parasite’s next-generation [60]. After a mosquito has become infectious, how many 808 humans (in each stratum) would it infect? In these models, the answer to that question 809 is n× p matrix, denoted RZ, describing transmission from an infectious mosquito in 810 each patch to humans in each strata: 811 RZ = bβ· f qΩ−1. (53) 31 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint How many infectious mosquitoes would arise from each human infection? The answer is 812 a p× n matrix, denoted RX, describing transmission from a human in each stratum to 813 mosquitoes: 814 RX = e−Ωτ· f qM· ( βT· diag (DH) ) . (54) The next-generation matrix by type is: 815 G = [ 0 RZ RX 0 ] (55) To describe reproductive success in terms of the parasite biology, we count reproductive 816 success through one full parasite generation, either from humans back to humans, or 817 mosquitoes back to mosquitoes. For the parasites, reproductive success through one full 818 generation requires two events, one of each type, so we square the matrix given by 819 Eq. 55 to get a new matrix in block form: 820 G2 = [ R 0 0 Z ] . (56) We thus get two diagonal block sub-matrices describing reproductive success in the 821 parasite’s next generation, denoted R andZ. Reproductive success from human 822 population strata back to human strata is described by an n× n matrixR = RZ· RX: 823 R = bβ·V· diag (W)· βT· diag (DH) . (57) Reproductive success from mosquito through the population strata back to mosquitoes, 824 described patch-by-patch is described by the p× p matrixZ = RX· RZ: 825 Z = e−Ωτ· diag (f qM W ) ·D· f qΩ−1 (58) We have also formulated the next-generation matrix for systems with multiple vector 826 species (Supplement 4). 827 The Spectral Average We can also compute RC as a spectral average through 828 simulation, which is one useful way of illustrating what a spectral average means. To do 829 so, we choose a vector describing the distribution of parasites in a founding generation, 830 X0 orY0, and iterate parasite infections across i successive parasite generations: 831 Yi+1 =ZYi or Xi+1 =RXi. (59) We define the vector: Ei =Xi+1 ∥Xi∥ or Ei =Yi+1 ∥Yi∥ . where∥X∥ or∥Y∥ is a scalar that denotes is magnitude. Over many generations, Ei 832 converges to the lead eigenvector, a scalar value also called the spectral average or RC: 833 RC = lim i→∞ ∥Ei∥ (60) 32 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint and it is interpreted as the asymptotic average reproductive success expressed as a 834 number of infected hosts per host, per generation. Note that it is asymptotic only for 835 the linearized system defined by Eq. 55 or Eq. 56. 836 Quasi-Thresholds for Endemic Malaria Without malaria importation, RC > 1 is 837 a threshold criterion. Analysis of models without malaria importation have consistently 838 demonstrated that malaria is either absent or that there is a single globally, 839 asymptotically stable equilibrium. When there is imported malaria, there are three 840 sufficient criteria for some local parasite transmission to occur within the area: 841 1. max{δ} > 0 and RC > 0; 842 2. max{(1− υ)Xδ} > 0 and RC > 0; 843 3. RC > 1. 844 If condition 1 or condition 2 is satisfied, then malaria will be present in an area, and if 845 RC > 0 then there will be some local transmission. If RC > 1, malaria transmission 846 would be sustained in the absence of importation. We thus call RC > 1 a 847 quasi-threshold for endemic transmission to occur within the spatial domain: endemic 848 describes places where RC > 1, and pseudo-endemic places where 0 < RC < 1 with 849 significant levels of transmission. 850 Quantifying Transmission in a Place 851 The framework, models developed within it, and the associated spatial metrics were 852 designed to have the skill required to describe and quantify heterogeneous spatial 853 transmission dynamics of malaria in a specific place at a particular time. We have not 854 explicitly defined algorithms for the observational processes that would map model 855 states onto observable quantities, which would be required to extend this mathematical 856 modeling framework into a state space modeling framework to rigorously fit models to 857 data. Instead, we have focused on the mathematics of these processes: time spent by 858 humans; other blood hosts; daily mosquito rhythms; mosquito host preferences, time at 859 risk; and mosquito mobility. Similarly, the models for mosquito ecology and population 860 dynamics describe the mathematics of mosquito mobility, in terms of explicit 861 assumptions about the locations of aquatic habitats, heterogeneous distributions of 862 resources, and mosquito mobility patterns that emerge from a search for resources. By 863 quantifying spatial patterns in terms of the underlying processes – including malaria 864 importation, mosquito ecology and spatial population dynamics, parasite transmission 865 dynamics, human mobility, and malaria epidemiology – the equations point towards a 866 general inferential framework. 867 Models developed within this framework involve substantially more parameters than 868 the Ross-Macdonald model. This is an inevitable consequence of a decision to model 869 transmission at a particular place and time. If any local features are important for 870 transmission, then a larger set of quantities must be estimated to understand and 871 quantify those features. This gives rise to an important but difficult practical question: 872 What is the relationship between the amount of local intelligence and the specificity of 873 33 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint the policy advice that can be offered? With minimal local information, it is possible to 874 offer generic policy advice, but it may not be necessary to know everything about a 875 place to tailor advice to context. With this framework, it is possible for models to 876 evolve as the amount information increases, and the models may be used to look ahead 877 to prioritize missing data: How can programs identify missing information that would 878 most rapidly improve the effectiveness of malaria control? These contextual factors and 879 the related questions are addressed below. 880 Malaria Landscapes While the Ross-Macdonald model describes parasite 881 transmission between abstractly defined mosquito and human populations, the 882 framework we have described was developed to understand and quantify malaria 883 importation and transmission among structured mosquito and human populations in a 884 well-defined geographical area. (Using a model dX /dt that describe the infection 885 dynamics of other pathogens and immunity in vertebrate host populations, and making 886 other appropriate choices, the framework could be used as a basis for modeling dengue, 887 West Nile virus, or other mosquito-borne pathogen transmission dynamics, as well.) 888 Since the models are developed to approximate malaria transmission in an actual place, 889 after defining an observational process, the model outputs would be verifiable 890 statements about real quantities over some specific period of time. 891 As a practical first step, model building starts by defining a set of structural 892 elements – patches, human population strata, and aquatic habits – that are appropriate 893 for the needs of a study (e.g. Fig 8 illustrates some options for simulating malaria on 894 Bioko Island, Equatorial Guinea). A geographical study area is usually defined by 895 projects, programs, or political boundaries. In planning interventions for a defined area, 896 an important concern is connectivity to surrounding areas. How much malaria is 897 imported by daily human movement or travel? Are the mosquito populations within the 898 area strongly connected to others nearby? 899 Using spatial metrics to identify differences in transmission patterns and the flow of 900 parasites across a landscape can help control programs prioritize drugs, outreach, and 901 medical attention to populations, and vector or larval control to places. Using our 902 differential equation framework to reconstruct the equilibrium analysis presented in [45], 903 we have generated spatial bulk transmission matrices (diag(H )·R) among areas for 904 Bioko Island, Equatorial Guinea. In Fig. 9 different patterns of pathogen transport are 905 readily apparent between persons who live in Malabo (left), the densely populated 906 capitol of the island and a sink for travellers, and Luba (right), a small settlement in the 907 Southern half of the island. The pattern of travel seen in Luba typifies most of the areas 908 outside of Malabo, where individuals most often travel to the capitol but not to the 909 other outlying settlements. These patterns affect transmission, where we see parasites 910 originating in Malabo tend to stay in the city. Parasites originating in Luba either tend 911 to stay highly local, or are transported to Malabo when those persons move. Because 912 malarial mosquitoes tend to fare less well in urban settings, these spatial metrics can 913 help understand how high prevalence can be sustained in otherwise unsuitable locations. 914 An equally important question is about heterogeneity in mosquito population 915 densities within the area and heterogeneity in the risk of exposure, which should inform 916 the definition of patches and the choice of a patch size. Patches, in this model, are 917 defined around adult mosquito activities, and each “patch” has a geographical location. 918 34 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint The patch is the spatial unit that defines the algorithms for time spent, blood feeding, 919 egg laying, adult mosquito survival and dispersal. The concept of a patch is flexible 920 enough to model blood feeding indoors and outdoors at the same geographical locations, 921 which may be useful to inform programmatic questions about the effectiveness of vector 922 control measures that target indoor biting (Fig. 10). Since the patch is the basis for 923 computing most aspects of blood feeding, the patches define the structure for human 924 time spent at risk, including (if required) quantifying time spent indoors vs. outdoors, 925 and mosquito movement rates from indoors to outdoors, from outdoors to indoors, or 926 from outdoors to other outdoor patches. 927 An important basic concern is the spatial granularity of the patches used for 928 simulation (see Fig. 8). Some questions remain unresolved about the appropriate spatial 929 scales and ways to define patches for describing and analyzing malaria transmission for 930 policy (e.g., to compute IRS coverage). One advantage of this framework is that it is 931 possible to build nested models with different spatial grains and compare them. Smaller 932 patches more accurately capture heterogeneity in a landscape while increasing the 933 number of parameters that need to be inferred during calibration to data.k 934 Aquatic habitats are located in patches, but the model was designed to assign 935 patches to habitats assuming the habitats had an actual location. Patches in this 936 framework need not have any human residents or any available hosts, so that mosquito 937 habitats in the uninhabited areas around human households are contributing to 938 transmission. Mosquito population dynamics are coupled through related equations 939 describing gravidity, egg laying and egg deposition. The framework thus does not 940 impose any constraints on either the method for constructing patches, or on the number 941 or arrangement of aquatic habitats within the spatial domain. Given the modular 942 nature of these models, the dynamics of immature mosquito populations in each aquatic 943 habitat depend only on its parameters and the egg deposition rates. The productivity of 944 any one aquatic habitat in an area is, however, coupled to other habitats through egg 945 laying by adult mosquitoes that could have emerged anywhere. 946 To improve the accuracy of models, human populations can be segmented into strata 947 to reduce heterogeneity in traits that affect malaria: the first segmentation is by 948 residency. In this framework, which is designed to quantify process affecting 949 transmission, heterogeneity in any trait affecting transmission is dealt with by 950 sub-dividing the population into homogeneous (or less heterogeneous) strata, such as by 951 age, travel habits or patterns, ITN usage, vaccination, care seeking, or any effects of 952 immunity affecting malaria epidemiology or transmission. 953 Notably, all this structural flexibility is achieved through membership matrices and 954 through the variables describing resource availability, which links search weights, 955 functional responses, and other functional forms to guarantee mathematical consistency 956 (e.g. avoiding problems when denominators are zero) despite structural changes. Suites 957 of models can be developed to address concerns about data gaps and uncertainty that 958 are appropriate for studies. Model complexity can be modified by changing dynamical 959 modules, by changing functional forms, by fixing or changing parameters, by splitting 960 and joining patches, by splitting or joining strata, or by adding and subtracting aquatic 961 habitats. With the ability to split and join patches or strata, any model can be mapped 962 onto simpler, nested models in a series of simple join operations until it is collapsed onto 963 a single-patch, single-stratum Ross-Macdonald model. This is functionally what is 964 35 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint meant by scalable complexity. 965 It is thus as easy to modify and evaluate the effects of model structure ( e.g. the 966 number of strata) as it is to vary parameters, to facilitate developing suites of models, 967 including models with nested patches or nested strata, to explore tradeoffs in building 968 and calibrating models at various levels of detail. 969 Mosquito Blood F eeding and Ecology Three constant parameters describing 970 mosquito behavior are a standard part of the Ross-Macdonald model [56,57]: the daily 971 death rate of mosquitoes ( g), the overall daily blood feeding rate ( f), and the human 972 blood feeding fraction ( q). Incorporating the possibility of dynamical feedback between 973 the future emergence of adults and current population size means we have added the 974 population egg-laying rate (Γ). Adding spatial complexity to the model means the daily 975 emigration rate ( σ), mosquito dispersal (K ), distribution of habitats (N ) and the 976 distribution of eggs among patches (U) are additional parameters which define how 977 populations may interact in space. While our analysis has focused on steady states, the 978 models were formulated with parameters that can vary over time in response to 979 changing availability of resources [33–36]. 980 In this framework, the values of all these parameters are computed with functional 981 responses based on resource availability, mosquito biology and innate preferences that 982 constrain the parameters within sensible ranges. This formulation emphasizes how 983 baseline mosquito bionomics for different species could respond to available resources 984 and how those responses would be modified by control. In particular, the same human 985 behaviors can give rise to very different blood feeding patterns for different vector 986 species, depending on the daily rhythms, host preferences, and aquatic ecology of 987 different vector species (Supplement 4). We thus have a basis for understanding 988 mosquito behaviors and ecology as a baseline that may have been modified by vector 989 control or weather. 990 Blood feeding in this model thus makes an important distinction between 991 anthropophily, or innate mosquito preferences for hosts of different types, and 992 anthropophagy, summarized by the human blood feeding rate ( f q). Models can also 993 consider a difference between the time of day when mosquitoes are actively searching for 994 blood (ξ ) and the blood feeding rates by time of day (f ), which vary with host 995 availability. Innate, species-specific host preferences are embodied in functional forms 996 and parameters, while the rates describing what has happened also depend on context. 997 Similarly, mosquito population dynamics are an emergent feature of a resource 998 landscape. Since searching for resources is also associated with resource availability, 999 adult mosquitoes will tend to aggregate in patches that have habitats and other 1000 required resources. In these models, egg-deposition rates in habitats by volant adult 1001 populations are spatially heterogeneous and only partially determined by the emergence 1002 rates of adults from a single habitat. The concept of a carrying capacity is, perhaps, not 1003 as useful as the concept of habitat productivity and the functional forms that determine 1004 how the number of adults emerging is related to the number of eggs laid [31]. A 1005 habitat’s carrying capacity only makes sense in the abstract – if adult mosquitoes 1006 emerging from a single habitat only laid eggs in that natal habitat. In this framework, 1007 the aquatic population dynamic module determines how adult mosquito emergence 1008 rates respond to egg laying by the adult population. 1009 36 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint The parameters describing these processes are both habitat-specific and 1010 time-dependent: density-independent mortality, density-dependent mortality, the 1011 response to crowding, maturation rates, and search weights could vary for every habitat. 1012 A habitat can thus disappear seasonally (which occurs when wν = 0), or weather could 1013 affect immature mosquito maturation and mortality rates. If a study called for 1014 modeling resource-based competition or stage-structured mosquito populations, the 1015 equations describing aquatic populations (dL/dt) can be modified as needed (Fig 1). 1016 The framework thus facilitates the construction of realistic models of mosquito ecology, 1017 insofar as it is justified by data available and the needs of a study. 1018 Local Exposure, Human Biting Rates and Mixing In defining the algorithms 1019 for blood feeding, we also developed a new model for the human biting rate (HBR) and 1020 by extension, the entomological inoculation rate (EIR), two basic metrics used to 1021 measure malaria transmission entomologically. 1022 The model emphasizes that for any population stratum, the risk of exposure to 1023 biting mosquitoes is distributed spatially. In these models, this is determined by a biting 1024 distribution matrix (β ). A similar matrix has appeared in other models for the spatial 1025 dynamics of mosquito-borne diseases for which human mobility is based on a concept of 1026 ”visitation” or time spent – classified as Lagrangian movement [7,8,10,12 –15,17,18,45]. 1027 Here, β is based on a concept of availability, the weighted, ambient population at risk. 1028 Availability is computed from observable quantities, and it is computed dynamically for 1029 arbitrarily defined human strata and changing availability (the denominator). The 1030 formulas guarantee consistency in blood feeding: the number of human blood meals 1031 taken by mosquitoes is equal to the number of blood meals received by the humans. 1032 In the new model, the HBR is defined as β· f qM and the EIR is β· f qZ, so that the 1033 number of bites received by each stratum depends on how they spend their time at risk. 1034 In studies that have reported a value for the HBR or EIR, the quantity reported is 1035 based on catch counts by a person or device in a place. In this model, the quantity that 1036 is closest to the quantities being estimated is pHBR or pEIR, the number human blood 1037 meals, or infectious human blood meals in a patch, per available person, per day 1038 (f qM/W or f qZ/W). A person who is in a patch at a particular time of day would 1039 experience the local biting rates at that time scaled by a search weight ( f qM ξ(t)ωf /W 1040 or f qZξ(t)ωf /W ). The quantity being estimated by human landing catches is a 1041 measure of the intensity of exposure in a place. 1042 Since other hosts are also available, the number of mosquitoes caught also depends 1043 on the biases of the trapping method. In this model, each method for trapping 1044 mosquitoes can be thought of as having its own “availability,” and it is competing for 1045 the attention of mosquitoes. Each method for catching mosquitoes is biased in some 1046 unknown way. We thus suggest that field methods designed to estimate the EIR are 1047 best interpreted as a location-specific measure of risk in a place, and that 1048 epidemiologically relevant measures of risk must acknowledge exposure occurring for a 1049 period of time, including all the places where a person spends time. The pEIR, weighted 1050 by total availability, is a good approximation of the EIR only if a person spends most of 1051 their time at risk in that place. The formulas presented here are useful to quantify how 1052 local measures of mosquito blood feeding in a place could differ from what the humans 1053 living in that place would experience. What is the difference between risk for a human 1054 37 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint who moves around compared to their counterfactual self who never leaves home? 1055 The Spatial Scales of T ransmission Important considerations for planning, 1056 monitoring, and evaluating malaria control are the spatial scales that characterize 1057 transmission, as defined by parasite dispersal in mosquitoes or humans. We have 1058 defined parasite dispersal rigorously in terms of the locations where blood meals 1059 occurred that transmitted parasites in dispersal chains. While these definitions are 1060 compelling, the distribution of distances separating every pair of infectious bites in a 1061 chain of malaria infections can only be approximated using other data. In practice, the 1062 framework we have described makes a distinction between local transmission and 1063 imported or exported malaria. The framework makes the most sense mathematically if 1064 most transmission is local, but the framework also defines quantities for malaria 1065 importation and exportation, making it possible to study connectivity using a frame 1066 that shifts among spatial domains and across spatial scales. 1067 After drawing a bounding box to define a spatial domain and a set of patches, we 1068 classify any pair of bites in a transmission chain where at least one occurred in the 1069 patch: either both bites occurred somewhere in the spatial domain, called local 1070 transmission; or the first bite occurred outside the spatial domain, called imported 1071 malaria; or the second bite occurred outside the spatial domain, called exported malaria. 1072 These measures of imported and exported malaria thus provide a basis for 1073 understanding and quantifying dispersal within and among defined geographical areas. 1074 These models weigh the consequences of imported malaria, but as a practical matter, 1075 the importance of exported malaria is difficult to quantify because the expected number 1076 of subsequent bites depends on conditions somewhere else. Importantly, the fraction 1077 that stays local may differ depending on whether the parasite is moving in a mosquito 1078 or a human. Similar definitions and arguments would apply to transmission through a 1079 parasite, a full parasite generation encompassing three bites and two jumps. The 1080 metrics we have developed describe transmission within a defined geographical domain, 1081 but if there is a need, the models can be reformulated for a larger spatial domain. 1082 The models and metrics provide a way of characterizing the spatial scales of 1083 transmission by computing the cumulative fraction of all transmission occurring within 1084 a circle of a given radius. Sensible points on that curve can be compared by patch: 1085 What distances contain 80%, 90%, 95%, or 98% of all transmission? These estimates 1086 are, out of necessity, based on estimated quantities – models of mosquito mobility, 1087 human mobility, and modeled mosquito population density – about which there is 1088 substantial uncertainty. 1089 Despite the overall uncertainty, these spatial scales are constrained by limits on time 1090 and travel. Some quantities are known from census data ( e.g. population distributions). 1091 Most mosquito dispersal distances are short. Mosquitoes can move large distances, but 1092 most stay within 1 km of a natal habitat [62]. For humans, the fraction of time spent 1093 declines sharply with distance away from home. A large fraction of time is spent at 1094 home, especially at night, and a larger fraction of the time is spent within roughly 10 1095 km of home. The fraction of time spent drops off sharply with log 10 distance. The 1096 spatial scales also depend on transmission intensity. In places with highly heterogeneous 1097 transmission, places with the highest transmission intensity, will have the greater the 1098 fraction of transmission that occurs at short distances. 1099 38 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Mosquitoes, T ravel, and T ransmission Highly spatially resolved data describing 1100 the EIR is rarely available. It is often cheaper, albeit less accurate, to use cross-sectional 1101 blood survey data describing malaria prevalence ( i.e. the parasite rate, PR) to estimate 1102 local transmission. Spatial models and spatial metrics described herein provide some 1103 guidelines about how patterns in the PR can be used to identify areas with the most 1104 mosquitoes, particularly given the enormous heterogeneity in human population density. 1105 It is commonly assumed that local clustering of cases implies that there is local 1106 transmission. For models developed in this framework, the vectorial capacity matrix 1107 (Eq. 39) describes parasite dispersion by mosquitoes, and evidence suggests that the 1108 spatial scales describing parasite dispersal by mosquitoes could vary by context [62]. 1109 Importantly, imported malaria can confound the relationship between local 1110 transmission by mosquitoes and prevalence. Travel habits and other traits describing 1111 humans often cluster spatially, partly because human neighborhoods are organized by 1112 socio-economic status. Spatial clustering of cases could arise if travel habits and thus 1113 malaria importation rates are spatially clustered, giving the appearance of local 1114 transmission. 1115 Measuring Reproductive Success The most complete measure of transmission in 1116 an area is a reproductive number – the number of malaria cases arising from each 1117 malaria case after one complete parasite generation. We have defined reproductive 1118 matrices in several ways as matrices describing reproductive success among patches 1119 within a spatial domain, which can be used to define local reproductive numbers as 1120 cases arising from a patch. These reproductive matrices form a basis for investigating 1121 the appropriate spatial scales to measure and model transmission, for estimating 1122 contamination in randomized control trials, and for understanding the spatial effect 1123 sizes of control. These can put other data into a context that is relevant for 1124 transmission. For example, mosquito counts data and measures of malaria can vary over 1125 very short distances [28,62]. The functional relevance of local heterogeneity in mosquito 1126 catch counts or in malaria prevalence can be critically examined by examining a matrix 1127 that integrates the effects of parasite movement in both mosquitoes and humans. After 1128 fully considering the uncertainty, it may be possible to determine the relevant spatial 1129 scales of transmission and thus the relevant spatial units for estimating reproductive 1130 numbers for malaria dynamics and control. 1131

1132 The simplicity of the Ross-Macdonald model can be contrasted with Hackett’s 1133 description of the elaborate and context-dependent nature of malaria that he observed 1134 in the field [27]: 1135 . . .malaria is so moulded and altered by local conditions that it becomes a 1136 thousand different diseases and epidemiological puzzles. Like chess, it is 1137 played with a few pieces, but is capable of an infinite variety of situations. 1138 The Ross-Macdonald model clearly identified enough chess pieces to develop basic 1139 concepts and theory to describe and measure malaria transmission [1], such as vectorial 1140 39 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint capacity, the basic reproductive numbers, daily human biting rates, sporozoite rates, 1141 entomological inoculation rates, and malaria parasite rates ( i.e. prevalence). These 1142 basic metrics have formed the basis for quantitative studies of malaria transmission, but 1143 they ignored heterogeneity and complexity. In particular, the metrics and associated 1144 concepts describing parasite dispersal in infected mosquitoes and humans were missing. 1145 Parasite dispersal is defined by the locations where infecting bites occurred in chains 1146 of transmission, tracing dispersal events backwards through alternating jumps in 1147 moving, infected humans and mosquitoes. It is practically impossible to study 1148 transmission directly, but this framework has established a quantitative basis for 1149 studying transmission through a set of constructs describing closely related processes 1150 that can be observed. We have established a basis for describing dispersal rigorously, 1151 and for analyzing dispersal and simulating transmission. The metrics and concepts we 1152 have proposed here are designed to quantify transmission (and uncertainty about 1153 transmission) through the study of patterns and the processes that generated them. 1154 The metrics provide a rigorous way of quantifying parasite dispersal and spatial 1155 transmission intensity. 1156 In developing models of a specific place for monitoring and evaluating malaria, it is 1157 important to understand where and when transmission occurs as well as the local 1158 contextual factors that shape transmission. In the Ross-Macdonald model, the basic 1159 notions of reproductive success, transmission, and community effect sizes of control were 1160 based on the abstract notion of a population, but it was never clear how to define a 1161 population for purposes of quantifying malaria transmission dynamics: “What, if 1162 anything, is a malaria population?” Focal transmission has been described [63], but 1163 without a quantitative basis for quantifying malaria spatial heterogeneity and spatial 1164 dynamics, there was no basis for a nuanced quantitative discussion about “What, if 1165 anything, is a focus?” Without defining explicit boundary conditions, it was easy to 1166 ignore malaria importation: “What fraction of malaria in a defined area was 1167 attributable to local transmission?” Without modeling structured populations, it was 1168 impossible to understand how differences in human behaviors would affect 1169 transmission [64]. Who is responsible for most local transmission or malaria 1170 importation? In malaria control, these discussions have focused on the issue of 1171 stratification, but it remains unclear whether those strata should define sub-populations, 1172 spatial areas, or both. Without a framework for understanding malaria transmission 1173 spatially in heterogeneous populations, it was difficult to develop a consistent 1174 methodology for quantifying transmission in a specific place and time. 1175 We have synthesized a set of old, new, and revised models to fully develop concepts, 1176 constrain parameters, and update basic concepts and metrics in a spatial context. New 1177 algorithms have filled a need to connect model parameters with data and remove bias 1178 while guaranteeing mathematical consistency. The new framework and spatial metrics 1179 make model complexity scalable, and it provides a way to study the role of context in 1180 mosquito ecology and malaria transmission. How and why do bionomic parameters vary 1181 over space and time? What spatial scales characterize mosquito populations? What are 1182 the appropriate spatial scales to measure transmission and intervention coverage as a 1183 spatial average? What are some appropriate methods for dealing with population 1184 heterogeneity, including heterogeneity arising from differences in behavior, exposure, or 1185 immunity? 1186 40 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint The framework emphasizes the way we organize our knowledge about malaria into 1187 bins of expertise. Given the complexity of the problem, this means modellers can build 1188 models that adapt over time as more information about transmission in a place 1189 accumulated. The first models can focus on components whose dynamics are better 1190 known, and use simpler, pragmatic approaches to parts of the model whose mechanistic 1191 foundations are more uncertain. Our framework can make building ensembles of 1192 plausible models to cover this uncertainty easier. Model building and model comparison 1193 makes it possible to weigh the importance of various factors in context. In asking where 1194 transmission is occurring, we are concerned about mosquito populations, human 1195 behaviors, and human blood feeding. In asking who is responsible for malaria, we are 1196 not just concerned about differences in infectiousness, but also populations who import 1197 malaria, and strata who play an out-sized role in moving malaria around an area. These 1198 are the basic quantities that play a role in spatial targeting and in tailoring 1199 interventions to context. 1200

1201 The goal of this study was to develop and present framework that could support model 1202 building for planning, monitoring and evaluating malaria control programs. Suites of 1203 models developed in this framework can be used to synthesize data, to quantify the 1204 major factors affecting transmission in a particular place, to identify critical data gaps 1205 and prioritize new data collection, to propagate uncertainty through analyses, and to 1206 support policy. We plan to use the framework to synthesize evidence and to give robust 1207 policy advice about malaria control on Bioko Island, and elsewhere, iteratively as part 1208 of adaptive malaria control. The spatial metrics and concepts describe an important 1209 dimension of malaria transmission that can help tailor interventions and spatially target 1210 interventions. In future studies, we plan to address concerns about the temporal 1211 dimensions of transmission, including threshold conditions, forcing by weather, and the 1212 spatial dimensions of malaria control. In adaptive management, the goal is to support 1213 monitoring and evaluation by developing rigorous methods that quantify malaria 1214 transmission as a changing baseline (e.g., forced by weather) that has been modified by 1215 control. In other settings, the framework can be used to enhance the design of 1216 randomized control trials or to help programs implement and interpret ad hoc 1217 experiments to fill local knowledge gaps. Simulation-based analytics in this framework 1218 can be updated using evidence collected by malaria programs to update models and 1219 analysis and revise policy recommendations, to target and tailor interventions, and to 1220 use evidence to adapt to changing local conditions. 1221 Acknowledgments 1222 DLS and SLW were supported by a grant from the National Science Foundation, 1223 Directorate for Technology, Innovation, and Partnerships (TIP) as part of the 1224 Convergence Accelerator Program (NSF 2040688). DLS, JMH, and AC were supported 1225 by a grant from the National Institute of Allergies and Infectious Diseases (R01 1226 AI163398). DLS was supported by grant from the National Institutes of Health 1227 41 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint (2U19AI089674). DLS, DMS, and JNN, were supported by a grant from the Bill and 1228 Melinda Gates Foundation (INV 030600). This project is partly a product of discussions 1229 with the mosquito working groups of working groups over several years RAPIDD 1230 (Research and Policy for Infectious Disease Dynamics). Over that time, many of us 1231 benefitted from the unwavering support and inspiration of F. Ellis McKenzie. OJB was 1232 supported by a UK Medical Research Council Career Development Award 1233 (MR/V031112/1). 1234 42 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Bioko Island: Sectors, Areas, Regions W. Malabo Malabo E. Malabo East CoastWest Coast Luba Riaba South Fig 8. An important practical concern is spatial granularity of patches for simulation-based studies. For Bioko Island, Equatorial Guinea, for example, we could define patches at several scales: the whole island; or approximately 240 occupied areas (1km× 1km, the squares); approximately 4 , 400 occupied 100m× 100m sectors (points); or 8 distinct regions (the colors of the squares); or clusters of contiguous sectors (the colors of the points); or approximately 70,000 individual households. An important concern is that the weight of evidence – the number of observations per patch – declines sharply as granularity of the simulations increases. This framework makes it possible to define a set of nested (or partially nested) studies that modify the number and size of patches, which requires modifying the human and mosquito mobility sub-models, but that holds other aspects of the model constant. 43 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Fig 9. (Left): bulk transmission metric describing transmission from the most densely populated area in Malabo, the capitol city, seen as the bright cell in the Northern tip of the island, to all other populated areas. (Right): bulk transmission from the most highly populated area in the south of the island (Luba), seen as the bright cell in the small harbor on the Western coast of the island. Fig 10. Structural Elements of the framework are flexible to facilitate building models that are appropriate for various settings. These diagrams illustrate two examples. left) A forest malaria model with seven patches (including 3 villages and 2 campsites), 6 population strata, and 5 aquatic habitats. The village residents are stratified into loggers and other residents. Loggers from different villages spend time at home or in campsites, which have no permanent residents. Aquatic habitats (the moons) can be in villages, in campsites, or in patches near villages. Some villages ( e.g. village 3), could lack mosquitoes but still have populations at risk. Right) It is also possible to model indoor and outdoor blood feeding with indoor and outdoor patches that share the same place. In these models, movement indoors vs. outdoors in the same place is modeled differently from movement among outdoor patches. 44 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint Supplements 1235 List of Supplements 1236 1. Supplement 1 - a github repository https://dd-harp.github.io/exDE/ 1237 2. Supplement 2 - Glossary 1238 3. Supplement 3 - Modular Notation 1239 4. Supplement 4 - Human Travel and Mobility 1240 5. Supplement 5 - Vector Dynamics 1241 45 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint

1. Smith DL, Battle KE, Hay SI, Barker CM, Scott TW, McKenzie FE. Ross, Macdonald, and a theory for the dynamics and control of mosquito-transmitted pathogens. PLoS Pathog. 2012;8(4):e1002588. 2. Reiner RC Jr, Perkins TA, Barker CM, Niu T, Chaves LF, Ellis AM, et al. A systematic review of mathematical models of mosquito-borne pathogen transmission: 1970-2010. J R Soc Interface. 2013;10(81):20120921. 3. Smith DL, Perkins TA, Reiner RC Jr, Barker CM, Niu T, Chaves LF, et al. Recasting the theory of mosquito-borne pathogen transmission dynamics and control. Trans R Soc Trop Med Hyg. 2014;108(4):185–197. 4. Levin S. The Problem of Relevant Detail. In: Busenberg S, Martelli M, editors. Differential Equations Models in Biology, Epidemiology and Ecology: Proceedings of a Conference held in Claremont California, January 13–16, 1990. Berlin, Heidelberg: Springer Berlin Heidelberg; 1991. p. 9–15. 5. Radcliffe J. The initial geographical spread of host-vector and carrier-borne epidemics. J Appl Probab. 1973;10(4):703–717. 6. Bailey NTJ. The Biomathematics of Malaria. Oxford: Charles Griffin & Company Ltd.; 1982. 7. Dye C, Hasibeder G. Population dynamics of mosquito-borne disease: effects of flies which bite some people more frequently than others. Trans R Soc Trop Med Hyg. 1986;80(1):69–77. 8. Hasibeder G, Dye C. Population dynamics of mosquito-borne disease: persistence in a completely heterogeneous environment. Theor Popul Biol. 1988;33(1):31–53. 9. Sattenspiel L, Dietz K. A structured epidemic model incorporating geographic mobility among regions. Math Biosci. 1995;128(1-2):71–91. 10. Torres-Sorando L, Rodriguez DJ. Models of spatio-temporal dynamics in malaria. Ecol Modell. 1997;104(2-3):231–240. 11. Jansen VA, Lloyd AL. Local stability analysis of spatially homogeneous solutions of multi-patch systems. J Math Biol. 2000;41(3):232–252. 12. Auger P, Kouokam E, Sallet G, Tchuente M, Tsanou B. The Ross–Macdonald model in a patchy environment. Math Biosci. 2008;216(2):123–131. 13. Adams B, Kapan DD. Man bites mosquito: understanding the contribution of human movement to vector-borne disease dynamics. PLoS One. 2009;4(8):e6763. 14. Cosner C, Beier JC, Cantrell RS, Impoinvil D, Kapitanski L, Potts MD, et al. The effects of human movement on the persistence of vector-borne diseases. J Theor Biol. 2009;258(4):550–560. 15. Gao D, Ruan S. A multi-patch malaria model with logistic growth populations. SIAM J Appl Math. 2012;72(3):819–841. 46 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint 16. Perkins TA, Scott TW, Le Menach A, Smith DL. Heterogeneity, mixing, and the spatial scales of mosquito-borne pathogen transmission. PLoS Comput Biol. 2013;9(12):e1003327. 17. Ruktanonchai NW, DeLeenheer P, Tatem AJ, Alegana VA, Caughlin TT, Zu Erbach-Schoenberg E, et al. Identifying Malaria Transmission Foci for Elimination Using Human Mobility Data. PLoS Comput Biol. 2016;12(4):e1004846. 18. Citron DT, Guerra CA, Dolgert AJ, Wu SL, Henry JM, S´ anchez C HM, et al. Comparing metapopulation dynamics of infectious diseases under different models of human movement. Proc Natl Acad Sci U S A. 2021;118(18). 19. Guerra CA, Citron DT, Garc´ ıa GA, Smith DL. Characterising malaria connectivity using malaria indicator survey data. Malar J. 2019;18:440. 20. Smith NR, Trauer JM, Gambhir M, Richards JS, Maude RJ, Keith JM, et al. Agent-based models of malaria transmission: A systematic review. Malar J. 2018;17(1):299. 21. Grillet ME, Hern´ andez-Villena JV, Llewellyn MS, Paniz-Mondolfi AE, Tami A, Vincenti-Gonzalez MF, et al. Venezuela’s humanitarian crisis, resurgence of vector-borne diseases, and implications for spillover in the region. Lancet Infect Dis. 2019;19(5):e149–e161. 22. Jaramillo-Ochoa R, Sippy R, Farrell DF, Cueva-Aponte C, Beltr´ an-Ayala E, Gonzaga JL, et al. Effects of political instability in Venezuela on malaria resurgence at Ecuador-Peru border, 2018. Emerg Infect Dis. 2019;25(4):834–836. 23. Cohen JM, Smith DL, Cotter C, Ward A, Yamey G, Sabot OJ, et al. Malaria resurgence: a systematic review and assessment of its causes. Malar J. 2012;11:122. 24. Plowe CV, Kublin JG, Dzinjalamala FK, Kamwendo DS, Mukadam RAG, Chimpeni P, et al. Sustained clinical efficacy of sulfadoxine-pyrimethamine for uncomplicated falciparum malaria in Malawi after 10 years as first line treatment: five year prospective study. BMJ. 2004;328(7439):545. 25. Epstein A, Maiteki-Sebuguzi C, Namuganga JF, Nankabirwa JI, Gonahasa S, Opigo J, et al. Resurgence of malaria in Uganda despite sustained indoor residual spraying and repeated long lasting insecticidal net distributions. PLOS Global Public Health. 2022;2(9):e0000676. 26. Hackett LW, Missiroli A. The natural disappearance of malaria in certain regions of Europe. Am J Epidemiol. 1931;13(1):57–78. 27. Hackett LW. Malaria in Europe. An Ecological Study. London, Oxford University Press; 1937. 28. Bejon P, Williams TN, Nyundo C, Hay SI, Benz D, Gething PW, et al. A micro-epidemiological analysis of febrile malaria in Coastal Kenya showing hotspots within hotspots. Elife. 2014;3(2-3):e02130. 47 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint 29. Ross R. Some quantitative studies in epidemiology. Nature. 1911;87(2188):466–467. 30. Macdonald G. The analysis of equilibrium in malaria. Trop Dis Bull. 1952;49(9):813–829. 31. Smith DL, Perkins TA, Tusting LS, Scott TW, Lindsay SW. Mosquito population regulation and larval source management in heterogeneous environments. PLoS One. 2013;8(8):e71247. 32. Brady OJ, Godfray HCJ, Tatem AJ, Gething PW, Cohen JM, McKenzie FE, et al. Adult vector control, mosquito ecology and malaria transmission. Int Health. 2015;7(2):121–129. 33. Killeen GF, McKenzie FE, Foy BD, Bøgh C, Beier JC. The availability of potential hosts as a determinant of feeding behaviours and malaria transmission by African mosquito populations. Trans R Soc Trop Med Hyg. 2001;95(5):469–476. 34. Ma BO, Roitberg BD. The role of resource availability and state-dependence in the foraging strategy of blood-feeding mosquitoes. Evol Ecol Res. 2008;10(8):1111–1130. 35. Yakob L. How do biting disease vectors behaviourally respond to host availability? Parasit Vectors. 2016;9:468. 36. Thongsripong P, Hyman JM, Kapan DD, Bennett SN. Human-mosquito contact: A missing link in our understanding of mosquito-borne disease transmission dynamics. Ann Entomol Soc Am. 2021;114(4):397–414. 37. Service, M W. Mosquito (Diptera: Culicidae) dispersal–the long and short of it. J Med Entomol. 1997;34(6):579–588. 38. Smith DL, Dushoff J, McKenzie FE. The risk of a mosquito-borne infection in a heterogeneous environment. PLoS Biol. 2004;2(11):e368. 39. Wu SL, S´ anchez C HM, Henry JM, Citron DT, Zhang Q, Compton K, et al. Vector bionomics and vectorial capacity as emergent properties of mosquito behaviors and ecology. PLoS Comput Biol. 2020;16(4):e1007446. 40. Brady OJ, Godfray HCJ, Tatem AJ, Gething PW, Cohen JM, McKenzie FE, et al. Vectorial capacity and vector control: reconsidering sensitivity to parameters for malaria elimination. Trans R Soc Trop Med Hyg. 2016;110(2):107–117. 41. Smith DL, Drakeley CJ, Chiyaka C, Hay SI. A quantitative analysis of transmission efficiency versus intensity for malaria. Nat Commun. 2010;1(8):108. 42. Le Menach A, Tatem AJ, Cohen JM, Hay SI, Randell H, Patil AP, et al. Travel risk, malaria importation and malaria transmission in Zanzibar. Sci Rep. 2011;1:93. 48 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint 43. R Core Team. R: A Language and Environment for Statistical Computing; 2022. 44. Soetaert K, Petzoldt T, Setzer RW. Solving Differential Equations in R: Package deSolve; 2010. 45. Citron DT, Guerra CA, Garc´ ıa GA, Wu SL, Battle KE, Gibson HS, et al. Quantifying malaria acquired during travel and its role in malaria elimination on Bioko Island. Malar J. 2021;20:359. 46. Rund SSC, O’Donnell AJ, Gentile JE, Reece SE. Daily rhythms in mosquitoes and their consequences for malaria transmission. Insects. 2016;7(2). 47. Monroe A, Moore S, Okumu F, Kiware S, Lobo NF, Koenker H, et al. Methods and indicators for measuring patterns of human exposure to malaria vectors. Malar J. 2020;19(1):207. 48. Monroe A, Msaky D, Kiware S, Tarimo BB, Moore S, Haji K, et al. Patterns of human exposure to malaria vectors in Zanzibar and implications for malaria elimination efforts. Malar J. 2020;19(1):1–14. 49. Aron JL, May RM. The population dynamics of malaria. In: Anderson RM, editor. The Population Dynamics of Infectious Diseases: Theory and Applications. Boston, MA: Springer US; 1982. p. 139–179. 50. Le Menach A, McKenzie FE, Flahault A, Smith DL. The unexpected importance of mosquito oviposition behaviour for malaria: Non-productive larval habitats can be sources for malaria transmission. Malar J. 2005;4(1):23. 51. N˚ asell I. Hybrid Models of Tropical Infections. vol. 59 of Lecture Notes in Biomathematics. 1st ed. Berlin: Springer-Verlag; 1985. 52. Johnston GL, Smith DL, Fidock DA. Malaria’s missing number: calculating the human component of R 0 by a within-host mechanistic model of Plasmodium falciparum infection and transmission. PLoS Comput Biol. 2013;9(4):e1003025. 53. Guerra CA, Kang SY, Citron DT, Hergott DEB, Perry M, Smith J, et al. Human mobility patterns and malaria importation on Bioko Island. Nat Commun. 2019;10:2332. 54. Darroch JN, Seneta E. On quasi-stationary distributions in absorbing continuous-time finite Markov chains. J Appl Probab. 1967;4(1):192–196. 55. Wu SL, Bennett JB, S´ anchez C HM, Dolgert AJ, Le´ on TM, Marshall JM. MGDrivE 2: A simulation framework for gene drive systems incorporating seasonality and epidemiological dynamics. PLoS Comput Biol. 2021;17(5):e1009030. 56. Macdonald G. The analysis of the sporozoite rate. Trop Dis Bull. 1952;49(6):569–586. 57. Smith DL, McKenzie FE. Statics and dynamics of malaria infection in Anopheles mosquitoes. Malar J. 2004;3:13. 49 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint 58. Garrett-Jones C. Prognosis for interruption of malaria transmission through assessment of the mosquito’s vectorial capacity. Nature. 1964;204:1173–1175. 59. Smith DL, Musiime AK, Maxwell K, Lindsay SW, Kiware S. A New Test of a Theory about Old Mosquitoes. Trends Parasitol. 2021;37(3):185–194. 60. van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. 2002;180(1-2):29–48. 61. Diekmann O, Heesterbeek JAP, Metz JAJ. On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations. J Math Biol. 1990;28(4):365–382. 62. Carter R, Mendis KN, Roberts D. Spatial targeting of interventions against malaria. Bull World Health Organ. 2000;78(12):1401–1411. 63. Bousema T, Griffin JT, Sauerwein RW, Smith DL, Churcher TS, Takken W, et al. Hitting hotspots: spatial targeting of malaria for control and elimination. PLoS Med. 2012;9(1):e1001165. 64. Sturrock HJW, Hsiang MS, Cohen JM, Smith DL, Greenhouse B, Bousema T, et al. Targeting asymptomatic malaria infections: active surveillance in control and elimination. PLoS Med. 2013;10(6):e1001467. 50 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted November 15, 2022. ; https://doi.org/10.1101/2022.11.07.22282044doi: medRxiv preprint