Spatial close-kin mark-recapture methods to estimate dispersal parameters and barrier strength for mosquitoes

Spatial close-kin mark-recapture methods to estimate dispersal parameters and barrier strength for mosquitoes | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Spatial close-kin mark-recapture methods to estimate dispersal parameters and barrier strength for mosquitoes View ORCID Profile John M. Marshall , Shuyi Yang , Jared B. Bennett , Igor Filipović , Gordana Rašić doi: https://doi.org/10.1101/2025.05.11.653364 John M. Marshall 1 Divisions of Biostatistics and Epidemiology, School of Public Health, University of California , Berkeley, CA, USA 2 Innovative Genomics Institute, University of California , Berkeley, CA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for John M. Marshall For correspondence: john.marshall{at}berkeley.edu Shuyi Yang 1 Divisions of Biostatistics and Epidemiology, School of Public Health, University of California , Berkeley, CA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Jared B. Bennett 1 Divisions of Biostatistics and Epidemiology, School of Public Health, University of California , Berkeley, CA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Igor Filipović 3 Mosquito Genomics, QIMR Berghofer Medical Research Institute , Brisbane, Australia Find this author on Google Scholar Find this author on PubMed Search for this author on this site Gordana Rašić 3 Mosquito Genomics, QIMR Berghofer Medical Research Institute , Brisbane, Australia Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF Abstract Close-kin mark-recapture (CKMR) methods have recently been used to infer demographic parameters for several aquatic and terrestrial species. For mosquitoes, the spatial distribution of close-kin pairs has been used to estimate mean dispersal distance, of relevance to vector-borne disease transmission and genetic biocontrol strategies. Close-kin methods have advantages over traditional mark-release-recapture (MRR) methods as the mark is genetic, removing the need for physical marking and recapturing that may interfere with movement behavior. Here, we extend CKMR methods to accommodate spatial structure alongside life history for mosquitoes and comparable insects. We derive kinship probabilities for parent-offspring and full-sibling pairs in a spatial context, where an individual in each pair may be a larva or adult. Using the dengue vector Aedes aegypti as a case study, we use an individual-based model of mosquito life history to test the effectiveness of this approach at estimating parameters such as mean dispersal distance, daily staying probability, and the strength of a barrier to movement. Considering a simulated population of 9,025 adult mosquitoes arranged on a 19-by-19 grid, we find the CKMR approach provides unbiased and precise estimates of mean dispersal distance given a total of 2,500 adult females sampled over a three-month period using 25 traps evenly spread throughout the landscape. The CKMR approach is also able to estimate parameters of more complex dispersal kernels, such as the daily staying probability of a zero-inflated exponential kernel, or the strength of a barrier to movement, provided the magnitude of these parameters is greater than 0.5. These results suggest that CKMR provides an insightful characterization of mosquito dispersal that is complementary to conventional MRR methods. Author summary Close-kin mark-recapture (CKMR) is a genetic analogue of mark-release-recapture (MRR) in which the frequency of marked individuals in a sample is used to infer demographic parameters such as census population size and mean dispersal distance. These methods have been widely applied to aquatic species; however their application to mosquitoes is yet to be rigorously explored. Previous theoretical work demonstrated the potential for CKMR to infer parameters such as population size and mortality rate for randomly-mixing mosquito populations, and close-kin-based methods have been used to infer movement patterns for Aedes aegypti mosquitoes in Singapore and Malaysia. Here, we use simulations to explore the potential for formal CKMR methods to characterize mosquito dispersal patterns. We find that formal CKMR methods are able to accurately estimate mean dispersal distance, and to estimate additional parameters, such as the strength of a landscape barrier and the probability that a mosquito remains within its population node each day. CKMR and other close-kin-based methods provide insights into mosquito dispersal complementary to commonly-used alternatives such as MRR, as they capture displacement across several generations and are not compromised by the marking process. 1. Introduction Malaria, dengue, chikungunya and other mosquito-borne diseases continue to pose a major burden throughout much of the world [ 1 , 2 ]. Novel biological and genetics-based interventions, such as releases of mosquitoes infected with Wolbachia or engineered with gene drives, offer much promise to complement traditional control tools such as insecticide-treated nets, vaccines and antimalarial drugs. A common feature of these novel tools is the need for a detailed understanding of mosquito movement in order to design effective field trials and interventions, and to address biosafety concerns. In a recent randomized controlled trial (RCT) of Wolbachia -based population replacement of the dengue vector, Aedes aegypti , in Yogyakarta, Indonesia, Wolbachia was observed to spread significantly from intervention to control areas within one year of release [ 3 ]. This highlights the importance of quantifying mosquito movement to determine optimal spatial units for vector control RCTs. Predicting intentional geographic spread of self-propagating interventions such as Wolbachia and gene drive is also crucial, as is assessing the potential for confinement and logistics of reversibility during a trial [ 4 ]. A handful of methods are available to characterize mosquito movement patterns. The most direct of these is mark-release-recapture (MRR), hundreds of which studies have been conducted for Ae. aegypti and the malaria vector Anopheles gambiae in recent decades [ 5 ]. In MRR, a portion of a population is captured, marked and released, and subsequent collections are checked for recaptures. The fraction of recaptures over time can be used to infer population size and daily mortality, while times and distances between release and recapture events can be used to infer dispersal patterns. A major shortcoming of MRR is that inferred dispersal patterns may be modified by the process of marking and capturing, and for mosquitoes specifically, releasing females may increase the risk of local disease transmission. Several genetic methods are available to characterize mosquito movement on a larger spatial scale, effectively estimating dispersal averaged over several generations. Wright’s fixation index, F ST , can be calculated using genetic markers such as single nucleotide polymorphisms or microsatellites [ 6 , 7 ], and population assignment tests can be used to infer movement between well-structured populations at a scale beyond the mean dispersal range of a species [ 8 ]. Here, we theoretically explore the role that close-kin mark-recapture (CKMR) methods could play in characterizing the dispersal patterns of Ae. aegypti mosquito vectors. In CKMR, the detection of a close-kin pair (parent-offspring, siblings, etc.) in a sample is analogous to the recapture of a marked individual in the MRR method [ 9 ]. Detection of several close-kin pairs separated by a given distance provides information about movement that occurred over a small number of generations, thus informing dispersal patterns on an intermediate spatial scale, between that of MRR methods and most other genetic approaches. Advantages of CKMR methods stem from the mark being a genetically-inferred close-kin relationship, removing the need for physical marking and recapturing. Two recent studies have used the spatial distribution of close-kin pairs to characterize Ae. aegypti movement patterns [ 10 , 11 ], both using approaches inspired by the CKMR formalism [ 9 ]. Jasper et al . [ 10 ] estimate relatedness across three orders of kinship and estimate life stage-specific dispersal by considering possible life history events between each kinship pair. Filipović et al . [ 11 ] consider up to three orders of kinship and use coordinates of close-kin pairs to infer the distribution of distances between the birth and ovipositing sites of breeding females, and fit these distances to a variety of dispersal kernels. Here, we apply formal CKMR methods with spatial structure to estimate movement parameters of Ae. aegypti . Formal CKMR methods are based on explicit calculations of kinship probabilities, and the likelihood of observing a given number and category of close-kin pairs across space and time [ 9 ]. While the availability of spatial kinship information has been discussed as having the potential to inform dispersal parameters within a formal CKMR framework [ 9 , 12 , 13 ], very few studies have explored this through simulation [ 14 , 15 ]. Our approach builds upon CKMR methods specific to the life history of mosquitoes, previously used to estimate demographic parameters such as census population size, and adult and larval mortality rates [ 16 ]. As in that study, we use an in silico model of mosquito life history, this time with spatial structure, to generate kinship data and validate our inference methods. Using this approach, we determine optimal sampling schemes (sampled life stages, sample size and spatial distribution of traps) to accurately and efficiently estimate dispersal parameters, including mean dispersal distance, daily staying probability, and the strength of a barrier to movement. 2. Methods 2.1. Mosquito life history and spatial structure As per our previous mosquito CKMR study [ 16 ], we base our analysis on a discrete-time version of the lumped age-class model [ 17 , 18 ], applied to mosquitoes [ 19 ] ( Figure 1A ). This model considers discrete life history stages - egg (E), larva (L), pupa (P) and adult (A) - with sub-adult stages having defined durations - T E , T L and T P for eggs, larvae and pupae, respectively. We use a daily time-step, since mosquito samples tend to be recorded by day, and this is adequate to model the organism’s population dynamics [ 20 ]. Daily mortality rates vary according to life stage - µ E , µ L , µ P and µ A for eggs, larvae, pupae and adults, respectively - and density-dependent mortality occurs at the larval stage. Sex is modeled at the adult stage - half of pupae emerge as females (F), and the other half as males (M). Females mate once upon emergence, and retain the genetic material from that mating event for the remainder of their lives. Males mate at a rate equal to the female emergence rate which, for a population at equilibrium, is equal to the female mortality rate, µ A . Females lay eggs at a rate, β , which is assumed to be independent of age. Download figure Open in new tab Fig 1. Mosquito life history and spatial structure. In the lumped age-class model, mosquitoes are divided into four life stages: egg, larva, pupa and adult (A) . The durations of the sub-adult stages are T E , T L and T P for eggs, larvae and pupae, respectively. Sex is modeled at the adult stage, with half of pupae developing into females and half developing into males. Daily mortality rates vary by life stage - µ E , µ L , µ P and µ A for eggs, larvae, pupae and adults, respectively. Density-dependent mortality occurs at the larval stage and is a function of the total number of larvae, N L . Females mate once upon emergence, and retain the genetic material from that mating event for the remainder of their lives. Males mate at a rate equal to the female emergence rate. Females lay eggs at a rate, β . In the spatial extension of the lumped age-class model, mosquito populations are distributed in space, with movement between them defined by an exponential (solid line) or zero-inflated exponential dispersal kernel (dashed lines) (B) . The daily probability of remaining in the same population, p 0 , is varied while preserving the mean dispersal distance. This value is trimmed from the plot, but specified in the key. Mosquito populations are distributed according to a 19-by-19 grid of households (circles), with mosquito traps distributed in select households (black circles) according to the sampling scheme (C) . In some simulations and analyses, a barrier to movement is included (solid line) (D) . In extending the lumped age-class model to space, we consider a spatial distribution of mosquito populations, with movement between them defined by a dispersal kernel ( Figure 1B-D ). Discrete populations in the resulting metapopulation are considered to be randomly mixing populations to which the lumped age-class model applies. The resolution of the individual populations (in terms of size) should be chosen according to the dispersal properties of the species being considered. For Ae. aegypti , populations on the scale of households may be appropriate, as this species tends to disperse locally, often remaining in the same household for the duration of their lifespan [ 21 ]. For An. gambiae , dispersal occurs over larger distances and villages may be an appropriate population unit [ 6 ]. By default in this paper, daily per-capita movement probabilities between populations are derived from an exponential dispersal kernel ( Figure 1B ). For populations i and j a distance d ij apart, the rate of movement between them is: Here, 1/ λ d represents the mean daily dispersal distance, conditional upon movement, and n represents the number of populations in the landscape. For a given origin, i , the dispersal kernel entries, m ( i ,), sum to 1. Computing m ij for all combinations of origins and destinations produces the movement matrix, M . We also consider a zero-inflated exponential kernel ( Figure 1B ) which includes an additional parameter, p 0 , to represent the daily probability that a mosquito remains in the same population. For this kernel, we have: Here, 1 /λ c represents the mean daily dispersal distance, conditional upon movement, and the diagonal elements of the movement matrix, m ii , equal p 0 . Default life history, demographic and dispersal parameters for Ae. aegypti are listed in Table 1 . Given the difficulty of measuring juvenile stage mortality rates in the wild, these are chosen for consistency with observed population growth rates in the absence of density-dependence (see Sharma et al . [ 16 ] for formulae and derivations). Larval mortality increases with larval density and, according to the lumped age-class model, reaches a set value when the population is at equilibrium. Although mosquito populations vary seasonally, we assume a constant adult population size, N A , for this CKMR analysis, and restrict sampling to a period of three months, corresponding to a season. Minor population size fluctuations occur in the simulation model due to sampling and stochasticity. For dispersal, we consider movement within a 19-by-19 grid of households, based on a suburban setting such as Cairns in Queensland, Australia [ 21 ]. Mosquito traps are distributed in select households according to a specified sampling scheme ( Figure 1C ). In some simulations and analyses, a barrier to movement is included, which could represent a road, freeway or open park space ( Figure 1D ). View this table: View inline View popup Download powerpoint Table 1. Demographic, life history and dispersal parameters for Aedes aegypti mosquitoes. 2.2. Kinship probabilities Following the CKMR methodology of Bravington et al . [ 9 ] and its application to mosquito populations by Sharma et al . [ 16 ], we now derive spatial kinship probabilities for mother-offspring and full-sibling pairs based on the lumped age-class mosquito life history model. Each kinship probability is calculated as the reproductive output consistent with that relationship divided by the total reproductive output of all adult females in that population. In each case, we consider two individuals (adult or larva) sampled at known locations, x 1 and x 2 , and times, t 1 and t 2 , with probability symbols and references to equations listed in Table 2 . Note that mosquito sampling is lethal. View this table: View inline View popup Download powerpoint Table 2. Kinship categories, sampled life stages, sampling times, locations, and probability symbols used in spatial close-kin mark-recapture analysis. 2.2.1. Mother-offspring Let us begin with the simplest possible kinship probability, P MOL ( x 2 , t 2 | x 1 , t 1 ), which represents the probability that, given an adult female sampled at location x 1 on day t 1 , a larva sampled at location x 2 on day t 2 is her offspring. This can be expressed as the relative larval reproductive output at location x 2 on day t 2 of an adult female sampled at location x 1 on day t 1 : Here, E MOL ( x 2 , t 2 | x 1 , t 1 ) represents the expected number of surviving larval offspring at location x 2 on day t 2 from an adult female sampled at location x 1 on day t 1 , and E L ( x 2 ) represents the expected number of surviving larval offspring at location x 2 from all adult females in the population at times consistent with the time of larval sampling. Note that, since we are assuming a constant population size, E L ( x 2 ) is independent of time and is given by: Here, N F ( x 2 ) represents the equilibrium adult female population size, and y 2 represents the day of egg-laying. Considering day 0 as the reference day (in place of t 2 ), the egg must have been laid between days (0 − T E − ( T L − 1)) and (0 − T E ). Equation 4 therefore represents the expected number of offspring laid by all adult females at location x 2 that survive the egg and larva stages up to the time of sampling (day 0). E MOL ( x 2 , t 2 | x 1 , t 1 ), on the other hand, is specific to the sampled adult female, the location, x 2 , and the day of larval sampling, t 2 . This is given by: Here, the day of egg-laying, y 2 , is summed over days ( t 2 − T E − ( T L − 1)) through ( t 2 − T E ), for consistency with the larva being present on the day of sampling ( Figure 2 ). The first term in the summation represents the probability that the adult female sampled on day t 1 is alive on the day of egg-laying ( y 2 ), and the second term (in large brackets) represents the expected surviving larval output of this adult female at location x 2 on day t 2 . This latter term is equal to the probability that the larval offspring is sampled at location x 2 on day t 2 given the mother is sampled at location x 1 on day t 1 and the egg is laid on day y 2 , ψ MOL ( x 2 , t 2 | x 1 , t 1 , y 2 ), multiplied by their daily egg production, β , multiplied by the proportion of eggs that survive the egg and larva stages from the day they were laid up to the day of sampling. An indicator function is included to limit consideration to cases where the day of egg-laying lies within the adult female’s possible lifetime - i.e., between days t 1 and ( t 1 − T A ), where T A represents the maximum possible age of an adult mosquito. Although adult lifetime is exponentially-distributed, a value of T A may be chosen that captures most of this distribution and leads to accurate parameter inference. Download figure Open in new tab Fig 2. Schematic representation of spatial mother-larval offspring kinship probability. Parameters and state variables are as defined in Table 1 and §2.1. Subscript 1 refers to the parent, and subscript 2 refers to the offspring (the perspective from which probabilities are calculated). Circles represent living individuals, squares represent sampled individuals, and colors represent their locations: blue for the sampled parent, x 1 , and purple for the sampled offspring, x 2 . Parents are sampled on day t 1 , eggs are laid on day y 2 , and offspring are sampled on day t 2 . Offspring kinship probabilities are the ratio of the expected number of surviving offspring from a given adult at location x 2 on day t 2 (shown), and the expected number of surviving offspring from all adult females for this location and day (not shown). Calculating the expected number of surviving larval offspring at location x 2 on day t 2 from an adult female requires considering days of egg-laying, y 2 , consistent with maternal ages at sampling in the range [0, T A ), and larval offspring ages at sampling in the range [0, T L ). The only movement to consider is that of the mother (orange arrow). Calculating ψ MOL ( x 2 , t 2 | x 1 , t 1 , y 2 ) requires considering a single movement type - the movement of the mother between egg-laying and sampling. Given the mother is sampled at location x 1 at time t 1 , the probability that her previous location at the time of egg-laying, y 2 , is x 2 is calculated by normalizing over all possible egg-laying locations, i.e.: In general, ρ ( x j , t j | x i , t i ) represents the probability that an adult mosquito is at location x j on day t j , given its location is x i on day t i . This is given by: This is the ( x i , x j ) entry ( x i th row and x j th column) of the daily movement matrix, M , raised to the power, (1 + t j − t i ). The power accounts for the possibility that the adult female may move on each day between egg-laying (or in later calculations, emergence) and sampling, inclusive of these two days. Extending the mother-larval offspring kinship probability to the mother-adult offspring case is described in S1 Text §2.1. Extensions to father-offspring cases are described in S1 Text §2.2. 2.2.2. Full-siblings Next, we consider the full-sibling kinship probability for adult-adult pairs, P FSAA ( x 2 , t 2 | x 1 , t 1 ), which represents the probability that, given an adult sampled at location x 1 on day t 1 , an adult sampled at location x 2 on day t 2 is their full-sibling. This can be expressed as the relative adult reproductive output at location x 2 on day t 2 of the mother of a larva sampled at location x 1 on day t 1 : Here, E FSAA ( x 2 , t 2 | x 1 , t 1 ) represents the expected number of surviving adults at location x 2 on day t 2 that are full-siblings of an adult sampled at location x 1 on day t 1 , and E A ( x 2 ) represents the expected number of surviving adult offspring at location x 2 from all adult females at times consistent with the time of adult offspring sampling. Assuming a population at equilibrium, E A ( x 2 ) is independent of time and is given by: For convenience, let us refer to the adult sampled on day t 1 as individual 1. To calculate E FSAA ( x 2 , t 2 | x 1 , t 1 ), we treat the day that egg 1 is laid, y 1 , as a latent variable and take an expectation over it: Here, the expectation over the day that egg 1 is laid, y 1 , is taken over days ( t 1 − T E − T L − T P − ( T A − 1)) through ( t 1 − T E − T L − T P ), for consistency with the day that larva 1 is sampled ( Figure 3 ). The term E FSAA ( x 2 , t 2 | x 1 , t 1 , y 1 ) represents the expected number of surviving adults at location x 2 on day t 2 that are full-siblings of adult 1, conditional upon egg 1 being laid on day y 1 . Additionally, p A ( t 1 − y 1 − T E − T L − T P ) represents the probability that egg 1 is laid on day ( t 1 − y 1 − T E − T L − T P ). In general, p A ( t ) represents the probability that a given adult in the population has age t which, following from the daily adult survival probability, (1 − µ A ), is given by: Download figure Open in new tab Fig 3. Schematic representation of spatial adult-adult full-sibling kinship probabilities. Parameters and state variables are as defined in Table 1 and §2.2. Subscript 1 refers to the reference sibling, and subscript 2 refers to the sibling from whose perspective the probabilities are calculated. Circles represent living individuals, squares represent sampled individuals, and colors represent their locations: blue for sibling 1, x 1 , purple for sibling 2, x 2 , grey for the location of egg-laying for sibling 1, and pink for the location of egg-laying for sibling 2. The location of egg-laying for the firstborn sibling (here, sibling 1) is denoted by x * , and x † denotes the egg-laying location for the other sibling. Sibling 1 is sampled on day t 1 and laid on day y 1 . Sibling 2 is sampled on day t 2 and laid on day y 2 . Sibling kinship probabilities are the ratio of the expected number of surviving siblings of a given individual at location x 2 on day t 2 (shown), and the expected number of surviving offspring from all adult females for this location and day (not shown). Calculating the expected number of surviving full-siblings at location x 2 on day t 2 requires considering days of egg-laying, y 1 and y 2 , consistent with adult ages at sampling in the range [0, T A ). There are three movements to consider: those of the mother and two adult siblings (orange arrows). Movement probabilities consider both orders of egg-laying. The term E FSAA ( x 2 , t 2 | x 1 , t 1 , y 1 ) represents the expected number of surviving adults at location x 2 on day t 2 that are full-siblings of adult 1, conditional upon egg 1 being laid on day y 1 . This is given by: Here, the day of sibling egg-laying, y 2 , is summed over days ( y 1 − ( T A − 1)) through ( y 1 + ( T A − 1)), for consistency with the mother’s potential lifespan ( Figure 3 ), and considering that either of the larval offspring may have been laid first. Since we are summing two equally-weighted scenarios regarding offspring order, we include a multiplier of 1 / 2 in the expectation. The first term within the summation then represents the probability that the mother is alive on the day of sibling egg-laying, with the absolute value, | y 2 − y 1 |, accounting for both offspring orders. The second term (in larger brackets) represents the expected adult offspring output of the mother of adult 1 at location x 2 on day t 2 . This is the same equation as for the mother-larval offspring case with three exceptions: i) daily egg production is multiplied by the proportion of eggs that survive the egg, larva, pupa and adult stages up to the day of sampling to reflect the fact that adults rather than larvae are being sampled, ii) the indicator function limits consideration to cases where the day of sibling egg-laying, y 2 , is between days ( t 2 − T E − T L − T P − ( T A − 1)) and ( t 2 − T E − T L − T P ), for consistency with an adult sibling being sampled on day t 2 , and iii) there are now three movements captured by the composite movement term, ψ FSAA ( x 2 , t 2 | x 1 , t 1 , y 1 , y 2 ), which represents the probability that adult offspring 2 is sampled at location x 2 on day t 2 given that adult offspring 1 is sampled at location x 1 on day t 1 , egg 1 is laid on day y 1 , and egg 2 is laid on day y 2 . Calculating ψ FSAA ( x 2 , t 2 | x 1 , t 1 , y 1 , y 2 ) requires considering the mother’s movement between egg-laying events, in addition to the movement of both adult offspring for both offspring orders, i.e.: Here, we take an expectation over a latent egg-laying location for the firstborn sibling, x * , and multiply the probability that the firstborn sibling is laid at location x * given that adult sibling 1 is sampled at location x 1 by the probability that adult sibling 2 is sampled at location x 2 (the former probability requires normalizing over all egg-laying locations). For both adult siblings, movement begins after development through the egg, larva and pupa life stages (i.e., after ( T E + T L + T P ) days), and the mother’s movement between egg-laying events is incorporated by effectively adding | y 2 − y 1 | days of movement either backwards or forwards in time, depending on the offspring order. Extending full-sibling kinship probabilities to other life stage pairs is relatively straightforward. We consider the case of larva-larva, larva-adult and adult-larva full-sibling pairs in S1 Text §2.3. 2.3. Pseudo-likelihood calculation The goal of this spatial CKMR analysis is to make inferences about dispersal parameters given data on the frequency, timing and location of observed close-kin pairs. Here, we calculate the likelihood of parent-offspring and full-sibling pairs in a manner that takes advantage of the nature of the kinship probabilities and the sampling process. The kinship probabilities for each pair of individuals are assumed to be independent of each other, even though they are not. For this reason the combined likelihood is referred to as a “pseudo-likelihood” [ 9 ] - an approach that has been shown to produce accurate parameter and variance estimates provided the size of each sampling event is sufficiently low relative to the total population size [ 28 , 29 ]. 2.3.1. Parent-offspring pairs Let us begin by considering the mother-larval offspring kinship probability, p MOL ( x 2 , t 2 | x 1 , t 1 ), which represents the probability that, given an adult female sampled at location x 1 on day t 1 , a given larva sampled at location x 2 on day t 2 is her offspring. Now consider n F ( x 1 , t 1 ) adult females sampled at location x 1 on day t 1 . The probability that a larva sampled at location x 2 on day t 2 has a mother amongst the n F ( x 1 , t 1 ) sampled adult females, p MOL ( x 2 , t 2 | x 1 , t 1 ), is equal to one minus the probability that none of the n F ( x 1 , t 1 ) sampled adult females are the larva’s mother, i.e.: Here, P MOL ( x 2 , t 2 | x 1 , t 1 ) is as defined in Equation 3 . Now consider n L ( x 2 , t 2 ) larvae sampled at location x 2 on day t 2 , and let k MOL ( x 2 , t 2 | x 1 , t 1 ) be the number of larvae sampled at location x 2 on day t 2 that have a mother amongst the adult females sampled at location x 1 on day t 1 . The pseudo-likelihood that k MOL ( x 2 , t 2 | x 1 , t 1 ) of the n L ( x 2 , t 2 ) larvae sampled at location x 2 on day t 2 have a mother amongst the adult females sampled at location x 1 on day t 1 follows from the binomial distribution: The full log-pseudo-likelihood for mother-larval offspring pairs, Λ MOL , follows from summing the log-pseudo-likelihood over all adult female sampling days, t 1 , over consistent larval offspring sampling days, t 2 , and over all adult female and larval sampling locations, x 1 and x 2 , respectively: Note that, for the purpose of parameter interference, we can drop the first term in the pseudo-likelihood equation, and for the purpose of efficient computation, we consider consistent adult sampling days from ( t 1 + T E + T L + T P − ( T A − 1)) through ( t 1 + T E + T L + T P + ( T A − 1)). The earliest adult sampling day (relative to t 1 ) corresponds to the case where the mother laid the offspring at the beginning of her life, was sampled at the end of her life, and the adult offspring was sampled at the beginning of its life. The latest adult sampling day (relative to t 1 ) corresponds to the case where the mother was sampled on the day they laid their offspring, and the adult offspring was sampled at the end of its life. Parent-offspring pseudo-likelihood equations for other sampled sexes and life stages follow an equivalent formulation. The main point to note is that consistent offspring sampling days are specific to the kinship and sampled life stages being considered (these can be deduced from event history diagrams like those in Figure 2 ). For adult offspring cases where x 1 = x 2 and t 1 = t 2 , the number of sampled adults, n A ( x 2 , t 2 ), is reduced by one to account for the fact that an adult cannot be its own parent. The joint log-pseudo-likelihood for parent-offspring pairs is then given by: Here, Λ MOA , Λ FOL and Λ FOA denote the log-pseudo-likelihoods for mother-adult offspring pairs, father-larval offspring pairs and father-adult offspring pairs, respectively. 2.3.2. Full-sibling pairs For full-siblings, we begin with the adult-adult full-sibling kinship probability, P FSAA ( x 2 , t 2 | x 1 , t 1 ), defined in Equation 8 , which represents the probability that, given an adult sampled at location x 1 on day t 1 , an adult sampled at location x 2 on day t 2 is their full-sibling. We consider a given adult, indexed by i and sampled at location x 1 ( i ) on day t 1 ( i ), and n A ( x 2 , t 2 ) adults sampled at location x 2 on day t 2 . Let k FSAA ( i, x 2 , t 2 ) be the number of adults sampled at location x 2 on day t 2 that are full-siblings of adult i . The pseudo-likelihood that k FSAA ( i, x 2 , t 2 ) of the n A ( x 2 , t 2 ) sampled adults at location x 2 on day t 2 are full-siblings of adult i follows from the binomial distribution: Note that, for cases where x 1 ( i ) = x 2 and t 1 ( i ) = t 2 , the number of sampled adults at location x 2 on day t 2 , n A ( x 2 , t 2 ), is reduced by one to account for the fact that an adult cannot be its own sibling. Additionally, when counting siblings, we only consider siblings with indices > i to avoid double-counting. The full log-pseudo-likelihood for adult-adult full-sibling pairs, Λ FSAA , follows from summing the log-pseudo-likelihood over all sampled adults, i , over all sampled locations, x 2 , and over consistent adult sampling days, t 2 : Consistent adult sampling days for this case are from ( t 1 ( i ) − 2( T A − 1)) through ( t 1 ( i ) + 2( T A − 1)). The earliest adult sampling day (relative to t 1 ( i )) corresponds to the case where the mother laid individual 2 at the beginning of her life and individual 1 at the end of her life, adult 1 was sampled at the end of its life, and adult 2 was sampled soon after emergence. The latest adult sampling day (relative to t 1 ( i )) corresponds to the reverse case. Full-sibling pseudo-likelihood equations for other life stage pairs follow an equivalent formulation, with consistent sampling days specific to the kinship and sampled life stages being considered (these can be deduced from event history diagrams like those in Figure 3 ). The joint log-pseudo-likelihood for full-sibling pairs is then given by: Here, Λ FSLL , Λ FSLA and Λ FSAL denote the log-pseudo-likelihoods for larva-larva, larva-adult and adult-larva full-sibling pairs, respectively. 2.3.4. Parameter inference Despite parent-offspring and full-sibling kinship probabilities not being independent, the pseudo-likelihood approach enables us to combine these likelihoods, provided the size of each sampling event is sufficiently low relative to the total population size [ 9 ]. As we will see later, our simulation studies suggest this to be the case. We therefore combine these log-pseudo-likelihoods to obtain a log-pseudo-likelihood for the entire data set: Parameter inference can then proceed by varying a subset of the dispersal and/or demographic parameters in Table 1 in order to minimize −Λ. We used the nlminb function implemented in the optimx function in R [ 30 ] to perform our optimizations. 2.4. Individual-based simulation model We used a previously-developed simulation package, mPlex [ 16 ], to model mosquito life history and to test the effectiveness of the CKMR approach at estimating mosquito dispersal and demographic parameters. The model is an individual-based adaptation of a previous model, MGDrivE [ 31 ], which is a genetic and spatial extension of thelumped age-class model applied to mosquitoes by Hancock and Godfray [ 19 ] and Deredec et al . [ 20 ] ( Figure 1A ). Our previous application of mPlex to CKMR-based inference problems considered only panmmictic populations; however, the functionality had already been included to account for spatial population structure, with mosquitoes being distributed across populations in a metapopulation [ 31 ]. Each population has an equilibrium adult population size, , and exchanges migrants with the other populations. Populations are partitioned according to discrete life stages - egg, larva, pupa and adult - with sub-adult stages having fixed durations as defined earlier. See Sharma et al . [ 16 ] for more details. 3. Results We used simulated data from the individual-based mosquito model to explore the feasibility of spatial CKMR to infer dispersal parameters for Ae. aegypti . Our simulated metapopulation consisted of a 19-by-19 grid of households ( Figure 1C ) each inhabited by 25 adult mosquitoes at equilibrium with bionomic parameters listed in Table 1 . Landscape dimensions were chosen to accommodate the trap arrangements described below, as well as a buffer width of at least three non-trap nodes along each landscape edge (e.g., Figure 1D ) to reduce boundary effects. Open questions concern the optimal sampling scheme to estimate dispersal parameters for Ae. aegypti using spatial CKMR methods, and the range of dispersal parameters that can be accurately estimated. To address these questions, we first explored logistically-feasible sampling schemes to estimate mean dispersal distance by varying: i) sampled life stage (larva or adult), ii) total sample size (1,000-3,000 sequenced individuals), and iii) the number and spacing of trap nodes (arranged in 4-by-4, 5-by-5 or 6-by-6 grids with zero, one or two population nodes separating each trap node). Based on a previous analysis to estimate mosquito demographic parameters using CKMR [ 16 ], our adult samples were of females (since mosquito traps are often tailored to this sex), sampling frequency was biweekly (as is often the case for mosquito surveillance programs [ 32 ]), and sampling duration was for three months (corresponding to a season). Our likelihood calculations were based on parent-offspring and full-sibling pairs. Half-sibling pairs could be included for individual data analyses; however, computational burden prevented us from including these for exploratory analyses. Initially, we focused on estimating mean daily dispersal distance, 1 /λ d , and for subsequent analyses, also estimated daily staying probability, p 0 , and barrier strength, δ . 3.1. Optimal sampling scheme to estimate daily dispersal distance To estimate mean daily dispersal distance, our default sampling scheme consisted of a total of 1,000 sequenced individuals sampled biweekly over a three-month period spread over a 6-by-6 grid of trap nodes with one population node separating each trap node ( Figure 1C ) (i.e., ca. 1 individual sampled twice per week from each trap node, for a total of 1,000 individuals across all trap nodes after three months of sampling). We first explored the optimal distribution of sampled life stage to estimate 1 /λ d , exploring three scenarios: all larvae, all adult females, and half larvae/half adult females. Results of 100 simulation-and-analysis replicates for each scenario are depicted in Figure 4A . These suggest that all three life stage scenarios result in adequate parameter inference, in terms of both accuracy of the median and tightness of the interquartile range (IQR), albeit consistently underestimating 1 /λ d by ca. 10-15%. Since most mosquito traps are designed to target adult females, we will proceed with modeling samples of this life stage. While larval samples also produce adequate inference, either exclusively or together with adult female samples, larvae are more difficult to sample in the environment due to breeding sites sometimes being hidden or inaccessible. The problem of larval sampling is exacerbated for CKMR studies due to the requirement that individuals be sampled independently. Therefore, if multiple larvae are collected in a larval dip from a single breeding site, only one can be used in the analysis to prevent biasing the number of collected sibling pairs upwards, thus further increasing the effort required to achieve a substantial larval sample. Download figure Open in new tab Fig 4. Sampling schemes to estimate 1/ λ d for Ae. aegypti . Violin plots depict estimates of 1/ λ d for sampling scenarios described in §3.1. The simulated metapopulation consists of a 19-by-19 grid of households each inhabited by 25 adult Ae. aegypti at equilibrium with bionomic parameters listed in Table 1 . Boxes depict median and interquartile ranges of 100 simulation-and-analysis replicates for each scenario, thin lines represent 5% and 95% quantiles, points represent outliers, and kernel density plots are superimposed. The default sampling scheme consists of a total of 1,000 individuals sampled as ca. 1 individual collected twice weekly over a three-month period for each trap node, considering a 6-by-6 array of trap nodes with one population node separating each trap node ( Figure 1C ). In panel (A) , three life stage proportions are explored: all larvae, all adult females, and half larvae/half adult females. In panel (B) , all sampled individuals are adult females, and total sample sizes of 1,000, 1,500, 2,000, 2,500 and 3,000 are explored. In panel (C) , a sample size of 2,000 is adopted, and the number and spacing of trap nodes is varied (arranged in 4-by-4, 5-by-5 or 6-by-6 grids with zero, one or two population nodes separating each trap node). In panel (D) , trap nodes are arranged in a 5-by-5 grid with two population nodes separating each trap node ( Figure 1D ), and the optimal total sample size is revisited. In panel (E) , a sample size of 2,500 adult females is adopted, and the sampled life stage is revisited. The optimal sampling scheme consists of 2,500 adult females collected biweekly over a three-month period spread over a 5-by-5 grid of trap nodes with two population nodes separating each trap node. Next, we explored the optimal sample size to estimate 1/ λ d for Ae. aegypti . We performed 100 simulation-and-analysis replicates for each of five total sample sizes - 1,000, 1,500, 2,000, 2,500 and 3,000 adult females - depicted in Figure 4B . Results suggest that estimates of 1/ λ d become more precise for larger sample sizes (as measured by the IQR); but there are diminishing returns in precision for sample sizes larger than 2,000. The median estimate of 1/ λ d remains underestimated by ca. 10-15% regardless of sample size. We therefore proceed with an optimal sample size of 2,000 adult females, collected biweekly over a three-month period. This equates to ca. 2 individuals sampled twice per week from each trap node over the three-month period. Next, we explored the optimal number and spacing of trap nodes to estimate 1/ λ d for Ae. aegypti . We performed 100 simulation-and-analysis replicates for trap nodes arranged in 4-by-4, 5-by-5 or 6-by-6 grids with zero, one or two population nodes separating each trap node (except for the case of a 6-by-6 grid of trap nodes where the simulated landscape could only accommodate zero or one population node separating each trap node). Results, depicted in Figure 4C , suggest that when traps are optimally spread out throughout a landscape, inference of mean dispersal distance is very accurate (as measured by the difference between the median and true value). Notably, estimates of 1 /λ d are unbiased for trap nodes arranged in a 5-by-5 grid with two population nodes separating each trap node, and are relatively close to the true value for a smaller number of trap nodes (e.g., arranged in a 4-by-4 grid with two population nodes separating each trap node), or for a larger number of trap nodes that are closer together (e.g., arranged in a 6-by-6 grid with one population node separating each trap node). Conversely, estimates of 1 /λ d are highly biased when traps are clustered together, as can be seen for grids of 4-by-4, 5-by-5 or 6-by-6 adjacent trap nodes. In these cases, the median estimate of 1 /λ d is ca. 30% less than the true value. This is an intuitive result, as when trap nodes are clustered together, they are more likely to capture nearby movement but less likely to capture distant movement, hence leading to underestimates of mean daily dispersal. The optimal sampling scenario where trap nodes are separated by two population nodes equates to a distance between trap nodes of 49.8 meters (since population nodes are separated by 16.6 meters in these simulations), which is approximately equal to the mean lifetime dispersal distance of Ae. aegypti mosquitoes used to parameterize this model of 45.2 meters [ 11 ]. Finally, the 5-by-5 grid of trap nodes with each trap node separated by two population nodes represents the largest number of traps that can be accommodated on the simulated landscape for this degree of separation. Finally, given the significance of the number and spacing of trap nodes for estimating 1/ λ d evident in Figure 4C , we revisited the optimal sample size and life stage distribution given a 5-by-5 grid of trap nodes with each trap node separated by two population nodes. Results for sample size again suggest that samples of 2,000 adult females produce adequate estimates of 1/ λ d ; but also suggest an improvement in precision (as measured by IQR) for a population size of 2,500 ( Figure 4D ). We therefore proceed with this total sample size, which equates to ca. 3-4 individuals sampled twice per week for each trap node over the three-month period. Regarding the optimal distribution of sampled life stage, results again suggest that all three life stage scenarios result in comparable parameter inference (as measured by the accuracy of the mean and the tightness of the IQR) ( Figure 4E ). We therefore continue with an adult female sample, as per our previous reasoning. With all of these considerations in mind, the optimal sampling scheme therefore consists of a total of 2,500 adult females collected biweekly over a three-month period spread over a 5-by-5 grid of trap nodes with two population nodes separating each trap node. 3.2. CKMR-based estimates of barrier strength and daily staying probability Given the optimal sampling scheme, we employ the flexibility of the formal CKMR approach to estimate additional parameters describing more complex dispersal patterns - the strength of a barrier to movement, δ , and the daily staying probability for a zero-inflated exponential kernel, p 0 . To begin, we consider a barrier to movement, which could represent a road, freeway or open park space for Ae. aegypti , as has been documented as being important for dispersal of this species [ 21 ]. We depict the barrier as a line through the landscape ( Figure 1D ), whereby movement to the other side of the barrier is reduced by a factor, δ , and movement on the same side of the barrier is unaltered. We explore δ values in the range [0.1, 0.9] and adopt the optimal sampling scheme determined in §3.1 - a total sample of 2,500 adult females sampled from a 5-by-5 grid of trap nodes with two population nodes separating each trap node. Results in Figure 5C suggest that estimates of δ are very accurate (as measured by the distance between the median of 100 simulations and the true value) for δ ≥ 0.5, and barrier strength estimates become more precise (as measured by IQR) for δ ≥ 0.75. Estimates of barrier strength tend to be slightly overestimated for δ ≤ 0.25; but the true value still falls within the IQR. Estimates of 1/ λ d are reasonably accurate when barrier strength is simultaneously estimated ( Figure 5A ), with the true value consistently falling within the IQR; however, the median estimate from 100 simulations falls below the true value by ca. 3-5% for barrier strengths ≥ 0.75. This is perhaps not surprising, considering that the presence of a strong barrier in a landscape reduces the mean distance traveled by an individual over their lifetime. Download figure Open in new tab Fig 5. Estimates of barrier strength and daily staying probability using spatial CKMR methods. In the first (left) analysis, violin plots depict estimates of mean daily dispersal distance, 1 /λ d (A) , and barrier strength, δ (C) , obtained using the spatial CKMR approach for the optimal sampling scheme determined in §3.1, and considering a barrier to movement as depicted in Figure 1D whereby movement to the other side of the barrier is reduced by a factor, δ , in the range [0, 0.9]. In the second (right) analysis, violin plots depict estimates of mean daily dispersal distance conditional upon movement, 1/ λ c (B) , and daily staying probability, p 0 (D) , in the range [0.1, 0.9], again obtained using the spatial CKMR approach for the optimal sampling scheme determined in §3.1, and considering a zero-inflated dispersal kernel as described in Equation 2 and depicted in Figure 1B . The simulated metapopulation consists of a 19-by-19 grid of households each inhabited by 25 adult Ae. aegypti at equilibrium with bionomic parameters listed in Table 1 . Boxes depict median and interquartile ranges of 100 simulation-and-analysis replicates for each scenario, thin lines represent 5% and 95% quantiles, points represent outliers, and kernel density plots are superimposed. True parameter values are depicted by triangles, or by a dotted line when they are consistent across all cases. Next, we consider a zero-inflated exponential dispersal kernel and estimate the daily staying probability, p 0 , alongside the mean daily dispersal distance conditional upon movement, λ c . Exploring this kernel is motivated by the fact that several mosquito species, such as Ae. aegypti , obtain most of their resources from a small area such as a household and tend to remain at this location for the majority of their lifetime [ 5 ]. The zero-inflated exponential dispersal kernel is formulated in Equation 2 , and dispersal kernels having p 0 ∈ [0.1, 0.9] are depicted in Figure 1B . For each value of p 0 , we maintain a mean daily dispersal distance of 15.3 m ( Table 1 ), which translates to a daily dispersal distance conditional upon movement, 1/ λ c , of 16.1 m for p 0 = 0.1, 17.7 m for p 0 = 0.25, 21.6 m for p 0 = 0.5, 30.6 m for p 0 = 0.75, and 48.4 m for p 0 = 0.9. Again, we assume the optimal sampling scheme determined in §3.1. Results in Figure 5B suggest that estimates of 1/ λ c are accurate and precise for all explored values of p 0 , as measured by the difference between the median of 100 simulations and the true value of 1/ λ c , and by the IQR. On the other hand, p 0 can be accurately and precisely estimated for true p 0 ≥ 0.75; but smaller true values of p 0 are overestimated by increasingly large degrees (as measured by the median of 100 simulation replicates) - a true value of p 0 of 0.5 is overestimated by ca. 0.09, a true value of 0.25 is overestimated by ca. 0.20, and a true value of 0.1 is overestimated by ca. 0.36 ( Figure 5D ). Overestimates of p 0 may be related to the calculation of movement probabilities on a grid landscape. For the 19-by-19 landscape depicted in Figure 1C , an exponential kernel without zero-inflation has diagonal entries of its transition matrix (effectively, staying probabilities) between 0.18 and 0.38, and hence staying probabilities of a zero-inflated kernel with p 0 ∈ {0.1, 0.25} may be difficult to discern. 4. Discussion We have demonstrated that the CKMR formalism with spatial structure described by Bravington et al . [ 9 ] can be used to estimate dispersal parameters of Ae. aegypti , a major vector of arboviruses such as dengue and Zika virus, as a case study. Using a spatial individual-based simulation framework [ 16 ] based on the lumped age-class model [ 17 ] applied to mosquitoes [ 19 ], we have shown that these methods accurately estimate mean daily dispersal distance, 1 /λ d , and can also estimate parameters of more complicated dispersal kernels, such as the strength of a barrier to movement, δ , and the daily staying probability, p 0 , for a zero-inflated exponential kernel. The optimal sampling scheme inferred in this study is consistent with standard Ae. aegypti field sampling protocols, provided the distribution of traps is carefully selected. As for a previous in silico mosquito CKMR analysis [ 16 ], we found that sampling adult females biweekly over a period of three months is adequate, as is commonplace for mosquito sampling [ 32 ] and consistent with the length of a season. The simulated 19-by-19 grid landscape with nodes (in this case, households) spaced 16.6 meters apart was designed to resemble a suburban setting such as that of Cairns, Australia, as a case study (the neighborhood dimensions were chosen to resemble those of the Cairns suburb, Yorkeys Knob). In this setting, the optimal explored sampling scheme consisted of a five-by-five grid of trap nodes with two household nodes separating each trap node. Of note, traps in this configuration were separated by approximately the simulated lifetime dispersal distance (i.e., a separation distance of 49.8 meters c.f. a lifetime dispersal distance of 45.2 meters [ 11 ]), with the largest number of traps that could be accommodated on the simulated landscape being optimal. This optimal trap layout and the optimal sample size of 2,500 adult females is amenable to Ae. aegypti field protocols. Dividing this total sample size across the full network of traps and a three-month sampling period equates to ca. 3-4 adult females sampled twice-weekly for each trap node. This sample size is achievable; but were it to present a challenge, then additional traps could be placed throughout the landscape. It is worth noting that smaller sample sizes would likely be required for smaller populations (our simulations assumed 25 adult mosquitoes per household [ 22 ]). This formal CKMR approach to estimating mosquito dispersal parameters is complementary to the recently-proposed methods of Filipović et al . [ 11 ] and Jasper et al . [ 10 ], with each method having its own strengths and weaknesses. The Filipović et al . method is described in full in the Methods section of [ 11 ]; but in brief, adult females are captured while ovipositing, kinship categories determined, and a mean generational displacement calculated for each close-kin pair. This calculation considers the accumulated displacement between the close-kin individuals, and the set of possible movement events that led to it. A dispersal kernel is fitted to the set of mean generational displacements for all close-kin pairs. The Jasper et al . method is described in full in the Materials and Methods section of [ 10 ]; but in brief, eggs are collected from ovitraps, kinship categories determined, and “axial standard deviations” are calculated for each kinship category. The mean dispersal distance of adult females between emergence and egg-laying is then calculated from variance formulae that incorporate the axial standard deviations of observed kinship categories (see Equations 1 - 3 of [ 10 ]). This method can also be applied to sampled adult females, as demonstrated in [ 33 ], although it performs better when applied to sampled eggs. Key strengths of the Jasper et al . [ 10 ] and Filipović et al . [ 11 ] methods are that: i) they are simpler than the formal CKMR approach, and hence require less computational investment, and ii) they accommodate second and third-degree close-kin without computational burden. This however makes a systematic comparison of the performance of the three methods difficult, as the Jasper et al . and Filipović et al . methods are ideally performed on larger landscapes that accommodate displacement accumulated over 2-3 generations, while the formal CKMR approach can be implemented on a smaller landscape, inferring mean daily dispersal from displacement accumulated over shorter time periods spanning one day through two generations. Of note, the Filipović et al . method can be applied exclusively to first-degree close-kin, making a reduced landscape amenable to that approach; however, this comes at the expense of sacrificing second and third-degree close-kin data that the method otherwise benefits from. For the formal CKMR approach, the incorporation of temporal information enables higher-resolution inference of dispersal parameters from data collected over a smaller area; however, the computational burden associated with this additional detail limits the size of a landscape that can be considered using this method. A key benefit of the formal CKMR approach is its ability to estimate parameters of more complex dispersal kernels, such as the staying probability of a zero-inflated exponential kernel, and of more complex landscapes, such as the strength of a barrier to dispersal. Exponential and zero-inflated exponential kernels were explored in this analysis; but any number of dispersal kernels could be explored, provided available data is consistent with identifiability of their parameters. A related benefit of the CKMR approach is that it can be tailored to the life history and landscape of a specific mosquito species and location. This includes three-dimensional landscapes such as the multi-storey housing blocks analyzed in Singapore by Filipović et al . [ 11 ]. It should be noted that, with this capability comes a need to understand the local ecology of the species before parameter inference begins. E.g., when estimating dispersal parameters in this study, we incorporated a fully-specified life history model in addition to knowledge of the distribution of habitat patches and the functional form of the dispersal kernel. That said; the formal CKMR approach can also be used to estimate demographic parameters unrelated to dispersal, such as census population size, adult and larval daily mortality rates, and larval life stage duration, as shown in Sharma et al . [ 16 ]. As a preliminary exploration of the application of formal CKMR methods to estimate dispersal parameters of mosquitoes, this study has several limitations. First, the same life history and landscape model ( Figure 1 ) was used as a basis for both the population simulations and CKMR analysis. Additionally, other than the parameters being estimated, the same parameters were used in both the simulations and analysis. This represents an overly generous scenario as compared to the field, where life history is varied and complex and parameters are only approximately known. That said; this is an appropriate starting point to verify the utility of the method. Second, we have assumed perfect kinship inference throughout. Incorporating kinship uncertainty into the CKMR likelihood equations is theoretically possible [ 34 ], although this has produced little improvement in parameter inference at large computational cost when applied to data from fish species [ 28 ]. One approach to assess the importance of kinship assignment errors would be to introduce them at the simulation phase, and to then test the robustness of the methods to these errors. Here, there is an important distinction between type I (false positive) and II (false negative) error rates. Studies in fish species suggest that kinship inference must have a very low type I error rate in order for CKMR parameter inference to be informative [ 9 ]. The application of formal CKMR methods to spatial settings, while contemplated since their inception [ 9 ], has only been considered in a small number of cases for species with simpler life histories [ 14 , 15 ]. The extension of these methods to a metapopulation of insects having an egg-larva-pupa-adult life history is promising for insights this approach may provide for other species, alongside insights from simpler close-kin-based approaches [ 10 , 11 ]. Potential applications to An. gambiae are of particular interest, given the importance of this species for malaria transmission and the importance of quantifying dispersal patterns for planning vector control interventions and field trials. The increased dispersal range of this species is important to acknowledge [ 5 ], as is the potential for age-grading techniques [ 35 ] to enhance parameter inference. Several species of insect agricultural pests share a similar life history, and close-kin-based approaches such as CKMR should be explored to provide insight into their dispersal patterns. 5. Conclusions We have theoretically demonstrated the application of spatial CKMR methods to estimate dispersal parameters for mosquitoes, with Ae. aegypti as a case study. Close-kin-based methods have advantages over traditional MRR methods, as the mark is genetic, removing the need for physical marking and recapturing. Encouragingly, we find that optimal spatial CKMR sampling schemes are consistent with Ae. aegypti ecology and field studies, provided the spatial distribution of traps is carefully chosen. In our in silico case study, we found that traps distributed in a grid layout separated at a distance approximately equal to the mean lifetime dispersal distance of this species were optimal. The formal CKMR approach is complementary to two simpler close-kin-based methods at estimating mean dispersal distance, with each approach having its own strengths and weaknesses. The CKMR approach is particularly computationally intensive, restricting its ability to be applied to second and third-degree close-kin and to larger landscapes; but its ability to be tailored to a specific landscape and dispersal kernel enable it to estimate additional parameters such as the daily staying probability and strength of a barrier to movement. Close-kin-based methods such as CKMR promise to provide further insight into the dispersal patterns of other insects of epidemiological and agricultural significance. Supporting information S1 Text. Supplemental model equations . Additional equations describing mosquito dispersal dynamics, and spatial kinship probabilities for parent-offspring and full-sibling pairs that, for brevity, were not included in the manuscript. Author contributions Conceptualization: John M. Marshall, Gordana Rašić Data curation: John M. Marshall, Shuyi Yang, Jared B. Bennett, Igor Filipović, Gordana Rašić Formal analysis: John M. Marshall, Shuyi Yang, Jared B. Bennett, Igor Filipović, Gordana Rašić Funding acquisition: John M. Marshall, Gordana Rašić Investigation: John M. Marshall, Shuyi Yang, Igor Filipović, Gordana Rašić Methodology: John M. Marshall, Shuyi Yang, Jared B. Bennett, Igor Filipović, Gordana Rašić Project administration: John M. Marshall Resources: John M. Marshall, Gordana Rašić Software: John M. Marshall, Shuyi Yang, Jared B. Bennett Supervision: John M. Marshall Validation: John M. Marshall Visualization: John M. Marshall, Shuyi Yang, Gordana Rašić Writing - original draft preparation: John M. Marshall Writing - review & editing: John M. Marshall, Shuyi Yang, Jared B. Bennett, Igor Filipović, Gordana Rašić Funding statement This work was supported by a National Institutes of Health R01 grant (1R01AI143698) awarded to J.M.M. and G.R. Competing interests The authors have declared that no competing interests exist. Data availability The source code for the individual-based mosquito simulation model is available at https://github.com/GilChrist19/mPlex . Documentation, including vignettes, are included for all simulation functions. The source code for inferring parameters based on the likelihood of the kinship data is available at https://github.com/MarshallLab/CKMR . Both sets of code are available under the GPL3 License and are free for other groups to modify and extend as needed. Acknowledgments We thank Dr. Yogita Sharma, Dr. Rachel Fewster and Dr. Mark Bravington for discussions regarding the application of CKMR methods to mosquito populations. Funder Information Declared National Institutes of Health, https://ror.org/01cwqze88 , 1R01AI143698 Footnotes https://github.com/MarshallLab/CKMR References 1. ↵ Bhatt S , Gething PW , Brady OJ , Messina JP , Farlow AW , Moyes CL , Drake JM , Brownstein JS , Hoen AG , Sankoh O , Myers MF , George DB , Jaenisch T , Wint GRW , Simmons CP , Scott TW , Farrar JJ , Hay SI . The global distribution and burden of dengue . Nature . 2013 ; 496 : 504 – 507 . OpenUrl CrossRef PubMed Web of Science 2. ↵ World Health Organization . World Malaria Report 2022 . Geneva, Switzerland . 3. ↵ Utarini A , Indriani C , Ahmad RA , Tantowijoyo W , Arguni E , Ansari MR , Supriyati E , Wardana DS , Meitika Y , Ernesia I , Nurhayati I , Prabowo E , Andari B , Green BR , Hodgson L , Cutcher Z , Rancès E , Ryan PA , O’Neill SL , Dufault SM , Tanamas SK , Jewell NP , Anders KL , Simmons CP . Efficacy of Wolbachia-infected mosquito deployments for the control of dengue . N Engl J Med . 2021 ; 384 : 2177 – 2186 . OpenUrl CrossRef PubMed 4. ↵ James S , Collins FH , Welkhoff PA , Emerson C , Godfray HCJ , Gottlieb M , Greenwood B , Lindsay SW , Mbogo CM , Okumu FO , Quemada H , Savadogo M , Singh JA , Tountas KH , Touré YT. Pathway to deployment of gene drive mosquitoes as a potential biocontrol tool for elimination of malaria in sub-Saharan Africa: Recommendations of a scientific working group . Am J Top Med Hyg . 2018 ; 98 : 1 – 49 . OpenUrl 5. ↵ Guerra CA , Reiner RC , Perkins TA , Lindsay SW , Midega JT , Brady OJ , Barker CM , Reisen WK , Harrington LC , Takken W , Kitron U. A global assembly of adult female mosquito mark-release-recapture data to inform the control of mosquito-borne pathogens . Parasites & Vectors . 2014 ; 7 : 1 – 5 . OpenUrl CrossRef PubMed 6. ↵ Taylor CE , Touré YT , Carnahan J , Norris NE , Dolo G , Traoré SF , Edillo FE , Lanzaro GC . Gene flow among populations of the malaria vector, Anopheles gambiae, in Mali, West Africa . Genetics . 2001 ; 157 : 743 – 750 . OpenUrl Abstract / FREE Full Text 7. ↵ Chen H , Minakawa N , Beier J , Yan G. Population genetic structure of Anopeheles gambiae mosquitoes on Lake Victoria islands, west Kenya . Malar J . 2004 ; 3 : 48 . OpenUrl CrossRef PubMed 8. ↵ Schmidt H , Lee Y , Collier TC , Hanemaaijer MJ , Kirstein OD , Ouledi A , Muleba M , Norris DE , Slatkin M , Cornel AJ , Lanzaro GC . Transcontinental dispersal of Anopheles gambiae occurred from West African origin via serial founder events . Commun Biol . 2019 ; 2 : 473 . OpenUrl CrossRef PubMed 9. ↵ Bravington MV , Skaug HJ , Anderson EC . Close-kin mark-recapture . Statistical Science . 2016 ; 31 : 259 – 274 . OpenUrl 10. ↵ Jasper M , Schmidt TL , Ahmad NW , Sinkins SP , Hoffmann AA . A genomic approach to inferring kinship reveals limited intergenerational dispersal in the yellow fever mosquito . Molecular Ecology Resources . 2019 ; 19 : 1254 – 1264 . OpenUrl CrossRef PubMed 11. ↵ Filipovíc I , Hapuarachchi HC , Tien WP , Razak MA , Lee C , Tan CH , Devine GJ , Rašíc G. Using spatial genetics to quantify mosquito dispersal for control programs . BMC Biology . 2020 ; 18 : 1 – 5 . OpenUrl CrossRef PubMed 12. ↵ Conn PB , Bravington MV , Baylis S , Ver Hoef JM . Robustness of close-kin mark-recapture estimators to dispersal limitation and spatially varying sampling probabilities . Ecology and Evolution . 2020 ; 10 : 5558 – 5569 . OpenUrl CrossRef 13. ↵ Trenkel VM , Charrier G , Lorance P , Bravington MV . Close-kin mark-recapture abundance estimation: practical insights and lessons learned . ICES Journal of Marine Science . 2022 ; 79 : 413 – 422 . OpenUrl CrossRef 14. ↵ Akita T. Estimating contemporary migration numbers of adults based on kinship relationships in iteroparous species . Molecular Ecology Resources . 2022 ; 22 : 3006 – 3017 . OpenUrl CrossRef PubMed 15. ↵ Sévêque A , Lonsinger RC , Waits LP , Morin DJ . Spatial close-kin mark-recapture models applied to terrestrial species with continuous natal dispersal . Methods in Ecology and Evolution . 2025 ; doi: 10.1111/2041-210X.14490 . OpenUrl CrossRef 16. ↵ Sharma Y , Bennett JB , Rašíc G , Marshall JM . Close-kin mark-recapture methods to estimate demographic parameters of mosquitoes . bioRxiv . 2022 ; doi: 10.1101/2022.02.19.481126 . OpenUrl Abstract / FREE Full Text 17. ↵ Freedman HI , Strobeck C Gurney WS , Nisbet RM . The systematic formulation of delay-differential models of age or size structured populations . In: Freedman HI , Strobeck C , editors. Population Biology . Springer , Berlin ; 1983 . pp. 163 – 172 . 18. ↵ Gurney WS , Nisbet RM , Lawton JH . The systematic formulation of tractable single-species population models incorporating age structure . Journal of Animal Ecology . 1983 ; 52 : 479 – 495 . OpenUrl CrossRef 19. ↵ Hancock PA , Godfray HC . Application of the lumped age-class technique to studying the dynamics of malaria-mosquito-human interactions . Malaria Journal . 2007 ; 6 : 98 . OpenUrl CrossRef PubMed 20. ↵ Deredec A , Godfray HC , Burt A. Requirements for effective malaria control with homing endonuclease genes . Proceedings of the National Academy of Sciences . 2011 ; 108 : E874 – 880 . OpenUrl Abstract / FREE Full Text 21. ↵ Schmidt TL , Filipovíc I , Hoffmann AA , Rašíc G. Fine-scale landscape genomics helps explain the slow spatial spread of Wolbachia through the Aedes aegypti population in Cairns, Australia . Heredity . 2018 ; 120 : 386 – 395 . OpenUrl CrossRef PubMed 22. ↵ Hoffmann AA , Montgomery BL , Popovici J , Iturbe-Ormaetxe I , Johnson PH , Muzzi F , Greenfield M , Durkan M , Leong YS , Dong Y , Cook H , Axford J , Callahan AG , Kenny N , Omodei C , McGraw EA , Ryan PA , Ritchie SA , Turelli M , O’Neill SL . Successful establishment of Wolbachia in Aedes populations to suppress dengue transmission . Nature . 2011 ; 476 : 454 – 457 . OpenUrl CrossRef PubMed Web of Science 23. Focks DA , Haile DG , Daniels E , Mount GA . Dynamic life table model for Aedes aegypti (Diptera: Culicidae): analysis of the literature and model development . Journal of Medical Entomology . 1993 ; 30 : 1003 – 1017 . OpenUrl CrossRef PubMed 24. Otero M , Solari HG , Schweigmann N. A stochastic population dynamics model for Aedes aegypti : formulation and application to a city with temperate climate . Bulletin of Mathematical Biology . 2006 ; 68 : 1945 – 1974 . OpenUrl CrossRef PubMed 25. Eisen L , Monaghan AJ , Lozano-Fuentes S , Steinhoff DF , Hayden MH , Bieringer PE . The impact of temperature on the bionomics of Aedes (Stegomyia) aegypti, with special reference to the cool geographic range margins . Journal of Medical Entomology . 2014 ; 51 : 496 – 516 . OpenUrl CrossRef PubMed 26. Simoy MI , Simoy MV , Canziani GA . The effect of temperature on the population dynamics of Aedes aegypti . Ecological Modelling . 2015 ; 314 : 100 – 110 . OpenUrl CrossRef 27. Harrington LC , Edman JD , Kittayapong P , Coleman RC , Clark GG , Sithiprasasna R , et al. Dispersal of the dengue vector Aedes aegypti within and between rural communities . Am J Trop Med Hyg . 2005 ; 72 : 209 – 220 . OpenUrl Abstract / FREE Full Text 28. ↵ Bravington MV , Grewe PM , Davies CR . Absolute abundance of southern bluefin tuna estimated by close-kin mark-recapture . Nature Communications . 2016 ; 7 : 1 – 8 . OpenUrl 29. ↵ Hillary RM , Bravington MV , Patterson TA , Grewe P , Bradford R , Feutry P , et al. Genetic relatedness reveals total population size of white sharks in eastern Australia and New Zealand . Scientific Reports . 2018 ; 8 : 1 – 9 . OpenUrl PubMed 30. ↵ Nash JC , Varadhan R. Unifying optimization algorithms to aid software system users: optimx for R . Journal of Statistical Software . 2011 ; doi: 10.18637/jss.v043.i09 . OpenUrl CrossRef 31. ↵ Sánchez C. tHM , Wu SL , Bennett JB , Marshall JM . MGDrivE: A modular simulation framework for the spread of gene drives through spatially-explicit mosquito populations . Methods in Ecology and Evolution . 2020 ; 11 : 229 – 239 . OpenUrl CrossRef 32. ↵ Ndiaye EH , Diallo D , Fall G , Ba Y , Faye O , Dia I , et al. Arboviruses isolated from the Barkedji mosquito-based surveillance system, 2012-2013 . BMC Infectious Diseases . 2018 ; 18 : 642 . pmid: 30541472 OpenUrl CrossRef PubMed 33. ↵ Jasper M , Hoffmann AA , Schmidt TL . Estimating dispersal using close kin dyads: The kindisperse R package . Molecular Ecology Resources . 2021 ; 22 : 1200 – 1212 . OpenUrl PubMed 34. ↵ Skaug HJ . Allele-sharing methods for estimation of population size . Biometrics . 2001 ; 57 : 750 – 756 . OpenUrl CrossRef PubMed Web of Science 35. ↵ Johnson BJ , Hugo LE , Churcher TS , Ong OTW , Devine GJ . Mosquito age grading and vector-control programmes . Trends in Parasitology . 2020 ; 36 : 39 – 51 . OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted May 15, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Spatial close-kin mark-recapture methods to estimate dispersal parameters and barrier strength for mosquitoes Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Spatial close-kin mark-recapture methods to estimate dispersal parameters and barrier strength for mosquitoes John M. Marshall , Shuyi Yang , Jared B. Bennett , Igor Filipović , Gordana Rašić bioRxiv 2025.05.11.653364; doi: https://doi.org/10.1101/2025.05.11.653364 Share This Article: Copy Citation Tools Spatial close-kin mark-recapture methods to estimate dispersal parameters and barrier strength for mosquitoes John M. Marshall , Shuyi Yang , Jared B. Bennett , Igor Filipović , Gordana Rašić bioRxiv 2025.05.11.653364; doi: https://doi.org/10.1101/2025.05.11.653364 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Ecology Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17691) Bioengineering (13892) Bioinformatics (41937) Biophysics (21452) Cancer Biology (18588) Cell Biology (25504) Clinical Trials (138) Developmental Biology (13378) Ecology (19899) Epidemiology (2067) Evolutionary Biology (24320) Genetics (15609) Genomics (22506) Immunology (17736) Microbiology (40394) Molecular Biology (17181) Neuroscience (88605) Paleontology (666) Pathology (2832) Pharmacology and Toxicology (4824) Physiology (7641) Plant Biology (15156) Scientific Communication and Education (2045) Synthetic Biology (4294) Systems Biology (9825) Zoology (2271)