Early engagement of stakeholders with individual-based modelling can inform research for improving invasive species management: the round goby as a case study

Individual-based models (IBMs) incorporating realistic representations of key range-front processes such as dispersal can be used as tools to investigate the dynamics of invasive species. Managers can apply insights from these models to take effective action to prevent further spread and prioritize measures preventing establishment of invasive species. We highlight here how early-stage IBMs (constructed under constraints of time and data availability) can also play an important role in deﬁning key research priorities for providing key information on the biology of an invasive species in order that subsequent models can provide robust insight into potential management interventions. The round goby, Neogobius melanostomus , is currently spreading through the Baltic Sea, with major negative effects being reported in the wake of its invasion. Together with stakeholders, we parameterize an IBM to investigate the goby’s potential spread pattern throughout the Gulf of Gdansk and the Baltic Sea. Model parameters were assigned by integrating information obtained through stakeholder interaction, from scientiﬁc literature, or estimated using an inverse modeling approach when not available. IBMs can provide valuable direction to research on invasive species even when there is limited data and/or time available to parameterize/ﬁt them to the degree to which we might aspire in an ideal world. Co-development of models with stakeholders can be used to recognize important invasion patterns, in addition to identifying and estimating unknown environmental parameters, thereby guiding the direction of future research. Well-parameterized and validated models are not required in the earlier stages of the modeling cycle where their main utility is as a tool for thought.

Individual-based models (IBMs) incorporating realistic representations of key range-front processes such as dispersal can be used as tools to investigate the dynamics of invasive species. Managers can apply insights from these models to take effective action to prevent further spread and prioritize measures preventing establishment of invasive species. We highlight here how early-stage IBMs (constructed under constraints of time and data availability) can also play an important role in defining key research priorities for providing key information on the biology of an invasive species in order that subsequent models can provide robust insight into potential management interventions. The round goby, Neogobius melanostomus, is currently spreading through the Baltic Sea, with major negative effects being reported in the wake of its invasion. Together with stakeholders, we parameterize an IBM to investigate the goby's potential spread pattern throughout the Gulf of Gdansk and the Baltic Sea. Model parameters were assigned by integrating information obtained through stakeholder interaction, from scientific literature, or estimated using an inverse modeling approach when not available. IBMs can provide valuable direction to research on invasive species even when there is limited data and/or time available to parameterize/fit them to the degree to which we might aspire in an ideal world. Co-development of models with stakeholders can be used to recognize important invasion patterns, in addition to identifying and estimating unknown environmental parameters, thereby guiding the direction of future research. Well-parameterized and validated models are not required in the earlier stages of the modeling cycle where their main utility is as a tool for thought. Keywords: individual-based model, RangeShifter, round-goby, invasive, stakeholder, dispersal, management

INTRODUCTION Invasive Species and the Need for Ecological Modeling
Invasive species are one of the driving forces behind biodiversity loss, and their persistence in non-native areas can result in substantial environmental and economic costs (Pimentel et al., 2000(Pimentel et al., , 2005Molnar et al., 2008;Cardador et al., 2016). Once established, invasive species have the potential to alter local habitat quality, increase competition for resources, prey on native populations, and spread disease (Kwon et al., 2006;Karlson et al., 2007;Salo et al., 2007;Crowl et al., 2008;Gallardo et al., 2016). As a result, the management and control of invasive species has been a central research focus for many years, and a priority for biological conservation.
There is a continual need for the development and improvement of both new and existing conservation management strategies either to control the spread, reduce biomass or, if possible, to eradicate an invasive species from its non-native environment (Ojaveer et al., 2015). However, implementing management procedures can be costly, both economically and environmentally (Hulme, 2009). Therefore, techniques for forecasting the spread of species and assessing the likely impact of alternative management strategies are desirable (Uden et al., 2015;Kotta et al., 2016;Katsanevakis et al., 2017). One such way to evaluate potential management strategies is through ecological modeling (Uden et al., 2015;Goldstein et al., 2016;Kotta et al., 2016). For example, being able to model the spatial distribution of a species accurately can potentially provide numerous facilities, such as predicting future distributions or furthering our understanding of the original invasion process (Adams et al., 2015).

Forecasting Dispersal in Invasive Species through Spatially Explicit Models
The accuracy and utility of process-based models for ecological forecasting has vastly improved over the past few years (Cuddington et al., 2013;Evans et al., 2013;Urban et al., 2016), particularly as the understanding surrounding ecological processes such as dispersal dynamics has increased Goldstein et al., 2016). As dispersal is one of the key determinants of species spatial dynamics, understanding and accurately simulating the dispersal process is central to predicting species spread (Hastings et al., 2005;Bocedi et al., 2014;Brown et al., 2014). Numerous studies demonstrate that dispersal is key to species undergoing range expansion, and that there is selection for increased dispersal propensity at the range front (Travis et al., 2010;Brown et al., 2014;Huang et al., 2015;Myles-Gonzalez et al., 2015;Parry et al., 2015;Therry et al., 2015). For example, in the invasive Cane toad (Rhinella marina Linnaeus, 1758), individuals in the invasion front disperse further, more frequently and in straighter paths than those in established core populations (Brown et al., 2014;Hudson et al., 2015), and even possess physiology that facilitates their dispersal propensity (Phillips et al., 2006). As such, spatially explicit models that incorporate ecological and even evolutionary or physiological complexity can be vital tools in making predictions regarding range extent and the effectiveness of control regimes for invasive species (Higgins et al., 1996;Meekins and Mccarthy, 2002;Vuilleumier et al., 2011;Goldstein et al., 2016). Calibrating and validating such models with suitable data, if available, can provide an excellent opportunity to investigate species-specific invasions, assess invasion patterns and address concerns. However, very rarely (if ever) will all the data required to parameterize a model fully be readily available in the literature. One way of obtaining such information is through stakeholder interaction. Involving stakeholders in the modeling process additionally allows for the continual evaluation of model utility, accuracy and the development of future model applications.

Early Engagement of Stakeholders in the Ecological Modeling Process
Often stakeholders encounter a model only at the stage where it has been tightly parameterized and validated by ecological researchers. Traditional thinking tends to be that a model needs to be well-parameterized and validated before it can be useful in an applied context. Indeed, an often encountered view is that it can be dangerous for a modeler to demonstrate an "immature" model to stakeholders due to risks of losing credibility or of providing unsound advice. However, developing a well-tested model can be a time consuming process, and this is problematic especially when early intervention is often critical for successful management outcomes. It has been repeatedly highlighted that early involvement of stakeholders into ecological management efforts increases chances for success (Bayliss et al., 2013;Seidl et al., 2013) and we consider that models can provide an important tool for thought at this early stage, well before they reach the level of maturity that we would expect them to have reached prior to providing robust management advice. In assessing the potential risks posed by an invasive species, and scoping out potential control options, scientists and stakeholders must first objectively assess where their knowledge might be incomplete (Krueger et al., 2012) and a prototype model can provide an excellent tool for formalizing the process of establishing what is already known, what is not known and, critically, identifying what it is that isn't known that is likely to be most influential in determining the invasion dynamics. Understanding of where key knowledge gaps exist can inform future research and data collection (Voinov and Bousquet, 2010). Furthermore, the importance of clear communication of the results and implications of risk assessment exercises to stakeholders and authorities has been emphasized, increasingly so in recent years, with particular relevance to the management of our marine ecosystems (Katsanevakis et al., 2017;Stelzenmüller et al., 2018). Here, we put this into practice, and emphasize that it can be extremely valuable to engage stakeholders with an early prototype model and use their input to tailor the modeling process to practical needs. We additionally emphasize the value that an early stage model can provide as a means for horizon scanning for potential threats due to the invasive species, and can be used to provide some initial risk assessments Frontiers in Ecology and Evolution | www.frontiersin.org 2 November 2017 | Volume 5 | Article 149 of particular threats (Parrott et al., 2012;Reed et al., 2013;Parrott, 2017).

Case Study: The Round Goby in the Baltic
As a case study, we use our experience of developing an earlystage model for the round goby's spread through the Baltic Sea in order to facilitate stakeholder engagement. The round goby (Neogobius melanostomus Pallas, 1811) is a species, for which ecological modeling can be valuable, firstly for formalizing the process of establishing what we know and what we still need to know and, subsequently, for developing well-tested models that can be used to provide robust management recommendations. This species is native to the Ponto-Caspian region, and has invaded the Great Lakes in North America and multiple locations throughout Europe, most likely as a result of transport through shipping routes via ship ballast water (Kornis et al., 2012;Kotta et al., 2016). The species has been termed "one of Europe's 100 worst invaders" and has in a recent evaluation of 18 taxa of non-indigenous species in the Baltic Sea region been found to be amongst those with the greatest impact (Kotta et al., 2016;N'Guyen et al., 2016). For the past 25 years, the species has been in the process of spreading throughout the Baltic Sea (Sapota, 2004;Schrandt et al., 2016). The first reported sighting was in 1990 in the Gulf of Gdansk, and since then, sightings of the species have been recorded in various areas of the Baltic (Kotta et al., 2016). Whilst some stages of the goby's spread have been well-documented, such as the introduction and invasion of the Gulf of Gdansk (Sapota, 2004) and the inner Danish waters (Azour et al., 2015;Carl et al., 2016), there are other stages of the invasion that are substantially lacking in information.
Here, we highlight how a spatially explicit ecological simulation platform, RangeShifter  can rapidly be used to develop an initial prototype model for early engagement of stakeholders with the process, and subsequently calibrated using spatial data available from the literature and input from stakeholders. We then demonstrate how this intermediate-stage model can be applied to further research in order to identify key data gaps that would need to be filled before a well-tested model could be used to robustly inform management actions.

OVERVIEW OF THE PROCESS
The work described in this paper has been designed to be consistent with the adaptive modeling approach for ecological forecasting outlined in Urban et al. (2016). The overall process of developing the model is outlined in Figure 1. A prototype model of the goby's spread throughout the Baltic was developed and parameterized within a 6 week period through an iterative process (Grimm et al., 2005;Grimm and Railsback, 2012) to present to stakeholders in a symposium context. This period of initial model development was by necessity short in our case, as we had been invited to a round goby symposium to discuss the potential utility of the RangeShifter software in the context of managing the round goby. The description of this initial model development will be kept brief, as it was predominantly an iterative process of altering parameter values and comparing the model output to that of the HELCOM round goby distribution (Michalek et al., 2012). The rapid production of a prototype model allowed demonstration of the potential utility of the model to stakeholders, especially for use in the future after more rigorous assessment. Furthermore, it also provided an overview of what the model could do, which opened the way for suggestions on scenarios and improvements that the model can be used to explore.

Stakeholder Collaboration
Overview of the Symposium A symposium centered around the spread and impact of the round goby in the Baltic Sea was held in Kalmar, Sweden in June 2016. The organization of the symposium was headed by the Swedish Agency for Marine and Freshwater Management 1 , and there were an estimated 30 attendees. The main stakeholder groups consisted of representatives of different levels of local and regional environmental administration, people that participated in a private capacity, and representatives of other groups interested or affected by round goby spread, such as recreational and professional fishermen. The symposium was followed by a workshop focusing on solutions to manage and impede the spread of the round goby throughout the Baltic Sea. During the symposium the overall research project and the model was presented in a 30 min power point presentation (Supplementary Figure 1). The presentation had two main components. First, RangeShifter was presented to the participants along with examples of how the software had already been used to address conservation relevant questions, including invasive species. This was key as a means for establishing our credibility as modelers. Second, the prototype goby model, implemented using RangeShifter, was presented to the stakeholders to demonstrate the potential utility of the model within the Baltic Sea and hence within the geographical focus of the participants' interests. Throughout this second part, we repeatedly stressed both the prototype nature of the model and the fact that while we were in a room full of goby experts, the modelers who had rapidly developed a prototype for demonstration were certainly not.
At the end of the presentation a specific call for input was issued: a slide stating "What we hope to get from you. . . " followed by six suggested inputs: Specific parameters (e.g., demographic and dispersal), The estimated introduction sites (and when), Patterns for comparing model outputs with spatial and temporal patterns of density and sediment type and habitat, Proposed management techniques. Following the presentation there was an open discussion with a call for feedback and input. In transdisciplinary projects it is important that both scientists and non-academic partners contribute on an equal footing N'Guyen et al., 2016). This is especially relevant in the case when inputs are qualitative rather than quantitative. In our case, we were interested in qualitative inputs, and we therefore designed the interaction with stakeholders as open and did not follow a standardized procedure. We felt this would ensure an FIGURE 1 | A schematic representation of the modeling process, from initial literature search to the proposed next steps. Model refinement and evaluation is iterative, reflecting the alterations that are constantly made to the model during the calibration process. Once refined, this model should then be reintroduced to stakeholders for further co-development.
atmosphere that encouraged stakeholders to contribute even anecdotal but possibly relevant information which they might be less inclined to share when e.g., filling out a questionnaire.

Outcomes of the Symposium
The interaction with stakeholders identified essential knowledge gaps, which would have gone unnoticed by us as scientists alone. Crucially the interactions also provided a clear focus in terms of what a useful model would need to include and would need to be able to predict in order that it was most useful to the stakeholders. Also, personal communications with multiple researchers and stakeholders present at the meeting provided an insight into the current understanding of the round goby's spatial presence in the Baltic Sea that was not obvious from searching the literature, including information on new studies that will yield high quality data. Three essential qualitative outcomes of the symposium that were derived from the interactions between modeling team and stakeholders provided strong focus for future work. These related to model building such that key processes driving the spread dynamic are properly represented and parameterized and to developing the model to ensure its relevance for informing key management decisions: First: A knowledge gap regarding the depth of goby dispersal was highlighted as potentially crucial. Prototype model results shown at the workshop included one suggesting that the invasion dynamic is likely to be very sensitive to the depth range over which gobies can disperse. At the workshop attendees noted that adult gobies are sometime caught in deeper water. However, it was suggested that this occurs during winter months and may reflect some adults exhibiting seasonal migration to deeper waters. It became obvious that whether gobies disperse through deep water or disperse solely in shallow areas is currently unknown. Understanding the depth range of goby dispersal may be of great importance to those involved with the round goby invasion for a number of reasons. Depth acting as a barrier to dispersal may be utilized in numerous management protocols to impede or inhibit goby spread into undesirable areas. Furthermore, understanding goby dispersal depth helps to predict future areas that may be under threat of round goby invasion, even without a human-mediated element to the dispersal. Identifying the potential importance of the depth sensitivity of dispersal for patterns of goby spread was a novel outcome of the workshop that will motivate new empirical work.
Second: Threats of the round goby's invasion of the freshwater systems that connect to the Baltic Sea, particularly with regards to Salmonids were identified, as the round goby may devastate their populations through egg consumption (Chotkowski and Ellen Marsden, 1999;Marentette et al., 2011;Ladago et al., 2016). This potential impact of the round goby was a key issue for many of the stakeholders present and highlighted the importance that to be useful for management a model would need to be able to effectively operate into riverine systems and potentially also account for salinity gradients and tolerances.
Third; The threat that the round goby poses to the long-tailed duck (Clangula hyemalis Linnaeus, 1758) was emphasized (Hearn et al., 2015). The Baltic Sea is the key wintering destination for the majority of the western Siberian and northern European populations of the long-tailed duck (Hearn et al., 2015), which currently faces a multitude of threats such as predation, competition, oil spills, gillnets, hunting, habitat destruction, and water traffic (information available on the BirdLife International website 2 ). The round goby and the long-tailed duck share a diet of mussels and crustaceans. Hence, the spread of the round goby to the overwintering habitats may result in competition for food. As the Baltic Sea is the main overwintering area, a reduction in food availability for the long-tailed duck in this area may prove devastating. Consequently, this area is recommended to be a crucial area to protect from invasion by the round goby (Hearn et al., 2015). Currently there are no effective means for estimating the risk that these areas will be invaded. Hence, estimations for whether, and if so when, the goby will reach the overwintering areas from their current distribution would be valuable, to estimate time-scales for conservation efforts for the long-tailed duck, and to design measures to protect the area from the spread of the round goby. We note here that while there was existing information in the literature highlighting the potential impact of round goby on long tailed duck in the duck's key overwintering sites (Hearn et al., 2015), it would have been unlikely that the modeling team would have easily found it. Thus, the stakeholder workshop provided a means for those fully familiar with the system to direct the modeling team to literature relating to the focal species and its potential impacts that isn't primarily about the focal species.

Modeling Population Dynamics and Dispersal in RangeShifter
We used a spatially explicit, individual-based model (IBM), RangeShifter  to simulate the spread of the round goby throughout the Baltic Sea. RangeShifter was developed in response to the demand for integrated dynamic models, and as such, provides a platform with which to model complex population dynamics and dispersal behaviors, at the individual scale (Franklin, 2010;Huntley et al., 2010;Thuiller et al., 2013;Lurgi et al., 2015).
To represent the Baltic Sea, a gridded seascape was created in ArcGIS 10.3.1 using raster data extracted from the EMODnet Bathymetry portal 3 Each cell was 2.5 by 2.5 km and characterized by depth. Population dynamics were modeled at the cell scale. The numbers of individual fish in the Baltic, or even in a local area, at reported densities (Vélez-Espino et al., 2010) would be far too large to be explicitly represented in the model, and therefore we treated a modeled "individual" as representing a localized established sub-population of unspecified size (hereafter "individual" for consistency with RangeShifter terminology), which was regarded as female in a single-sex model. It was not necessary to represent the overlapping generations of the species, but sufficient to model the reproductive rate of such "individuals, " i.e., the rate at which "daughter" sub-populations were produced, some of which would disperse to expand the range of the species.
At model initiation, individuals were assigned to cells within species introduction points at half carrying capacity. In each year, the overall dynamics consists of reproduction, death of adults, and offspring dispersal. Reproduction by each individual is determined by a stochastic draw from a Poisson distribution having a mean set by the maximum growth rate at low density and subject to density-dependent reduction following Maynard- Smith and Slatkin (1973). Carrying capacity, K, was set to 10 3 http://www.emodnet-bathymetry.eu/ individuals/ha for all cells. However, this limitation is unlikely to be critical for the pattern of range expansion on which we were focused, given that densities at the range front are expected to be much lower than in long-established areas Groen et al., 2012;Azour et al., 2015).
Once reproduction has taken place, individuals could emigrate away from their natal cell, an action dependent on the local density within the cell. If an individual left the cell, its trajectory was modeled using the Stochastic Movement Simulator (SMS; Palmer et al., 2011). SMS simulates an individual's path throughout the landscape, in which the direction of movement between cells is based on the relative cell "costs" to movement and on a tendency to follow a correlated path (directional persistence). The perceptual range, in which costs were evaluated, was set at 1 cell (i.e., no more than 2.5 km).

Incorporating the Stakeholder Input into the Model
A key issue that emerged from the stakeholder workshop was a lack of knowledge relating to the depths of water through which gobies can disperse. This issue was, in part, highlighted by some of the runs of the prototype model, demonstrated at the workshop, in which it was clear that including a depth threshold resulted in very different spread patterns than omitting one. Accordingly, cell cost was set in relation to a threshold depth for movement: the cost of traversing a cell of the depth threshold and deeper was set to a very high value, and the cost of traversing a cell above the depth threshold was set to a very low value. In doing so, individuals were much less likely to travel into deeper water than that set by the threshold. For all depths, each step an individual took had an associated spatially and temporally constant mortality risk.
Upon reaching a new cell, an individual had the opportunity to settle or continue movement to a different cell. The decision to settle was density-dependent. If the population density was too high in a cell, then the individual would not settle but continue to disperse to a neighboring cell .

Parameter Calibration and Assessing Model Performance
The majority of the parameters required for the model were not widely available in the literature or through online resources. Consequently, in order to calibrate the model parameters, the Gulf of Gdansk was chosen, as detailed spatial information regarding the goby's spread through the area was available. This spatial information was primarily obtained from the NOBANIS fact sheet, produced by Sapota (2012). NOBANIS is the European Network on Invasive Alien Species, and the project produces information and fact sheets on invasive alien species. The fact sheet, written and referenced by experts, provides a range of information including recommendations for management, species ecology and information regarding its historical introduction and spread. This temporal spatial presence information available in the NOBANIS fact sheet was used as a baseline to calibrate the model.
Parameter values were calibrated using a pattern-oriented modeling approach (POM) (Grimm et al., 2005;Grimm and Railsback, 2012;Bauduin et al., 2016), in which simulations Frontiers in Ecology and Evolution | www.frontiersin.org 5 November 2017 | Volume 5 | Article 149 were run for a variety of values for four key parameters, namely the maximum growth rate, the depth threshold, the perstep mortality risk, and the maximum settlement probability (Table 1), in order to find a combination which most precisely matched that of the historical round goby spread throughout the Gulf of Gdansk. For each simulation, the final model distribution was compared to the actual distribution reported in Sapota (2012). Other more minor parameter values, such as the depth threshold cost, were assigned during the creation of the prototype model, using an iterative process. During this process, the values chosen were arbitrary, and altered until the model output started to match the goby distribution seen in the NOBANIS fact sheet. Therefore, these parameters were used more as values to tune the initial model, rather than parameters that were important to investigate. The model's predicted output was compared to the observed output for each year that data were available, in order to obtain the most accurate dispersal pattern throughout the Gulf of Gdansk.

Accuracy of Model Calibration
To assess the accuracy of the model for each parameter combination, four metrics were used, in addition to visually inspecting the model output. Model specificity (in which both the observed distribution and the model's predicted distribution do not have individuals present in a cell), sensitivity (in which both the actual distribution and the model have individuals present in a cell), the receiver operator characteristic (ROC) curve with the associated area under the curve (AUC), and Cohen's Kappa, κ. The κ statistic represents a way to measure reliability, or precision, and compares the model prediction accuracy with the accuracy expected to occur by chance (Allouche et al., 2006). Sensitivity and specificity both vary from −1 to +1, in which a score of 0 represents no better than chance, and +1 would represent a perfect score. κ can vary between 0 and 1, where 0 represents an agreement no better than chance, and 1 represents a perfect agreement. An accurate model with an AUC score of an ROC plot would be close to 1. A score close to 0.5 would represent a poor model. Whilst the AUC is threshold independent, the other measurements are threshold dependent. The threshold during analysis is the cut-off value used to translate predicted probabilities into a presence or an absence. Consequently, for a predicted probability to be classed as a presence under a high threshold (such as 0.9), a cell would need to be colonized by individuals in 90% or more of replicated model runs given the specified combination of model parameters.
In order to calculate the sensitivity, specificity, AUC and κ for each parameter combination, each combination was repeated over 100 simulations. These metrics were calculated using the PresenceAbsence package in RStudio 3.3.1 (Freeman and Moisen, 2008).

Ports of Introduction and Modeled Population Initiation in the Gulf of Gdansk and the Baltic Sea
The species introduction points of the Baltic Sea, where populations were initiated, were estimated using information available in the literature (Kotta et al., 2016), species presence information from the symposium and shipping port and traffic information available on Baltic Transport Maps 4 Ports of introduction were assumed to be the closest shipping port to a current goby distribution. The initiation of a population at the entry points was staggered in an attempt to replicate the introduction of the goby throughout the Baltic at various points in time. For example, populations were initiated in the Gulf of Gdansk entry points at year 0 (representing the year 1990), but populations initiated around Kalmar were not initiated until year 20 (2010). The timing of the staggered introductions at various points on the map were based on estimates from the literature (reviewed by Kotta et al., 2016). The staggered introductions were carried out using a customized version of RangeShifter that allowed populations to be initialized in individual cells at specified times. Parameter values applied were informed by the results from fitting the model to the Gulf of Gdansk (Table 1).

Model Calibration: Role of Depth Threshold in the Gulf of Gdansk
Model accuracy was most strongly influenced by the depth threshold: 76% of the variance in κ was explained by depth, as compared with 11% by the maximum settlement probability, 4.6% by the per-step mortality risk, 3.3% by the maximum growth rate and negligible amounts by interactions. The model was most accurate for a depth threshold between 10 and 25 m, and accuracy increased slightly with decreasing settlement probability and mortality risk and with increasing growth rate (Figure 2). Similar conclusions regarding the importance of the depth threshold were drawn from the other accuracy metrics (  Figure 3, ranging from good (AUC and k scores close to 1) to poor (AUC and k scores close to 0.5 and 0, respectively). Through the process, a number of models with high accuracy were produced, with some models obtaining accuracy values of 0.8 for all four accuracy metrics, even when the model threshold was high (0.8) (Figure 4).

Model Output: Projections Based on the Role of Depth Threshold, across the Entire Baltic Sea
Despite obtaining a range of accurate parameter combinations for the Gulf of Gdansk, when they were applied to the entire Baltic, the overall model output was poor when compared to the extensive observed distribution spanning a substantial proportion of the Baltic coastline as reported in the literature (Figure 5). The accuracy scores calculated for The Baltic suggest that the model was not much better at predicting the goby distribution than chance (AUC scores close to 0.5, and other scores close to 0.1).

DISCUSSION
In this study, we rapidly developed a prototype model of round goby spatial dynamics that was used to facilitate early engagement with stakeholders. We subsequently combined data available in the literature and stakeholder input in order to calibrate the IBM such that it simulated the round goby's spread throughout the Gulf of Gdansk to a high level of accuracy. We then used the calibrated model to simulate its spread through the Baltic Sea, despite the limitation of imprecise and potentially inaccurate presence data. Our experience demonstrates the value of involving stakeholders early in the modeling process. Prototype model results had indicated that predicted spread was highly sensitive to the inclusion of a depth threshold for dispersal, and the subsequent stakeholder communication highlighted how little is currently understood about goby dispersal at various depths. Consequently, various depth thresholds were incorporated into the modeling, in order to assess the impact of depth on model accuracy and therefore goby dispersal. We demonstrated how, by using known spread patterns, it can be possible to use the model to infer details of the dispersal process, in this case related to the depth threshold of goby dispersal. In detail, we could show that that the limit to dispersal depth of the round goby lies between 10 and 25 m. Empirical data are now required on the depth sensitivity of dispersal such that a robustly parameterized model can be used by the stakeholder/modeler grouping in further steps toward identifying management options. The involvement of stakeholders as early as possible in the process and their regular inclusion throughout as co-developers of the modeling will facilitate a cooperation between scientists and stakeholders in putting possible management measures into practice.

Stakeholder Collaboration-Putting Theory to Practice
Research has identified that the long lag time between research and its publication hinders managers of biological invasions to make use of important results such as our models generated (Matzek et al., 2015). In addition, theory predicts that the success rate of management should be higher if stakeholders and scientists engage early on in the transdisciplinary process of managing an invasive species (Hirsch et al., 2016a;N'Guyen et al., 2016). The main reason behind this is that scientific results that were co-produced by relevant parties in a transdisciplinary process should have better social acceptance and higher compliance by decision makers . In our study, we put these theoretical predictions into practice and engaged in a modeling process that used stakeholder input as an essential component. Stakeholders provided two essential inputs regarding future model optimization: providing information on where higher quality distribution data would be available in the near future and on the priority of including depth in the model. Stakeholders contributed their knowledge and understanding on an equal footing. In an excellent recent contribution on how to co-develop models with stakeholders effectively to address pressing ecological problems, Parrott (2017) argues that it is important for the modelers to get to know the study system well before meeting with stakeholders. Parrott (2017) writes, "Knowing the system well is a key to gaining the trust and confidence of stakeholders in the ability of the modeler and the entire research team to contribute meaningfully to the issue. If the researchers are not from the area, they should spend time visiting and getting to know the region before initial meetings with stakeholders." We had been approached by stakeholders and asked to present the modeling software at a meeting on the threat posed by round gobies to illustrate what might be possible in terms of using RangeShifter to inform management of the species. We only had a few weeks ahead of the meeting in which we were able to build a prototype model for the goby and were thus unable to acquire substantial knowledge of the system prior to meeting stakeholders. However, at the meeting we were able to demonstrate our credibility as ecological modelers by first providing examples of how the RangeShifter was being used to address a range of other applied issues, including landscape management to conserve African forest birds (Aben et al., 2016), assisted colonization of butterflies in Finland  and the invasion of American mink (Neovison vison) in Scotland (Fraser et al., 2015).

Acknowledging the Different Roles of Scientists and Stakeholders
A potential advantage of the approach we took in this study is that the stakeholders naturally take the role as the species/system experts, and the potential risk whereby stakeholders perceive that the researchers assume the role of experts and tell them how their system works is reduced. One potential disadvantage of such an approach is that researchers cannot glean data from stakeholders in the form of quantitative assessments through e.g., specifically designed questionnaires. This disadvantage, however, is compensated by the fact that stakeholders can contribute their knowledge freely through unstructured interactions with researchers. For that, it is clearly critical that the modeling team gain the confidence of the stakeholders, but that need not be by having acquired detailed understanding of the particular study system in advance of a first meeting. Indeed, we suggest that the effective establishment of a model co-development group may be facilitated if this is actually not the case and at the start of the process there is a clear division of expertise between modelers and stakeholders. As the process of co-development of a model proceeds, both researchers and stakeholders can build upon this first interaction on an equal footing albeit with quite different expertise. Our study provides a practical example for future model building efforts on how to rapidly initiate transdisciplinary projects, which is absolutely vital if models are to be successfully used to inform early intervention against invasive species.

Model Calibration
Calibrating the model with precise spatial data produced a highly accurate model that simulated the spread of the goby throughout the Gulf of Gdansk over an 11 year period. The model outputs obtained from the calibration process highlighted the key role of the depth threshold to movement. However, when scaled up through space and applied to the whole of Baltic Sea, the model failed to predict a distribution similar to that observed in the literature. The failure to produce a model for the Baltic Sea with a high degree of accuracy has several implications. One of the main downfalls of the Baltic model seems to occur from uncertainty regarding introduction points. In order to obtain a predicted presence from the model that was similar to that of the observed presence, further introduction points would need to be added, if the parameters obtained from Gdansk were to be used. Although short-distance (∼30 km/year) active migration appears to occur in some local areas (Azour et al., 2015), this suggests that, at the scale of the Baltic sea, the goby did not disperse over long distances as a primary mode of invasion, but that human-mediated transport, for example via ships or other means, was the primary cause of invasion. As large ports were used in the model as the introduction points, this may also suggest that the goby was introduced to various areas that were not necessarily large commercial ports, but also small recreational ports. Subsequently, future efforts to manage the spread of the goby may benefit from focusing preventative measures on human-mediated transport, such as the cleaning of recreational boats (Hirsch et al., 2016c). This will be particularly important in protecting regions that would otherwise be likely to be out of the range of goby colonization due to their being effectively isolated by channels of deeper water.
Furthermore, the presence data used to produce the observed map for model calibration was at a coarse spatial scale. It may be that the goby's presence at various depths in the Baltic was not represented in the observed distribution at a fine enough resolution for accurate model assessment. Given more precise presence data, at a finer resolution, the accuracy of the models predicted goby presence in the Baltic Sea could improve substantially. One of the benefits of such models is the ability to identify on which future data collection efforts should be focused. This is in agreement with the recent call for mandatory catch records and citizen science programs in order to collect data on the round goby (Ojaveer et al., 2015). In the case of this modeling FIGURE 4 | An example of the plots used to assess the accuracy of parameter combinations, using kappa, specificity, and sensitivity measures. The accuracy measures vary from zero to one, in which a value of one represents a perfect accuracy measure and zero a poor one. The cutoff threshold represents the number of repeat simulations a cell was required to have been colonized, in order to be characterized as a presence in the final model evaluation, with 1.0 being 100% and 0.0 being 0% of repeats. (A) represents a set of parameter combinations that predict a goby presence close to that of the observed presence from the literature. In comparison, plots (B-D) demonstrate a decrease in model accuracy.
exercise, presence data over various depth distributions, and the identification and incorporation of the correct introduction points, have been identified as being critical for accurate model calibration.

Depth Sensitivity
In order to replicate the observed goby distribution throughout the Gulf of Gdansk, a dispersal depth limit of ∼20 ± 5 m produced the most accurate model. It is nevertheless important to note that this was calibrated using one area of the Baltic Sea. Thus, to obtain more accurate results, presence data spanning various depths over more locations in the Baltic Sea are required. Hitherto there have however not been any studies dedicated to investigating this aspect of the biology of the species. Furthermore, as round goby is not a commercial species, no catch-related depth information is available from the fishery. The sparse information that exists is from a Polish young fish surveys program, showing that, although generally considered a shallow water inhabitant, high catch rates occur at 50-60 m depth during winter months (November and January-March) (Grygiel, 2007). This suggests that during the cold season, the fish is wintering in deeper sea areas, but whether dispersal occurs during this period or when the fish resides in more shallow, coastal waters remains speculative. The present modeling exercise thus indicates that future research efforts should prioritize obtaining presence and absence data for round goby at various depths throughout the Baltic, and investigate whether dispersal to novel areas occurs during the warm or cold season. Although often expensive and time consuming to collect, this type of information has been achieved for several species though tagging studies (e.g., Boje et al., 2014). Furthermore, compilation of existing data from various national and international surveys and monitoring programs (e.g., the biannual Baltic International Trawl Surveys, BITS) could prove to be a cost-efficient way to obtain essential information. The depth threshold of round goby dispersal is an essential parameter not only for calibrating models, but also for FIGURE 5 | Comparison between (A), the observed goby distribution available from the literature, and (B), an example predicted distribution obtained from the model. The color of the cell represents the presence of a population in a cell, and therefore its colonization in a repetition. The presence varies from one to zero, with a value of one meaning the cell was colonized in every repetition and a value of zero meaning the cell was never colonized. The values along the axis represent the cell numbers of the landscape grid used in the modeling exercise, in which each cell size was 2.5 km.
incorporating into risk assessments of the species spread, both generally and for areas of special interest.

Salinity Tolerance and Ecological Parameters Influencing Spread
Although not identified by the stakeholders in the present study, parameters besides depth should be evaluated for their potential relevance for dispersal tendency. Charlebois et al. (2001) highlighted the need for research determining "dispersal mechanisms and environmental characteristics that limit dispersal." Round goby is considered a euryhaline species, which is able to adapt to salinities ranging from freshwater to brackish conditions. Previous studies have suggested that round goby will not endure oceanic conditions (i.e., high salinity) (Ellis and Macisaac, 2008;Karsiotis et al., 2012). A recent study acclimating round goby to salinities spanning from fresh to seawater has shown that slow increases in salinity (5 PSU per week) to salinities approaching oceanic conditions (30 PSU) severely affected the osmoregulatory capacity of round goby. Although survival was also reduced at oceanic salinities, still 61% of the fish survived at 30 PSU. So while salinity will likely not act as an effective barrier, it might still impede the ongoing spread of round goby through the salinity gradient from the brackish Baltic Sea and into the oceanic North Sea and this warrants its inclusion into dispersal models (Behrens et al., 2017). Further parameters which could turn out to be relevant depend on the study system and could include temperature (thermal limits in round goby are between 0.5 and 26 • C (Chekunova 1974 cited in Charlesbois et al., 2001) and, in running waters, flow velocity (round goby show a critical swimming speed of 35.5 cm s −1 ; Tierney et al., 2011). Recent research suggests that population niche modeling in combination with climatic parameters might benefit from the introduction of thresholds for certain environmental parameters (Almpanidou et al., 2016). Incorporating a minimum of climatic suitability might allow coupling of dispersal models with models of population establishment (Almpanidou et al., 2016). Understanding the interplay of population dynamics and dispersal is relevant for selecting population management options in newly identified populations .

Personality-Dependent Behavior as a Model Parameter and Management Options
Not only the abiotic environment, but also personalitydependent behaviors can be important at the invasion front, where local sub-populations consist mostly of bigger/older asocial individuals (Thorlacius et al., 2015). Recent research has found that personality traits can inform models of dispersal such that only individuals showing trait values above a certain threshold are predicted to disperse . In combination with the depth thresholds, such an approach can complement future models to achieve an even higher accuracy in predicting dispersal.
Until further information is available, our modeled depth trial results may be used as a preliminary guide to assess management regimes and prioritize management areas for vulnerable species. For example, from an applied perspective, the model results raise the prospect that artificial deep channels may stymie the spread of the species. Telemetry-based data on the spread of invasive crayfish in a Central European large lake has also suggested a spread along the shoreline down toward a certain depth isocline. This might make it plausible to slow the natural spread by barriers (Hirsch et al., 2016b). In Lake Tahoe, USA, invasive bivalves have been successfully controlled by the installation of gas impermeable benthic barriers (Wittmann et al., 2012). These examples demonstrate how knowledge of the spatial spread of an invasive species can directly inform its management.

Practical Model Application for Protecting the Long-Tailed Duck
In a practical application of this transdisciplinary approach, we designed a preliminary modeling experiment as an example of how detailed models developed with stakeholders can inform risk assessment of invasive species and help to identify priority areas for management. A key issue that emerged through the stakeholder interaction is the implication of the goby's invasion of the over-wintering habitat of the western Siberian and northern European wintering populations of the endangered long-tailed duck. The populations of the duck may be threatened by the round goby through exploitative competition for food (Hearn et al., 2015;Skabeikis and Lesutiene, 2015;BirdLife International website). In a preliminary trial, we used the calibrated model to demonstrate how an effectively parameterized model could be used to assess whether the long-tailed duck overwintering habitat was at risk of colonization from the round goby in the future. Again, we did this by running the model over a number of years and a number of depth thresholds. This produced a number of scenarios in which the overwintering habitat of the long-tailed duck was invaded, but the time it took for the invasion depended greatly on the depth threshold at which the round goby was able to disperse. Given the current uncertainty surrounding the results of these initial trials or the risk to long-tailed duck populations, and the potential influence results from them might have, we decided that it was premature to publish the results at this stage. However, whilst only being a preliminary experiment, this example reinforces generally (and very effectively reinforced to our co-development team of modelers and stakeholders) the importance of obtaining accurate spatial data regarding the presence of the round goby at various depths.

Future Modeling Perspectives
In this study, we made use of an IBM to simulate the spread of the invasive species. However, it is important to recognize that alternative approaches exist that could equally well be used in transdisciplinary work where models are co-developed to inform understanding and management of invasive species. Indeed, in future studies one valuable approach will be to utilize more than one of these modeling approaches in concert. For example, there can be considerable benefits of jointly developing a stochastic IBM and a typically deterministic integrodifference model to estimate rates of spread (e.g., Travis et al., 2011;Santini et al., 2016). Notably, while until recently integrodifference models have almost exclusively been used to project spread rates across homogenous landscapes, recent developments are enabling rapid simulation of integrodifference equations across spatially complex landscapes (e.g., Synes et al., 2016;Gilbert et al., 2017). One major potential advantage of the integrodifference approach is that the much faster speed of individual simulations will make inverse fitting of parameters through Bayesian approaches including approximate Bayesian computation much more readily achievable. A further important development will be to integrate environmental niche modeling with the population dynamic modeling approaches available.
A key challenge is to move beyond the approach most often taken in what are often termed hybrid species distribution models and to relate the environmental variables directly to the key demographic traits (e.g., reproduction, survival, and dispersal), rather than simply using the environmental niche model to demark suitable and unsuitable environments for a focal species. However, many such relationships have yet to be established in detail (see Zurell et al., 2016 for excellent discussion of key issues). We note here that regardless of the modeling approach taken, in order to engage with stakeholders effectively, it is extremely useful to have clear spatially realistic model output that enables individuals of different backgrounds to relate to the modeling process and its potential (as called for in Stelzenmüller et al., 2018). Thus, as we develop more sophisticated and complex models for predicting and managing spread, we need also to focus on how we develop effective approaches for presenting the results of these models (including associated uncertainties) in an accessible form for those stakeholders with whom we are jointly developing the models, and for others who are likely to find the models useful.

CONCLUDING REMARK
We calibrated an IBM for the round goby, using spatial presence information from the invasion of the Gulf of Gdansk. Stakeholder involvement with question design provided both a preliminary answer and future research directions. It is important that we encourage a culture of publishing work on the process of co-development of models, such that we can learn from one another's successes and failures. This will require more papers, such as this one, that are published at potentially earlier stages of model development and before models are necessarily ready for use to inform management action. In this instance, while short of being ready to inform management action, the model has helped to emphasize the requirement for investment in gathering greater empirical understanding of the depth at which round goby disperse. In the next part of the co-development modeling spiral (Parrott, 2017) this will be gathered, and models will be built using this information, together with higher quality information on human-mediated dispersal pathways, to improve our ability to capture Baltic-wide patterns of invasion and to enable improved forecasts of future distribution under alternative management options to be developed. In general, we promote the increased use of models as a heuristic device for horizon scanning and risk assessment of invasive species and suggest that this utility may be at least as influential as their more traditional usage for informing management at the stage when they are well-validated.  Fund Denmark. In addition, PH and JB would like to thank the Swedish Agency for Marine and Water Management for covering travel expenses to the stakeholder workshop at Kalmar. This study is part of the project PROBIS (BiodivERsA2013) and financially supported by the Swedish Research Council Formas.