Spatial scales in human movement between reservoirs of infection

2.1. Locations of infectious focal points and households

There are important ethical implications to the public availability of data on both the spatial locations and movement patterns of individuals over various scales in time, especially in low-and middle-income countries [30]. Among these are data privacy and the potential identifiability of individuals through the use of ‘big data’ by governments, companies and universities, and the potential use of this information for non-research purposes. Therefore, access to such high-quality data for study is scarce. This represents a challenge for those concerned with the various consequences of such movements and their methods of model validation. In this paper, we present mathematical models that can be used in the absence of access to high-quality movement data. Model dynamics are implemented on the spatial patterns of, e.g., building locations, obtained from real-world datasets, e.g., the high resolution settlement layer dataset generated by the Facebook Connectivity Lab [31]. Our methodology should allow researchers to draw useful conclusions on the impact of human spatial movement on NTD programmes from readily available geospatial data (coordinates of households, settlements and environmental features, such as water bodies). Other data sources, e.g., call data records (from mobile phones) and data from migration questionnaires in control studies, may be available under certain circumstances and can easily be integrated in our modelling framework. The models we introduce here should also be able to take this specific case-study information into account by appropriate parameter inference.

Let us first consider the embedding of spatial locations of households as points in a 2-dimensional flat space. In addition to these household points there will also exist focal points of infectious contact that depend on the helminth species, e.g., near buildings, water sources, parts of farmland with large concentrations of infectious material and households themselves. Choosing a particular focal point, we will use the observed spatial distributions of buildings in the high resolution settlement layer dataset [31] as a proxy for the spatial distribution of household locations around it. Therefore, the expectation of having n_<r nearby households or buildings within a radial distance of r from the focal point, i.e., E(n_<r), scales according to the following radially-dependent power law It is important to note that, since at first we consider only a single focal point from which household distances is measured, the distribution of focal points themselves is not relevant to the calculation. The spatial distribution of focal points will, however, become relevant later on when we cluster them into reservoirs of infection and investigate their spatial extent. Note, however, that households themselves may be a reasonable tracer of the underlying distribution for infectious focal points of STH [32, 33], though our calculations throughout will not depend on this possibility.

In the cases of STH and schistosomiasis transmission, although it is clear that the location and importance of their infectious focal points can vary over time — e.g., a chance decline in snail population for a particular location at the edge of a fresh water source leads to a drop in infectious spread — throughout this work, we shall assume that their spatial distribution is effectively temporally stationary over the timescales of interest. Note that for some NTDs, the infectious contact event locations can also be mobile, e.g., in the case of lymphatic filariasis transmission, where the larval stage of microfilariae enter the body by bites from infected mosquitoes — hence, the contact event locations will migrate with the human and mosquito populations.

By averaging over the observed point spread patterns in the high resolution settlement layer dataset, we obtain the following cumulative neighbouring point number distributions as a function of radial distance. This can be compared to the radial dependence of the power-law index α(r) we have quoted in Eq. (1). In order to give an indication of the effect of spatial heterogeneity in the distances between buildings for this averaged quantity, we have plotted some example lines (left column of plots) and the sample mean lines with root-mean-square deviation (RMSD) shaded regions of α(r) (right column of plots) from individual randomly sampled initial building locations for different countries in Fig. 2. In Central Malawi and Northern Benin, up to some substantial variance, below a transitioning radial scale r < r_tr, the index seems to have a mean of α ≃1 and above this scale r > r_tr the mean of the index begins a transition to α → 2. In the Central Ivory Coast, we see that the mean trend rapidly varies between α ≃ 1.5 and 0.5 below the transition scale r < r_tr, and then at radial distances much greater than this scale r» r_tr most of the power-laws appear to begin a transition to α →2. The substantial variation below the transition scale in the Ivory Coast is likely due to a greater degree of clustering on particularly small scales (note the drop towards α ≃ 0 before the transition scale indicates this as very few new points are added to the cumulative total for this range).

• Download figure
• Open in new tab

Figure 2.

Left column: The cumulative number of buildings (vertical axis) below a given radial distance (horizontal axis) from sample locations in Malawi, Benin and the Ivory Coast obtained from the high resolution settlement layer dataset generated by the Facebook Connectivity Lab [31]. Lines for different constant α and r_µ values are provided only as rough indicators for comparison. Right column: The sample mean lines and root-mean-square deviation (RMSD) shaded regions of the power law index α plotted against radial distance in each case. The RMSD was calculated by computing discrete derivatives between log-frequency bins as a function of log-radial distance. In each country, the collection of sample locations was randomly drawn from: the central region of Malawi, including Lilongwe; the wide northern region of Benin, spanning between Kandi and Djougou; and a wide central region of the Ivory Coast, including Yamoussoukro.

For all countries and initial locations in Fig. 2, we can confirm that the majority of the power-laws begin to converge to a mean of α → 2 above a particular transition scale r_tr which varies substantially between location. This convergence reflects the emergence of statistical spatial homogeneity only on very large spatial scales, but this is clearly not present below r_tr in all cases. The power-law signature of statistical homogeneity being α = 2 should become clear by noting that another way to generate it is to draw a Poisson point process with an intensity rdr ∝ r² (the flat Euclidean measure in 2 dimensions). Note that statistical isotropy is also implicitly confirmed by choosing different r = 0 locations for each line and observing the same behaviour up to some noise that is captured by the RMSD in the right column of plots. Deviations from this trend towards α = 2 (particularly in the Ivory Coast) are also present in many other locations in the same countries, due to the presence of substantial global heterogeneity, e.g., through proximity to large voids in the distribution of points, which exist for geographic reasons.

2.2. Obtaining a model using maximum entropy arguments

In this section we provide a simple argument to generate the broken power-law frequency distribution of work-home commute distances (as observed and described in, e.g., Kung et al. [24] and Zhao et al. [26]) using maximum entropy arguments. The distributions of movement distances we derive will be used to model the migration of individuals from their house-hold locations to a focal point of infection. In the absence of access high-quality site-specific data that we discussed in the previous section, we are using the work-home commuter distance distributions as a guide for the likely patterns of movement. Work-home commuter distances should be a reasonable proxy for the overall movement of individuals in many areas. However, the model that we develop in this section is adaptive and flexible to modification, even in circumstances where this proxy is inappropriate, we will be able to vary distribution parameters to quantify their impact on our conclusions.

For greater analytic insight in the expressions we derive in this section, and throughout most of the subsequent arguments made in this article, we shall now assume that Eq. (1) has a constant index α(r) = α. When considering the heterogeneity in the real power-law behaviour below r_tr in Fig. 2, this assumption will not be appropriate in all locations/situations. Note, however, the results we obtain in all cases may be easily numerically generalised to take into account complex transitions between different power-law behaviours, as are observed in Fig. 2. In addition, choosing several key values for constant α (which fit within the range of observed power law behaviours) in what follows will aid with quantifying its impact on our conclusions.

Assuming that each individual jumps to a single new location and then returns back to their starting position over the course of each day, we may follow a simple argument that incidentally generates a distribution over travel distances which is consistent with those observed and described in, e.g., Kung et al. [24] and Zhao et al. [26]. Let us take the distribution of neighbouring points to be isotropic. Defining p(r) as the probability density function (PDF) of an individual jumping to a given location as a function of its radial distance, one may then transform the flat radial probability measure p(r)dr by a Jaco-bian v(r) = dE (n_<r) / dr like so which accounts for the power-law density of points that may depart from linear scaling according to Eq. (1). Another way to view this rescaling is to consider the expectation of the probability measure for a radially-generated, statistically isotropic Poisson point process.

Note that if we additionally assume that the PDF of the jump distribution for each individual has a mean σ then the maximum entropy PDF is that of an exponential distribution p(r) = p(r; σ) = Exp(r; 1/σ), then by using Eq. (2) we can trivially demonstrate that the resulting normalised PDF is a Gamma distribution One possible interpretation for σ here is a number, with physical units of distance, which quantifies the access of that individual to varying degrees of transportation methods that are themselves capable of travelling varying distances. To clarify, without access to any motorised transportation, an individual may receive a low value for σ, and in contrast, if an indvidiual has access to a vehicle then they may receive a larger σ value. Despite this simple suggested explanation, the presence of variable scales for human travel distance distributions is not contingent on access to vehicular transportation at all. Indeed, human walking patterns themselves are known to exhibit multi-distance-scale characteristics in many contexts — see, e.g., Ref. [34].

Assuming that 1/σ values may vary across the population of individuals with a distribution of known mean E(1/σ) and expected scaling E(ln σ), the maximum entropy distribution for values of 1/σ is another gamma distribution. Choosing specific parameter values, we therefore have where the resulting marginal jump PDF is obtained by integration in the following way where Beta,(.;.,.,.) denotes the generalised beta prime distribution (or beta distribution of the second kind), which takes a form similar to that of a beta distribution. Note that, to avoid a divergent integral, one must specify that α, β, r_µ ∈ ℝ ⁺.

The main motivations for using maximum entropy probability distributions in the case above (and in general) include the beneficial properties that such distributions have when limited a priori information is available. In the present instance, these properties translate into distributions which account for the maximum amount of potential variability in outcomes while still satisfying the basic pre-defined classifying properties, e.g., the gamma distribution we have used above is maximally entropic under the condition that the mean and expected scaling of the distribution are known. Furthermore, it is often reasonable to expect that the steady states of many systems observed across the sciences are well-described by maximum entropy configurations.

We introduced another radial spatial scale r_µ (distinct from r_tr in the previous section) in Eq. (4) whose effect, in conjunction with α and β, we shall explore. In Fig. 3 we have plotted the marginal jump PDF given by Eq. (5) for a range of parameter values α, β and r_µ. By inspection of Eq. (5), in the limit where r « r_µ the distribution scales as ∼ r^α−1, whereas in the opposite limit r » r_µ the distribution exhibits a scaling ∼ r^−β−1. The scale of r_µ therefore plays the role of separating two regimes in the distribution of jump lengths for a given individual. For radial distances below r_µ, the geometric distribution of available spatial locations to jump to (which is encoded in the α parameter) dominates the behaviour, whereas, for radial distances above r_µ, the predisposition of individuals to travel a given distance combines with the geometry of points to give the distribution behaviour and this is encoded in the β parameter of the second gamma distribution in Eq. (4).

• Download figure
• Open in new tab

Figure 3.

Plots of the marginal probability density p(r; α, β, r_µ) as a function of radial distance (in units of km) generated using Eq. (5) for a range of α, β and r_µ values for comparison.

So far we have not assumed that the parameter r_µ in the distribution over inverse-jump scales 1/σ given by Eq. (4) is related to any purely geometrically-driven break in the power laws of Fig. 2. By comparing our observed power laws for the Ivory Coast with the observed work-home commuter distance distributions in Ref. [24], for instance, we see that r_µ « r_tr in general. This is perhaps unsurprising since the scale r_µ may arise due to the availability of easy access to vehicles that allow for longer-distance travel for either work or school, or water sources near a given settlement. Despite this fact, it is intriguing to note that the radial scale associated to the peak in clustering (α ≃ 1.5 well below r_tr) for the Ivory Coast dataset in Fig. 2 occurs at the same rough scale that is observed for r_µ in Ref. [24]. We leave the investigation into the potential link between these two scales for future work.

Despite statistical isotropy and homogeneity on large spatial scales, we have observed that the variation between the power-laws for cumulative building numbers observed in Fig. 2, and hence the α index of Eq. (1), exhibits a large degree of location dependence between regions. In Fig. 4, we see the effect that this small-scale heterogeneity in the distances between buildings has on the probability density for our one-jump model. We plot a comparison between the binned frequency of 10⁴ individuals (with exponential jump distributions each with a scale drawn from Eq. (4)) simulated to travel (solid lines) between two buildings on real-world map data. The data for building locations corresponds to the same regions of each country as in Fig. 2, and the marginalised jump probability densities derived from Eq. (5) for choices of α = 1, 2. In all cases we have fixed β = 2 to a single value for easier comparison of different lines; this value is approximately the same as the tails observed in the distance distributions of Kung et al. [24].

• Download figure
• Open in new tab

Figure 4.

The binned frequency of 10⁴ individuals simulated to travel a between two buildings (solid lines) whose locations correspond to data from the same regions of each country as in Fig. 2 (where some globally-applied uniform random point density thinning has been used to reduce computational load of exploring the distribution tails). This corresponds to the same marginalised jump probability density as Eq. (5) but instead of using a geometric power-law given by Eq. (1) we have implemented our single jumps on real-world map data (with exponential jump distributions each with a scale drawn from Eq. (4)). For comparison, we have also plotted some jump probabilities calculated using Eq. (5) with α = 1 (dashed lines) and α = 2 (dotted lines). In all cases we have fixed β = 2.

The agreement between the α = 1, 2 curves and the distances generated from single simulated jumps between buildings on the real-world map data appears to be variable between countries. Overall, it is unsurprising that the shapes of the curves all seem to agree for scales r > r_µ, because the asymptotic ∼ r^β−1 behaviour dominates in all cases, which is consistent with our model. In the opposing limit, however, the small-scale variability in power laws for each country dominates and the departure from our constant-α models appears more pronounced — in favour of more transitory distance-dependent behaviour. Note also that the range in variability that the region between our α = 1 and α = 2 lines can account for is quite substantial. Such a range in behaviour motivates the use of both α = 1 and α = 2 in quoting the most important results that follow so as to correctly assess the impact of small-scale region-specific building distributions.

3.1. Homogeneous randomly-directed jumps

Following along similar lines to Sec. 2.2, let us now consider an individual who completes multiple successive jumps with exponentially-distributed lengths r_i (where we have indexed each successive jump event with an i) with a fixed mean scale σ within a day. For a fixed jump rate 𝒥 in time, consider the following compound Poisson process for a two-dimensional vector x(t) which encodes the 2-dimensional Euclidean coordi-nate position of an individual over time where U(Δθ_i; −π, π) is the uniform distribution PDF (assuming statistical isotropy) over the change in angle and Θ(·) is a Heaviside function. In Eq. (6) above, note that we are using an indicator function 𝟙_A(t) which takes value unity if t ∈ A, else zero.

The process specified by Eq. (6) assumes isotropy of both the geometric distribution of the points and a uniform-random direction choice of the individual. The distribution p[x(t)] is not known as a closed-form expression, however, it clearly has both a vanishing first moment and a second moment which scales according to σ² Furthermore, with a scale distribution over jump lengths applied — as in Eq. (5) — the second moment given in Eq. (10) becomes In the left panel of Fig. 5 we have plotted the binned frequency of total daily distances travelled by a population of 5 × 10⁴ individuals following the process defined by Eq. (6) and drawing each individual’s inverse-jump scale 1/σ from Eq. (4). We have also fixed 𝒥 = 1, 20 per day (in red and black lines, respectively), α = 1 and β = 1, 2, 3 as indicated by the increasing opacity within each triplet of lines. Contrasting these distributions, it is immediately clear that by increasing the jump rate 𝒥, the effect on the distribution of daily distances travelled is similar to increasing the value of α — see the middle panel of Fig. 3 for comparison.

• Download figure
• Open in new tab

Figure 5.

Numerical plots of the distance distributions generated by the multi-jump processes introduced in Sec. 3.1 with randomly-directed (left panel) and in Sec. 3.2 with unidirectional jumps (right panel). In the left panel, a population of 5 ×10⁴ individuals following the process defined by Eq. (6) have been drawn (drawing from Eq. (4) for their jump scale predispositions) and the binned frequency of the total distance evaluated at the end of the day with jump rates 𝒥 of = 1, 20 per day (in red and black lines, respectively), α = 1 and β = 1, 2, 3 as indicated by the increasing opacity within each triplet of lines. In the right panel, the marginal probability density given by Eq. (15) is depicted with jump rates of 𝒥 = 1, 10 per day (in red and black lines, respectively), α = 1 and β = 1, 2, 3 as indicated by the increasing opacity within each triplet of lines.

3.2. Homogeneous unidirectional jumps

Human movement patterns are not truly random, and instead, one might anticipate a strong directional dependence, e.g., to long distance travel of individuals over the course of multiple days to an intended destination. Due to this fact, let us try to quantify the effect of this behaviour by assuming the logical extreme, i.e., that the angular direction of an individual’s longdistance travel is fixed to a particular value. Following multiple successive jumps while assuming no variation in angle, Eq. (6) then becomes a compound Poisson process for the total jump distance x(t) ≡ |x(t) | in a single dimension. This process has a known characteristic function where ϕ_r(s) denotes the characteristic function of the stationary increments given by Eq. (3). Hence, we may write which has no closed-form inverse transformation, but is useful as an expression to calculate the moments. Alternatively, in order to find the distribution over x(t), we may marginalise over the number of jumps performed (which is Poisson-distributed by construction), which are themselves drawn from the jump distance distribution in the following way Therefore, if α = 1, Eq. (14) yields where ℐ _n(·) is the modified Bessel function of the first kind.

In the right panel of Fig. 5 we have plotted the marginal probability density given by integrating Eqs. (15) and (4) over 1/σ for direct comparison with the one-jump model of Eq. (5). We have used values of 𝒥 = 1, 10 per day (in red and black lines, respectively), α = 1 and β = 1, 2, 3 as indicated by the increasing opacity within each triplet of lines. Comparing this distribution with the one observed for randomly-directed multiple jumps (see Sec. 3.1), a similar, but even more pronounced, effect by varying 𝒥 is observed — since only a value of 𝒥= 10 is required to achieve roughly the same change. This is to be expected, as the unidirectional jumpers will always travel at least the same (or more often) a greater distance in total in comparison with the randomly-directed ones, which will likely induce a more severe deformation of the overall distribution when using the former.

The parameter degeneracy between α of the one-jump model and 𝒥 in the multi-jump processes described above suggests a test can be performed, to either validate a model, or discard it in favour of a modified version for a given situation with real data. By combining collected data on the total distance travelled per day and the local geometry of distances between buildings travelled to in a given region, the jump rate for each individual (if the same) could be statistically inferred. Comparing this inferred value to real data would provide a test of the movement models we have suggested in this work or potentially provide insight into where they may be improved to better reflect real human daily movement in a given setting.

The simple multi-jump models we have considered in this section have assumed that the distribution for the cumulative number of nearby buildings is homogeneous for each new jump. This approximation is not likely to work well on real-world map data as we have shown there is substantial small-scale heterogeneity in the distances between buildings exhibited in, e.g., Fig. 2. Furthermore, if the multi-jump movement model is purely diffusive over the building locations themselves (itself an assumption which must be investigated further), then there will be correlations in the location distributions between successive jumps that can lead to the distribution of jump distances exhibiting localising effects akin to those exhibited by powerlaw random banded matrices [35]. In other words, when an individual jumps from a cluster of points to, e.g., a point on the edge of their cluster, then they may experience a stronger ‘pull’ back towards the centre of this cluster on the next jump due to the anisotropy in their apparent local point density. This effect should be investigated in future work to evaluate its relevance to movement patterns.

A more sophisticated method to modify the model and deal with (uncorrelated) small-scale heterogeneity could be to sample from the numerically-obtained power laws from Fig. 2 to emulate the spatial heterogeneity in building distances on the real-world map data directly in the processes of Sec. 3. This method would have significant advantages, in terms of computational complexity, to direct simulation methods over the real world map data.

4.1. A coarse-grained stochastic reservoir network

Up until this point, the movement patterns we have discussed have not been contextualised in an NTD disease model. Since our current (static focal point) framework is best-suited to schistosomiasis and STH transmission, we shall use their well-studied transmission models as a reference throughout this section [36].

Let us now consider a group of focal points that are clustered into a ‘reservoir of infection’. Such a cluster is a source of infectious material from which disease-free individuals are exposed to new infections and, for those already with infections, to potentially increase them in intensity. With each of these clusters, let us now also associate a population of individuals with household locations that are sufficiently close such that the majority of their infectious inputs must be into one of their associated reservoir’s focal points. Note that the vague-ness in definition of ‘sufficiently close’ above immediately exposes some of the difficulties associated to defining an EU for the purpose of assessing local prevalence or intensity. Our work in the previous sections of this paper will allow us to proceed with a new approach to quantifying the spatial scales associated with the proximity of households to these focal points of infection. We will then use this approach to understand how spatial scales may be incoporated into the effective dynamical description of STH and schistosomiasis transmission.

Let us now denote the set of indicies which indentify all individuals who live in households that are less than a radial distance r_Λ away from a focal point that is included in the reservoir cluster as S _Λ. Additionally, to make a reservoir cluster, let us use r_Λ as a radial separation distance threshold below which any two focal points are clustered together to be part of the same infectious reservoir. For reservoirs of infection of either schistosomiasis or STH, assuming that the associated human population number N_Λ (the number of elements in S _Λ) does not change, the force of infection (FOI) Λ(t) at time t is updated according to where d_res is the death rate of the infectious material, d_w is the worm death rate, d is the human death rate, R₀ is the basic re-production number and l_i(t) is a compound Poisson process asociated to the time-dependent infectious material input of the i-th individual. Since helminths are dioecious, the definition for R₀, differs from its standard interpretation in standard microparasite diseases since the life cycle of these macroparasites depends on fertilisation of female worms (and hence the presence of both sexes) within a host [36]. Note also that in the equation above, and throughout, we shall neglect age structure in our description — though our formalism can be adapted to include this. This choice is for simplicity in presentation and will not affect our main conclusions.

By clustering focal points into a reservoir of infection by using a radial distance threshold, one might reasonably question the scaling in the spatial extent of the infectious focal points themselves. As was discussed in Sec. 3, the distribution of focal points of infection is likely to be both helminth species-dependent and variable according to local geographic considerations. Due to this variability, the spatial scales associated with clusters of focal points themselves will likely vary across a map. However, this fact will not directly affect the conclusions of the present work since we will consider distances between each individual focal point and its neighbouring households in turn. For specific case studies, the implementation of our algorithm for binding focal points together should be well-defined in most cases for STH and schistosomiasis transmission settings, as long as the distribution of focal points can be inferred. Recall also our previously stated point that, for STH, our house-hold distribution itself could potentially be a tracer for many of the focal points of infection [32, 33].

A more important effect, which is indirectly related to the variability in spatial scale of focal point clusters, is how variable cluster sizes might induce internal variability in the exposure of each individual to new infections. In models of helminth transmission [37], it is most common to account for this variability in exposure through assigning each individual an additional predisposition factor λ_i which is drawn from a gamma distribution, Gamma(λ_i; k, k), where k is the ‘aggregation parameter’ which modifies the variance-to-mean ratio of the distribution of worms within hosts. Due to the epidemiological processes that it aims to capture — and assuming that the basic reproduction number R₀, which accounts for the transmission intensity, remains unchanged — this aggregation parameter will likely vary to capture a change in exposure according to the spatial scale r_Λ chosen for the reservoir of infection. One simple model to capture this variability in exposure might be to consider the value of k itself to be drawn from a gamma distribution, i.e., p(k) = Gamma(k; a_k, b_k) from which samples are drawn across the map for the individuals closest to each focal point of infection. When binding such focal points together to get the global behaviour of the reservoir of infection, one would then have a mixture model of separate behaviours to consider in the distribution of worms within hosts, affecting the summation term of Eq. (16). Since variability of exposure with spatial scaling is not the main focus of the present work, we shall not investigate this any further here. However, due to its obvious importance to the dynamical behaviour of the reservoir by association, it will be an important component to be cognisent of in future work.

In Ref. [14] it was shown that, to good approximation for moderate-to-large population numbers, the worm burden distributions of individuals p(w_i, t) evolve according to the following reservoir birth-death process master equation The solution to this equation is a Poisson distribution p(w_i, t) = Pois[w_i; I_i(t)] with time-dependent intensity in which we have inserted the solution to Eq. (16). Summing over an ensemble of these individuals to get an overall distribution of worms within hosts, many epidemiogical variables, e.g., the prevalence of infection or mean parasite burden of hosts associated to the reservoir can be calculated (up to specifying the additional gamma-distributed predisposition to infection λ_i for each individual) while maintaining the finite population variance neglected by deterministic disease models. Therefore, in order to understand how these variables which are associated with each reservoir of infection are affected by spatial movements of human hosts, it is essential to correctly assess how Λ(t) is modified by inward migration.

Due to the movement patterns we have discussed, there may be individuals who contribute to one of the focal points of the reservoir but come from a region further away such that they are not counted in S _Λ. There may be individuals who are counted in S _Λ but instead contribute to a focal point outside of their own associated cluster, e.g., on a particular day. Eq. (16) will hence be perturbed by these movements in the following way where are also compound Poisson processes which sum over the amount of infectious material generated by these migrating individuals that is either net entering (+) or leaving (−) the reservoir focal points.

Note here that Eqs. (16) and (19) currently assume that N_Λ is static over the time period of interest, which may not be the case if one considers changes that take place over the course of, e.g., a year. This is because long-term effects can cause the population number to vary, including periodic seasonal migrant labour or family visits at particular times of year. Such long-term changes are likely to not have too great a dynamical significance on the terms of Eqs. (16) and (19), however, since these terms are associated to variations in the reservoir of infection which typically occur on much shorter timescales (1/d_res is typically on the order of weeks or months depending on helminth species).

Eq. (19) implicitly describes a network of spatially coarsegrained reservoirs of infection, each of which varying in spatial extent according to a combination of the distribution of their infectious foci, as well as the spatial coarse-graining scale choice r_Λ. This network of coarse-grained reservoir ‘nodes’ is connected by migratory ‘pulses’, continually affecting each reservoir’s temporal stability (which we shall discuss later on) as well as potentially non-negligible statistical cross-correlations in local epidemiological indicators (e.g., prevalence or intensity of infection) with respect to other reservoirs. Note also that the choice of r_Λ can therefore control the strength of these migratory interactions between reservoirs in the network. In Fig. 6 we have ilustrated this point explicitly on a digram with real map data based on buildings in central Malawi. The black (red) arrows on the map and zoomed illustrations show migratory movements of a radial distance below (above or equal to) r_Λ, which has been arbitrarily chosen for clear visual effect.

• Download figure
• Open in new tab

Figure 6.

A diagram of real map data (based on buildings in central Malawi) and a zoomed illustration which indicate the coarse-graining procedure of Eq. (21) — in which the black arrows are ‘removed’. Arrows depict individual movements from households (grey hollow dots on the right hand side zoomed illustration) to focal points of infection (black hollow dots on the right hand side zoomed illustration) where black arrows on both the map and its zoomed counterpart correspond to jumps over a radial distance r < r_Λ (where r_Λ is the spatial coarse-graining scale) and red arrows correspond to jumps over a radial distance r≥ r_Λ. On the left hand side zoomed illustration we see the emergence of the spatially coarse-grained reservoir network ‘nodes’ (nodes are neither drawn to scale nor have geometrically-accurate reservoir spatial shapes) connected by red migratory ‘links’ in time.

Let us now define T_Λ as the average rate of individuals travel-ling greater than a distance of r_Λ to a focal point of the reservoir of infection and contributing material to it. Following Ref. [14], this translates into the following decomposition of the positive migratory process in Eq. (19) where is the contribution to the infectious material quantity of a focal point by an individual: a random variable with a distribution which varies between disease, but likely follows a negative binomial character in some cases. This distribution may be specifically derived by transforming the travelling individual’s worm burden w_i (e.g., drawn from a Poisson distribution and Eq. (18) or an individual-based simulation) into an expected count of fertilised eggs or larvae using a helminth-specific mating function [36]. Note that this compound Poisson process is exact since it tracks the contribution entering the reservoir directly. If one wishes to temporally coarse-grain over the timescale 1/d_res and see the effect on the mean worm burden deterministic ODE model (as is done in Ref. [13]), an additional non-Markovian component arises to account for the reservoir pulse decay.

When an individual contributes to an infectious focal point which is further than r_Λ away from their household, they do not just make an additional contribution to the latter’s reservoir of infection, but they also remove their contribution from one of the focal points which is within r_Λ of their household (their local reservoir). In this sense, each random pulse within the process of one reservoir is exactly correlated to a pulse within the process of another’s. So, for model completeness, it is sufficient to define only the process (attributing an in-In tended location of travel) for each reservoir of infection defined over a map. This is tantamount to defining ‘links’ between the coarse-grained nodes of the network illustrated by Fig. 6.

Let us now consider the influence that the coarse-graining procedure has on Eq. (20). By analogy to physical theories in nature, note that the fundamental dynamical theory should not change — e.g., individual contributions to each reservoir should not ‘disappear’ — simply because of an arbitrary choice for r_Λ, and so one must define a procedure whereby the effective dynamics in Eq. (20) are always describing the same fundamental phenomenon, for any choice of r_Λ. An obvious way to do this is to consider the variation in probability mass associated to individual movement rates as r_Λ is varied. Based on the arguments of Sec. 2, we hence point out that the following scaling should exist for T_Λ which accounts for this difference in probability mass arising from a rescaling of distances where we have defined T_µ as the average rate of individuals travelling greater than a distance of r_µto a focal point of the reservoir of infection and contributing material to it.

In order to derive Eq. (21) we have assumed that all individuals jump once per day, so that the marginalised jump PDF p(r; α, β, r_µ) of Sec. 2 may be used directly. By integration of Eq. (21), one generally obtains where ₂F₁(·, ·; ·; ·) is the (Gauss) hypergeometric function. Fig. 7 we use the solution given by Eq. (22) to compute the quantity T_Λ/T_µ as a function of r_Λ /r_µ for a range of parameter choices.

• Download figure
• Open in new tab

Figure 7.

The spatially coarse-grained average daily pulse rate T_Λ into a focal point of infection as a fraction of its value T_µ from individuals arriving a distance of r_µ or greater away. This value is plotted as a function of the radial coarse-graining scale ratio r_Λ /r_µ used. The relationship is given by Eq. (22) for a range of α and β power-law parameters.

The variation of T_Λ/T_µ as a function of r_Λ provides a mapping between the average pulse rate for Eq. (20) that is associated to each focal point within the reservoir of infection (whose pulse amplitudes can be easily summed over to get a total contribution to and from the reservoir) and the distribution of households at different spatial coarse-graining scales. The value of r_Λ may hence be fixed to identify movement of individuals between structures at different scales, e.g., buildings, villages, towns, cities, etc. We may visualise this coarse-graining of movement distances by referring back to Fig. 6 — there is an emergence of black clusters of connected locations which are connected by red arrows at longer distances. By then clustering the focal points on this map (note this is a mock illustration and these are not real focal points) according to whether or not they are within a radial distance r_Λ of one another, we can then build coarse-grained reservoirs of infection associated to the scale r_Λ.

It is important to note that infectious reservoir clusters may themselves have a significant spatial extent due to the particular geometry of the nearby focal points, and so T_Λ should be considered to apply to each and every focal point (or small clusters of focal points) within the reservoir individually. Hence, for reservoirs of infection which correspond to clusters of focal points which, when bound together, have a radial spatial extent that is much smaller than r_Λ then the results of this section may be used directly for the entire reservoir. In situations where this is not the case — e.g., long and thin geometric relationships binding focal points around the edge of a lake — then our results may still be applied to each focal point (or small cluster of focal points) in turn. The application of our more general formalism is then subject to change depending on the specific relationships between focal points and hence, indirectly, on the species of helminth.

4.2. Reservoir Wiener process approximation

Reconsidering Eq. (16), note that the following Gaussian sample mean approximation applies for the terms summing over all (non-migratory) infectious material input into the reservoir where W_t is a Wiener process, ℳ _res(t) is a time-dependent mean and 𝒱 _res(t) is a time-dependent variance for the reservoir inputs. Such an approximation is motivated by the central limit theorem and will be most accurate in the limit of large population number N_Λ. By inserting Eq. (23) into Eq. (19), one finds a drift-diffusion which satisfies the following Gaussian distribution for Λ(t) Calculating an exact analytic form for 𝒱_res(t) in the case of each helminth species is still an open research question since this is very likely to be quite complex to model. For some analytic insight, however, in Ref. [13] it was shown that, with a negative binomial distribution of worms within hosts, the distribution over typical egg outputs from hosts for hookworm is well-approximated by another negative binomial (with different mean and variance). Under these conditions, one may infer, from the sum of negative binomial variances, that where k_res(t) is the aggregation parameter for the reservoir input negative binomial, which can be time-dependent before the system relaxes to the endemic steady state configuration.

Eq. (24) provides a significant improvement in computational efficiency over a full stochastic simulation. In particular, for endemic steady-state regions where the constant values of 𝒱_res(t) and ℳ _res(t) can be computed quickly for each parameter configuration, it could be used to generate a simulation likelihood for statistical inference over an entire map affected by helminth transmission and migrating infected individuals.

Due to migratory inputs from (and outputs to) other reservoirs, applying Eq. (24) consistently to each node of our reservoir network will require the inclusion of spatial covariances (or higher-order statistics) between samples. The strength of these spatial covariances can, in principle, be computed by including the migratory terms of Eq. (19) into the approximation above. We can see this by recalling that the rate of these migratory compound Poisson processes (see Eq. (20)) is fundamentally connected to the spatial distance scales via Eq. (22). This ‘fundamentals-based’ method of spatial epidemiological inference for helminth transmission would be distinct from other methods since it would not only include important effects which are inherent to a simulation, such as finite population variance of the underlying stochastic process, but also a theoretical understanding of how the effective dynamical description itself changes with spatial coarse-graining into separate EUs. Note that the ‘scale-invariance’ of the theoretical model we have suggested here can therefore, in principle, be inferred with data collected at one EU scale and then interpolated/extrapolated for model predictions at another.

Although it is beyond the scope of this paper, in future work, it will be of interest to compare the results obtained from other NTD spatial epidemiological models — which typically employ deterministic models a particular spatial scale, e.g., Refs. [9, 10] — with the spatial correlations induced in the diagnostic output from our suggested approach. Such a comparison would help quantify the importance of finite population variance and EU choice on reservoir temporal stability. Other NTD spatial modelling approaches one might consider a comparison to are spatial GRF-based models [6] — these are to be constrasted with the temporally Gaussian, but spatially complex, power-law-like, reservoir network we have proposed in this section. In addition, it would be interesting to explore how the approximate temporal Gaussianity of the model we have proposed in Eq. (24) changes when considering: dynamics further from endemic steady states; configurations closer to trans-mission breakpoints; or systems with significantly smaller population sizes [14].

4.3 The critical spatial scale

In order to evaluate the success of control measures through the calculation of epidemiological observables, such as the prevalence or intensity of infection, it is common for an EU to correspond to a particular spatial scale for the reservoir of infection — see, e.g., [4, 38, 39]. In Sec. 1 we discussed the possibility of defining EUs such that the effect that human movement has on their implied reservoirs of infection was minimal. Note that this is not just of importance to the temporal stability of the reservoirs, but also advantageous because it nullifies any statistical inference bias that may arise through strong migration effects. Such biases may occur when there is a large migration pulse rate of infectious material into a given region where the local intensity and prevalence are increased, but the statistical models using these diagnostics for inference do not include migration rates in their description.

It was shown in Refs. [13, 14] that if the average rate at which reservoir pulses occur exceeds the average death rate of the infectious material in the reservoir, d_res, then the effect of infected human migration on the transmission dynamics becomes particularly strong. In particular, the authors of Ref. [14] demonstrate that, for a range of finite population sizes, migration rates at or above this critical rate into a reservoir of infection can restart transmission in regions where a mass treatment programme has already achieved transmission elimination. In addition, it was shown that this rate of migration can stabilise the stronger ‘fade-out’ effects exhibited near transmission breakpoints. By combining the concept of a critical rate with our proposed spatial model for migration between reservoirs, one may derive a new and important critical spatial scale around the focal points of infectious reservoirs which depends on: the type of helminth, the small-scale geometry of locations and the specific patterns of human movement. This spatial scale can be used to define the geographical size of EUs.

If the spatial scale of an EU is chosen which corresponds to a region smaller than the critical scale around each of its focal points of infection, the average daily rate of pulses will always be larger in magnitude than the value , calculated using the critical scale. This is because the curves shown in Fig. 7 are always decreasing with increasing radial scale. This implies that when regions around focal points are defined at this critical scale or smaller for an EU, pulses of the form given by Eq. (19) cannot be safely neglected and should become important to take into account when modelling the transmission dynamics to correctly assess reservoir stability and possible observation biases. This means that if EUs are at this critical scale of smaller NTD programmes need to have measures in place to mitigate the impact of human movement on programme targets. For example, individuals who regularly move beyond the boundaries of EUs should be specifically targeted for treatment. When defining an EU in the opposite case, however, effects from migration between reservoirs of infection may be safely neglected.

Using Eq. (22) the critical spatial scale at which becomes equal to the death rate of the infectious material in the reservoir per day, i.e., , has been plotted for hookworm (which has d_res ≃ 0.071 [40]) in Fig. 8 as a function of the average daily pulse rate T_µ from individuals travelling distances at or above r_µ. In this plot we note that nearly all parameter combinations indicate a sharp decline in the critical scale ratio if T_µ is found to be below d_res — where, in particular, for values of α = 2 the scale ratio appears to fall extremely sharply and hence one can no longer find a finite critical scale below which migration becomes important.

• Download figure
• Open in new tab

Figure 8.

The critical spatial scale as a ratio at which , for hook-worm — where d_res ≃0.071 [40] — plotted as a ratio of r_µ. The value of this scale is shown against the value of the average daily pulse rate T_µ at or above r_µ. We have used Eq. (22) to generate this relationship with a (bisection) root-finding algorithm.

In contrast, for values T_µ > d_res of increasing orders of magnitude, one infers from Fig. 8 that the value of α is nearly irrelevant (which follows from the construction of T_µ) and, particularly for β = 2, 3, the increase in the critical scale ratio is very gradual. Such a relationship is also consistent with our expectation as the decline in the average pulse rate from people travelling from distances r_Λ»r_µ is particularly sharp for β = 2, 3 — see, e.g., Fig. 7 — and so one requires a significant increase in the amplitude of T_µ to achieve a significant change in the critical scale.

We acknowledge that there are important caveats to the relationship above (and the one shown in Fig. 7) which arise from the assumptions made in obtaining Eq. (21), i.e., the statistical homogeneity and isotropy of people living in the surrounding households who all jump only once per movement event/day directly to an infectious focal point. By relaxing the assumptions made in obtaining Eq. (21), the rate of pulses will vary between different reservoirs of infection. In addition to these modifying effects, the inclusion of multiple unidirectional jumps with different initial preferred directions may no longer be isotropically configured when viewed collectively over the course of the day — which would be the case when multiple individuals all converge to a globally preferred location such as a place of work or a school. We propose to consider such modifications in future work on specific case studies.

Perhaps an even more important model extension to consider was noted in Sec. 4.1, where the potential for the aggregation parameter k (for the distribution of worms within hosts) was hypothesised to also vary with spatial scale due to differences in infectious exposure. Although including this potential effect will not directly influence our conclusions with regard to migration, it is possible that another critical spatial scale exists which may be found by considering the value of k required to make endemic steady state stability possible. We leave the investigation into the possibility of this additional critical scale to future work.