Integrating an individual-based model with approximate Bayesian computation to predict the invasion of a freshwater fish provides insights into dispersal and range expansion dynamics

Short-distance dispersal enables introduced alien species to colonise and invade local habitats following their initial introduction, but is often poorly understood for many freshwater taxa. Knowledge gaps in range expansion of alien species can be overcome using predictive approaches such as individual based models (IBMs), especially if predictions can be improved through fitting to empirical data, but this can be challenging for models having multiple parameters. We therefore estimated the parameters of a model implemented in the RangeShifter IBM platform by approximate Bayesian computation (ABC) in order to predict the further invasion of a lowland river (Great Ouse, England) by a small-bodied invasive fish (bitterling Rhodeus sericeus). Prior estimates for parameters were obtained from the literature and expert opinion. Model fitting was conducted using a time-series (1983 to 2018) of sampling data at fixed locations and revealed that for 5 of 11 model parameters, the posterior distributions differed markedly from prior assumptions. In particular, sub-adult maximum emigration probability was substantially higher in the posteriors than priors. Simulations of bitterling range expansion predicted that following detection in 1984, their early expansion involved a relatively high population growth rate that stabilised after 5 years. The pattern of bitterling patch occupancy was sigmoidal, with 20% of the catchment occupied after 20 years, increasing to 80% after 30 years. Predictions were then for 95% occupancy after 69 years. The development of this IBM thus successfully simulated the range expansion dynamics of this small-bodied invasive fish, with ABC improving the simulation precision. This combined methodology also highlighted that sub-adult dispersal was more likely to contribute to the rapid colonisation rate than expert opinion suggested. These results emphasise the importance of time-series data for refining IBM parameters generally and increasing our understanding of dispersal behaviour and range expansion dynamics specifically.

Abstract Short-distance dispersal enables introduced alien species to colonise and invade local habitats following their initial introduction, but is often poorly understood for many freshwater taxa. Knowledge gaps in range expansion of alien species can be overcome using predictive approaches such as individual based models (IBMs), especially if predictions can be improved through fitting to empirical data, but this can be challenging for models having multiple parameters. We therefore estimated the parameters of a model implemented in the RangeShifter IBM platform by approximate Bayesian computation (ABC) in order to predict the further invasion of a lowland river (Great Ouse, England) by a small-bodied invasive fish (bitterling Rhodeus sericeus). Prior estimates for parameters were obtained from the literature and expert opinion. Model fitting was conducted using a time-series (1983 to 2018) of sampling data at fixed locations and revealed that for 5 of 11 model parameters, the posterior distributions differed markedly from prior assumptions.
In particular, sub-adult maximum emigration probability was substantially higher in the posteriors than priors. Simulations of bitterling range expansion predicted that following detection in 1984, their early expansion involved a relatively high population growth rate that stabilised after 5 years. The pattern of bitterling patch occupancy was sigmoidal, with 20% of the catchment occupied after 20 years, increasing to 80% after 30 years. Predictions were then for 95% occupancy after 69 years. The development of this IBM thus successfully simulated the range expansion dynamics of this small-bodied invasive fish, with ABC improving the simulation precision. This combined methodology also highlighted that sub-adult dispersal was more likely to contribute to the rapid colonisation rate than expert opinion suggested. These results emphasise the impor-

Introduction
Biological invasions are a major aspect of global environmental change, responsible for pervasive changes to native biota and ecosystems (Simberloff et al. 2013;Gámez-Virués et al. 2015). Following the introduction of an alien species into a new range, its invasion success depends, at least in part, on its dispersal dynamics (Byers and Pringle 2006;Havel et al. 2015). In order to make more informed decisions on measures required to control and contain invasive species, managers need information on dispersal dynamics, along with establishment rates and ecological impacts (Gozlan et al. 2010a;Early et al. 2016). The extent of many invasions has proved difficult to predict due to a general lack of knowledge on dispersal dynamics and their relationship with population parameters, especially where the invader lacks data from their native range (Karakus et al. 2018), so here we demonstrate an analytical method that potentially overcomes these issues.
Predictive approaches that provide realistic representations of real-life invasions and enable scenario testing can develop understanding of the dispersal dynamics of invasive species Samson et al. 2017). A range of modelling approaches exist for predicting the dynamics of range expansions, including analytical methods such as integro-difference modelling (Gilbert et al. 2014(Gilbert et al. , 2017 and stochastic simulations, including individual-based models (IBMs). IBMs have the benefits of flexibility in model formulation, although they tend to be slower to run during formal model fitting approaches and so can be more challenging to use. However, improvements in computer performance in recent years have helped overcome some of these challenges, resulting in IBMs being increasingly applied to ecological issues (e.g. Hedger et al. 2013a, b;DeAngelis and Grimm 2014;Boyd et al. 2018), with approaches for fitting these models to empirical data now emerging (van der Vaart et al. 2018).
Applications of IBMs to invasive species have included simulations on how population control efforts affect the individual movement and population demographics of Eastern brook trout Salvelinus fontinalis (Day et al. 2018) and invasive sea lamprey Petromyzon marinus (Madenjian et al. 2003;Neeson et al. 2012). The performance of invasion IBMs can, however, be improved when empirical data are available that enable model fitting and enable the parameters that most strongly influence the predicted patterns to be identified (Phang et al. 2016), but these data are rarely available. For example, although Samson et al. (2017) investigated the spread of invasive round goby Neogobius melanostomus using the 'RangeShifter' platform , with model parameters developed from stakeholder interaction, scientific literature and inverse modelling approaches, the model could not be calibrated fully due to an absence of empirical data on their invasion. When empirical data are available, the modelling processes can utilise these to derive more robust estimates of model parameters by the application of inverse fitting techniques, such as approximate Bayesian computation (ABC; van der Vaart et al. 2015). ABC enables estimates of model parameters to be refined by combining information from empirical data (such as spatial and temporal distribution data), with prior probabilities derived from literature and/or expert knowledge (van der Vaart et al. 2015). ABC has been used to estimate model parameters within IBMs by Boyd et al. (2018), who developed a generic marine fish bioenergetics IBM for evaluating fish population dynamics. To our knowledge, an ABC process within an IBM has not been applied to modelling the dispersal dynamics of an invading species, despite the importance of estimating their dispersal and population parameters, and how these vary with time since their introduction and establishment (Alford et al. 2009).
The importance of using processes such as ABC to estimate dispersal parameters of invaders is that these parameters strongly influence the invasion process. Whilst introductions of alien species via long-range dispersal events (via transport or through passive dispersal) strongly influence the large-scale rate of range expansion of invaders (e.g. between countries and regions) (Hastings et al. 2005;Wilson et al. 2009), these events are then followed by range infilling via diffusive spread and/or active dispersal (e.g. within the country or region of the introduction) (Gozlan et al. 2010b). Here, diffusion theory suggests that there will be a symmetrical, radial expansion from the area of introduction (Skellam 1951), with the rate of expansion dependent on the interaction of the dispersal ability and reproductive rate of the population (Fraser et al. 2015). It will also be influenced by factors including spatial heterogeneity, temporal variability and biotic resistance from native species (Hastings et al. 2005). This dispersal of the alien species at the 'leading edge' of their invasion range is important for understanding the rate at which habitats in the new region are colonised (Wilson et al. 2009) and the population parameters that are driving this.
Within fluvial environments, dispersal events are also likely to be subject to directional bias, given the ease of downstream movement by individuals via passive drift (Byers and Pringle 2006). If upstream colonisation is to be achieved, then active dispersal is important, especially if obstacles are to be overcome that can impede movement (Vitule et al. 2012). This makes it especially important to understand the processes driving active dispersal at the upstream leading edge of the invasion range. Tracking the natural dispersion of invasive fishes in rivers can be difficult when anthropogenic activities occur, such as unregulated secondary release (stocking) events by anglers to increase angling opportunity, as these releases are likely to result in more rapid colonisation rates than possible by natural processes alone (Antognazza et al. 2016). However, these activities are less apparent when the invader has low recreational and socio-economic value, such as in many small-bodied alien fishes (especially if the species is rarely used as bait fish by anglers), meaning their colonisation rates are primarily due to natural dispersal alone (Davies et al. 2013;Davies and Britton 2016).
The aim of this study was to thus incorporate an ABC process into an IBM for simulating the 30 year invasion of an alien, small-bodied fish in a river basin, including identifying the population parameters that most strongly influenced their rates of dispersion, and then predicting their future range expansion. The model invader was bitterling Rhodeus sericeus and the modelled river basin was the Great Ouse in Eastern England, with the model developed on the Range-Shifter IBM platform. In the Great Ouse, bitterling has undergone a natural range expansion since the 1980s, with the species not considered to have been subject to multiple releases due to their low recreational and socio-economic value.

Study species and river system
The bitterling is a freshwater fish of the Cyprinidae family that was introduced into Britain in the 1920 s, probably for ornamental reasons (Davies et al. 2004;Damme et al. 2007). A small-bodied (\ 70 mm) littoral species, it shares many life history traits with other small-bodied invasive fishes, such as a limited lifespan (\ 4 years) and early sexual maturity (in the second or third year of life). It has no angling or aquaculture value, so is considered as rarely subject to secondary stocking events for fishery interests (Davies et al. 2004). Unlike other small-bodied invaders, however, its reproduction involves a parasitic relationship with freshwater mussels, where females lay their eggs within the mussel gills Reynolds 2002, 2003;Damme et al. 2007). The presence of eggs in gills can impact mussel performance through decreased ventilation, food intake and growth (Reichard et al. 2006). The quality of individual mussels as hosts also reduces with increased parasitism due to gill damage (Mills et al. 2005;Reichard et al. 2007;Smith 2017), suggesting some density-dependent regulation of bitterling reproduction.
The River Great Ouse rises in central England and flows in a generally north-easterly direction before entering the North Sea, and drains a catchment of approximately 8380 km 2 . In its lower reaches, the river flows through areas of low-lying land of low gradient (fenland). In these fenland areas, the river channel is characterised by anthropogenic alteration for land drainage and flood relief, including the presence of artificial drainage ditches and pumping stations (Mostert 2017). The bitterling is believed to have been introduced into its tributary, the River Cam, in the mid-1970s (Davies et al. 2004), although the reason for its release (such as whether it was accidental or intentional) is not known.

Bitterling time-series data in the River Great Ouse
The Environment Agency (the public regulatory body for inland fisheries in England) and its predecessors commenced monitoring of the fish assemblage of the Great Ouse catchment in 1984; this monitoring involves sampling up to 72 sites approximately every 3 years using a consistent seine netting methodology (Bayley and Herendeen 2000). In the initial surveys of the mid-1980s, bitterling was captured only in a series of small channels that connect to the River Cam ( Fig. 1). Surveys completed up to 2017 provide data that enables their spread throughout much of the lower catchment to be tracked. The minimum data available for each of these surveys are the site location, date of sampling, and the number of bitterling captured (Environment Agency 2018). Due to the relatively large mesh size of the seine nets used, the majority of captured bitterling were 55 to 70 mm in length, i.e. mature adults of generally 2 or 3 years old, with smaller individuals unrepresented in catches. It was this data time-series that was used to estimate key parameters of the IBM (Environment-Agency 2018).

IBM development to model bitterling range expansion
The platform RangeShifter was used to develop the IBM, as it allows the development of spatially explicit IBMs in which the three key dispersal phases of emigration, transfer and settlement are represented independently Samson et al. 2017). A customised version of RangeShifter v. 2.0 was used, which incorporated code to estimate model parameters by ABC given suitable high-level observed data. Development of the IBM required the following steps: mapping the catchment (including splitting the continuous river stretches into discrete sub-population patches), collating data on bitterling demography, dispersal and range expansion in the Great Ouse catchment, setting the prior distributions of parameters to be estimated, model fitting by ABC and finally, simulation of future spread. Fitting the ABC was examined by using the customised version of Range-Shifter, with the simulation of future bittering spread then completed in RangeShifter v 1.1 ) using the final set of parameter samples from the ABC.

Mapping the catchment
The initial catchment layers were extracted as tiles of geographic mark-up language from Ordnance Survey data (Ordnance Survey 2018) and converted into a shapefile using QGIS 2.14.20. Following removal of minor drainage channels in which fish were assumed not to be present, it was converted to a raster format at a resolution of 50 m using ArcGIS 10.3.1. However, as the channel width throughout most of the catchment was substantially less than 50 m, each river cell was assigned a quality score using the mean channel width of the features enclosed by each cell. Thus, we created a raster of habitat quality, with this quality measure based on our assumption that sections of greater river width would have larger and more heterogeneous littoral areas that provided enhanced habitat quality for bitterling.
The catchment was then divided into a series of non-overlapping contiguous patches, where a patch comprised of a set of cells that delimited an area of river. Although the Great Ouse is a continuous and linear system, the use of patches assumed that each patch delimited the range of a reasonably selfcontained sub-population, which was connected to neighbouring sub-populations by dispersal. The delimitations were made in a consistent manner, in which confluences, weirs and other anthropogenic features were used to delimit patches where they were present (Fig. 1); elsewhere, patches were based on segmentation in the Ordnance Survey data which in turn was based on stream width (n = 272, mean length = 2420 m, SD = 1270 m). The variation in patch length had a minimum effect in the model, as patches of contrasting length were well-mixed locally throughout the catchment.

Bitterling development stages and population parameters
Three stages of bitterling development were defined in the model, 'juveniles', 'sub-adults' and 'adults'. Juveniles were the fish that initially develop inside the mussel gills before emergence and were less than one year old (0 ? , young-of-the-year). Sub-adults were fish between 1 and 2 years old that were primarily immature and so not considered as reproductively active. Finally, adults were fish of 2 ? years old and were considered as mature fish, capable of reproduction. This enabled each of the stages to be treated separately within the modelled dispersal process .
For each year in the simulation, the probabilities of juveniles developing into sub-adults, sub-adults developing into adults, and adults reproducing were set to unity, i.e. the event occurred if the individual survived. Reproduction was modelled using a transition matrix for each of the three stages defined above (Caswell b Fig. 1 The lower catchment of the Great Ouse river in eastern England, showing bitterling sampling sites (black circles) and potential barriers to their upstream dispersal within the catchment (grey triangle; weirs, locks, pumping stations, sluices and aqueducts). The river channel is in light grey, other than the area where bitterling were initially captured in 1984 (indicated by a black open circle). Co-ordinates are of the Ordnance Survey national grid, and the inset figure shows the location of the catchment by a black square. Ó Crown copyright and database rights 2018 Ordnance Survey (100025252)  2001), with the parameters being survival, development to the next stage and fecundity (Table 1). The latter was set as the number of offspring per female that survive to one year old at quasi-zero population density, and was subject to density-dependence (lower rates for higher population densities; Neubert and Caswell 2000). Survival by stage was also assumed to be density-dependent, and was weighted such that juveniles had no effect on later stages, and sub-adults had only 10% of the effect on adults as adults had on each other and on sub-adults. The density-dependence attributed to the two previous parameters was due to the negative effect that a large number of parasitizing bitterling eggs can have on the quality of mussel gills, where mussel performance is reduced and bitterling would not be able to reproduce again on those mussels with reduced fitness (Mills et al. 2005;Reichard et al. 2007;Smith 2017).

Dispersal parameters
In RangeShifter, dispersal is modelled in three phases: emigration, transfer and settlement, so that dependencies can be added to each phase separately . Emigration was set to zero for juveniles, which stay within mussel gills for a proportion of their first year of life. For older stages, emigration was low if density was below a certain inflection point and higher if density was above the inflection point. We assumed that the dispersal probability of the sub- ABC* in the 'Value' column denotes parameters that were estimated by approximate Bayesian computation adults would be low due to high mortality risks at this life stage, and therefore the adult fish would be the main dispersers. The transfer of individuals was modelled using the stochastic movement simulator (SMS). This models how the individual moves on a cell-by-cell basis, as determined by relative costs (in the sense of the least cost path approach) and a tendency to maintain a correlated path (directional persistence) (Palmer et al. 2011;Coulon et al. 2015;Samson et al. 2017). Our habitat quality raster was incorporated, so that a dispersing individual would be more likely to move into the wider of two branches when reaching a confluence. A per-step mortality constant was applied so that individuals would not move indefinitely if they did not find a suitable patch for settlement. Finally, settlement probability in a nonnatal patch was considered as density-independent, but less than 1.0 ( Table 1).

Estimation of parameters by ABC
Prior distributions of the eleven parameters to be estimated by ABC (Table 1) were generated using information extracted from the literature and expert opinion. Given that bitterling life history traits have not been studied extensively, information extracted from literature was mainly based on the parasitic relationship between mussels and bitterling, so that the model should incorporate some density dependence in fecundity to account for this (Mills et al. 2005;Reichard et al. 2007;Smith 2017). Otherwise, priors were based on author opinion from their experience of working on other small-bodied invasive cyprinid fish, such as Pseudorasbora parva (e.g. Britton et al. 2007Britton et al. , 2008Britton et al. , 2010Figs. 2, 3 and 4; Table 3). A total of 250,000 parameter combinations was sampled independently from the prior distributions and, for each of the parameter combinations, five replicate simulations were run. These simulations were initiated using the bitterling distribution in 1984 (using 1983 as year of starting simulation; Fig. 1), and finishing in 2018. Predicted patch-level presence and pre-reproduction sub-population sizes averaged over the five replicates for each simulation were compared with observed presence data and fish density estimates. A distance metric was computed to determine how close the predictions of the model given the sampled parameters were to the actual time-series of range expansion. We adapted the distance metric q of van der Vaart et al.
(2015) by introducing a weighting for each observed value, so that the model fit for the ith sample set becomes where D j is the observed value for empirical data point j given weight w j , m i,j is the corresponding predicted value from the model for parameter sample set i and the standard deviation sd(m j ) is a scaling factor to allow for the observed data points to be made at different scales (sub-population count up to many thousands, presence constrained to lie between zero and one). We used 620 observed data points. These comprised 394 presence/absence observations and 226 sub-population estimates. We down-weighted 220 of the presence/absence observations (56%) for which we had assumed absence (e.g. that a patch was not occupied in the years immediately preceding the first observation of bitterling within it); weightings were reduced from 1.0 by 0.1 in successive years from the actual observation up to a maximum reduction of 0.5. A ranking of the 250,000 distance metrics was then generated and the best-fitting 250 retained to provide the posterior probability distributions. These 250 samples provided the credible intervals for the estimated parameters based on the observed range expansion data, enabling comparison of prior versus posterior distributions (Csillery et al. 2010).

Goodness-of-fit
We estimated the goodness-of-fit of the model in two ways: (1) predicted presence/absence was averaged over all 250 parameter sets and a single goodness-offit statistic was calculated on the basis of one set of predictions encompassing parameter uncertainty, and (2) owing to parameter uncertainty, a goodness-of-fit statistic for each parameter set was calculated and then an average statistic and uncertainty was determined around it. In both cases, the true skill statistic (TSS) was used, which has been recommended for evaluating the accuracy of species distribution models (Allouche et al. 2006). The TSS takes a value from -1 to ?1, where zero equates to a fit no better than random and ?1 indicates a perfect fit.

Prediction of future range expansion
Following this model-fitting process, prediction of the future expansion of bitterling was simulated, starting from their initial detection in samples in 1984 (using 1983 as year of starting simulation) and running for the next 100 years. The starting point of 1983 was used instead of the most recent observed distribution (2018), as RangeShifter requires all patches to be initialised at the same density. Thus, had 2018 been used as the starting point, then simulations would have been based on that year's mean patch density and ignoring the high spatio-temporal variance in population sizes at the range front, thereby altering the patterns of density-dependent emigration and settlement in the years following initialisation.
All the parameters used during this simulation were set up as in the posterior distribution, i.e. 250 simulations were run which would give predictions allowing for model parameter uncertainty. Standard error and confidence intervals for 90% and 95% of patch occupancy from the 250 predictions were also calculated. Those percentages correspond to when the catchment is considered to be fully colonised. The Fig. 2 Prior (black) and posterior (grey) distributions of the demographic parameters: a rate of density dependence (1/b), b fecundity (ø), c stage-dependent survival rates (c1 juveniles, c2 sub-adults, c3 adults). PDF probability density function bootstrapping resampling technique was used for this purpose, as the data did not meet parametric assumptions (Mooney and Duval 1993;DiCiccio and Efron 1996;O'Hagana and Stevens 2003). All statistical analyses were conducted using R 3.5.1 (R Core Team 2018).
The data used to develop and calibrate the model are available from the Environment Agency (2018).

Model fitting
The posterior distributions for six of the parameters were similar to their priors, i.e. the empirical data provided little additional information upon which to reduce parameter uncertainty. For the other five, however, there were varying degrees of difference   (Figs. 2, 3 and 4). For rate of density dependence (1/b), sub-adult survival and adult maximum emigration probability, the posterior estimates were somewhat higher than priors. For perstep mortality probability, posterior estimates were lower than priors. Most notably, for sub-adult maximum emigration probability, there was a substantial difference between the relatively high maximum emigration probability posterior estimates and the prior assumption of extremely low probability. These five posterior distributions showed some tendency to be inter-correlated, especially sub-adult maximum emigration probability, for which low values tended to be associated with high values of 1/b and low values of per-step mortality probability (Table 4).
The goodness-of-fit of the model as calculated by method 1 yielded a TSS value of 0.759, and method 2 yielded a mean TSS value of 0.728 (90% confidence interval 0.685 to 0.765). The applied goodness-of-fit revealed that the model presented here was able to reproduce the observed pattern of colonisation of the catchment with a relatively high degree of accuracy. Moreover, a more accurate fit to the observed pattern was obtained if the 250 samples of our posterior distribution were treated as a collective whole than as individual predictions.

Simulating bitterling range expansion
From the initial records of bitterling presence in 1984, their predicted early spread matched relatively well with the observed time-series occupancy (Fig. 5a), although the model was unable to replicate the period of near stasis in the sampled data between about 1992 and 2007. Maps of model deviance of bitterling presence at the patch level that were averaged over the whole modelled period and for the decade of near stasis (supplementary material, Figs. 7 and 8) both demonstrated spatially correlated patterns, which suggested some potential influence of spatio-temporal variation in sampling effort and possibly also the effect of a sluice acting as a barrier in the southwestern part of the catchment. Predicted patch occupancy showed a sigmoidal pattern; while it took 20 years for bitterling to occupy 20% of the catchment, they are then predicted to only require a further 30 years to achieve 80% occupancy (2030) (Fig. 5b). With 95% confidence, it was predicted that 90% of the patches would be occupied after 61 to 63 years (2044 to 2046) and 95% occupied after 69 to 71 years (2052 to 2054) ( Fig. 6; Table 2).

Discussion
We have demonstrated how approximate Bayesian computation can be applied to estimate the parameters of a mechanistic simulation model for predicting the future spread of an invading alien species. Our prior knowledge of some of the model parameters was imprecise, but by combining that knowledge with observed data on the spread of the case-study species to date, and with estimates of its local population density, we were able to refine the parameter estimates and, just as crucially, allow for inter-correlations between them. The refined parameter estimates enabled predictions of future occupancy of the catchment to be made with a relatively high degree of precision. The model revealed a sigmoidal pattern in temporal patch occupancy by bitterling in the catchment, with predictions of 95% patch occupancy after 69 years. The implications of this model are now discussed in relation to the insights gained on bitterling dispersal by the model, the performance of the model in predicting the temporal and spatial pattern of bitterling dispersal, and finally how the model provides important insights into invasion management.

Model insights into bitterling dispersal
Our prior distributions of model parameters were mainly based on expert opinion of the ecology of small-bodied cyprinid fishes, as bitterling population biology is relatively data poor because of their negligible fishery value and interest. The only exception was that their mode of reproduction involves parasitism of mussels, resulting in some density dependence in recruitment (e.g. Mills et al. 2005;Reichard et al. 2007;Smith 2017). That the prior bitterling population biology data were limited was not unusual, as low value aquatic species often lack empirical data on their populations (Karakus et al. 2018;Tarkan et al. 2018). Thus, when these species are introduced into a new region, whether intentionally or accidentally (Gozlan et al. 2010a), the data available for predicting their invasiveness are often limited , resulting in poorly informed model parameters Urban et al. 2016) and final models with high uncertainty (Parry et al. 2013). We showed that monitoring programmes can provide dispersal time series that, when coupled with ABC methods, can help overcome this lack of prior data, enabling more robust predictions of the dispersal dynamics and future invasiveness of the modelled species (Neeson et al. 2012;Barros et al. 2016;Samson et al. 2017).
The ABC routine thus enabled the values of the data-poor model parameters to be predicted in a more robust manner (van der Vaart et al. 2015(van der Vaart et al. , 2016. Comparison of the prior versus posterior distributions of these parameters revealed that some posteriors differed little from their priors, suggesting that expert opinion appropriately informed the model priors. There were, however, five model parameters that were strongly informed by the dispersal time series and thus were poorly informed by expert opinion. Of these, the relatively high posterior estimates of maximum emigration probability of sub-adult bitterling (D) were particularly interesting. The prior distribution of D for sub-adults was based on the assumption that the dispersal probability of the sub-adults would be low, due to high mortality risks at this life stage. This assumption was formed due to the combination of asocial, individual fish often being the dispersers at the invasion front (Cote et al. 2010) and smaller-bodied individuals having higher predation risks (Nilsson and Brönmark 2000). Thus, small-bodied, sub-adult bitterling were assumed to be relatively sedentary to  Standard error (SE) and 95% confidence intervals (CI) were calculated by the adjusted bootstrap percentile method across all 250 simulations from the posterior parameter distribution maximise their survival, and any that did disperse would have a high probability of predation. Indeed, in invasive round goby Neogobius melanostomus, upstream-directed range expansion was led by the movement of larger bodied individuals of high trophic positions, rather than juveniles that were leaving high density areas due to, for example, high competition (Brandner et al. 2013). However, when we used low emigration rates of sub-adults, patch occupancy did not match the observed occupancy time series. Instead, the ABC component of the IBM showed that there should be dispersive behaviours by both subadults and adults driving this bitterling invasion. In part, this may be due to an artefact of RangeShifter, whereby an individual which has dispersed may not do so again. If the model allowed dispersal by adults only, then there would inevitably be a delay of two years between the colonisation of a patch and the production of the next wave of dispersers. We could compensate by making patches longer, but that would result in reduced spatial precision. The age-specific dispersal preferences of bitterling detected by our model have also been detected in other fishes (Frank 1992;Stiver et al. 2007). For example, high levels of gene flow in Lethrinus nebulosus were assumed to result from high adult dispersal, yet models indicated that it was larval dispersal, not adult, that caused the gene flow patterns (Berry et al. 2012). Correspondingly, the use of ABC provided an important insight into bitterling stage-specific dispersal and suggests some counter-intuitive dispersal patterns at the invasion front that warrant further empirical investigation.

Predictions of the bitterling dispersal pattern
The comparison of simulated versus actual time series data enabled the future development of bitterling invasion to be predicted by the IBM with relatively high confidence. In periods of rapid range expansion, individuals at the range front often show rapid population growth, facilitated by individuals investing heavily in somatic growth and reproduction when density dependence process are rarely apparent . For example, in N. melanostomus, individuals at the invasion front gained dispersal advantages by attaining large body sizes relatively quickly, facilitated by low competition (Brandner et al. 2013). With 20% occupancy of the catchment after 20 years, the somewhat higher colonisation rates thereafter were likely to be due to the population gaining sufficient distribution that it was then able to increase its occupancy of the catchment relatively quickly in both upstream and downstream directions, i.e. it had reached a level beyond which dispersal into a larger number of sites was possible in a relatively short timeframe. However, there was a period of relative stasis in the observed rate of colonisation lasting about 10 years, which the model was unable to capture. Environmental factors acting on key dispersal processes and parameters may have been important in this and would require further exploration in future model development. At present, however, we lack a detailed understanding of the effects of environmental factors on demographic and (especially) dispersion rates, and such relationships have yet to be incorporated into RangeShifter. Notwithstanding, environmental factors such as water temperature and flow rates in the first summer of life are recognised as important determinants of annual recruitment rates in other riverine cyprinid fishes in England and so might also be important in bitterling population dynamics (Nunn et al. 2007;Beardsley and Britton 2012).
While the IBM was able to provide a series of important insights into how bitterling, as a model small-bodied invasive fish, might disperse through a lowland river catchment, it was also apparent that issues remain with the model that could potentially be improved. The ABC routine could, for example, be refined, especially with regard to the need for the averaging of predicted values over a number of replicates to allow for RangeShifter being a stochastic model . This makes it more difficult for the model to capture the wide variation in adult population densities observed at many of the sampling sites in the years immediately following colonisation. Also, the manner in which the river environment was represented could be improved. For example, the River Great Ouse has a relatively linear river channel whose primary purpose is the flood and drainage management of the surrounding agricultural land. Whilst its separation into a series of patches in the model that accounted for the presence of artificial barriers (weirs, locks, etc.), in the areas away from these barriers, the patches were based mainly on size (2 to 3 km of river length). This was mainly to assist the modelling process and underlying assumptions. This meant, however, that patch delimitation was based on modelling requirements rather than on knowledge of bitterling population demographics and dispersal abilities. For this to be overcome would, however, require data on the individual movements of bitterling across different stages, something that remains technically difficult due to their small body sizes that makes the use of some common telemetry methods highly challenging (e.g. Klinard et al. 2018). Nevertheless, the relatively accurate predictions of the bitterling dispersal pattern were similar to the empirical data. This suggests that the model has high applied utility for simulating the outcome of, for example, management interventions that aim to inhibit their invasion.

Implications for invasion management
The results of this study have highlighted that IBMs have high utility for gaining knowledge on the dispersal processes of aquatic invaders that are difficult to obtain from empirical data collection alone and can be used to help develop more informed management practices (Sakai et al. 2001;Grimm et al. 2006;Samson et al. 2017). Indeed, even without completing any further simulations, the IBM results suggest there were two opportunities for management interventions to have been implemented on the Great Ouse that could have inhibited the bitterling invasion. The first would have been immediately following their initial detection in a very restricted spatial area in 1984, as management interventions are easier and more effective when the extent of invasion is limited (Pyke et al. 2008;Britton et al. 2011). The second opportunity would have been the 20-year period of low colonisation rates, although this would have been more difficult than previously due to the greater spatial extent of catchment occupancy. It is, however, acknowledged that eradicating or even controlling populations of invasive fishes in open systems is highly challenging (Britton et al. 2011;Davies and Britton 2015). Management interventions on this scale also usually require rapid implementation (Pyke et al. 2008), supported by robust invasion risk assessment processes (Copp et al. 2009). Management interventions are then usually only implemented on those invaders assessed as relatively high risk (Britton et al. 2011). Correspondingly, the lack of initial management interventions in 1984 were likely to have resulted from the paucity of invasion assessment tools available at that time, coupled with no predictive assessment of the potential extent of their invasion. Given the extent of their range today, even if invasion risk assessments suggest some population control is required, it would most likely be prohibitively expensive and/or have a low likelihood of success (Britton et al. 2011).

Conclusions
A relatively complex IBM was developed here that enabled key invasion processes, such as dispersal, to be incorporated into model fitting using an ABC regime. The model revealed that whilst a range of different combinations of parameter values fitted the observed time series data, nevertheless, it delivered some important predictions into the dispersal dynamics of bitterling. Thus, the approach delivered novel insights into the ecological behaviours and dynamics of this invader, with the model improving our ability to predict, and ultimately manage, successful invasive species.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.  Symbols for the parameters are defined in Table 1 with sub-headings 0, 1 and 2 corresponding consecutively to the fish stagestructure: juveniles, sub-adults and adults Bold font denotes combinations that are significant (P \ 0.05), **denotes (P \ 0.01)