Modelling the spatial dynamics of Maui dolphins using individual-based models

The current anthropogenic impacts on nature necessitate more research for nature conservation and restoration purposes. To answer ecological and conservation questions concerning endangered species, individual-based modelling is an obvious choice. Individual-based models can provide reliable results that may be used to predict the effects of different future conservation strategies, once calibrated correctly. Here, we calibrate an individualbased model of Maui dolphin movement, which generates Maui dolphin probability distribution maps. We used sighting data for calibration of the chosen parameter combinations; for each simulation run, collected simulated data was compared to the empirical survey data, resulting in cost (Badness-of-Fit) estimates. Using costs of four different aspects of dolphin behaviour, we estimated the most likely parameter combinations. With optimized parameter values, Maui dolphin probability distribution maps were created, resulting in distributions that fall well outside of the current protection zones where either gillnets or trawling or both are prohibited. With these results, protected areas can be properly adjusted to the estimated distribution of this critically endangered species and so aid in their conservation.


Introduction
For the conservation of species vulnerable to extinction, it is useful to explore the effectiveness of different restoration scenarios. Predicting the effects of different protection measures is of great importance, as policy makers can adjust their actions accordingly and adopt the management strategies that will provide the best possible outcome (Cash et al., 2002). To deliver high quality predictions for possible scenarios, models are needed that represent the study system in a credible, salient, and legitimate fashion (Cash et al., 2002). Yet, many ecological questions are too complex for simple deduction or successful modelling with spatially implicit population models. Custom-made individual-based models (IBMs) are often the best analysis tool for these kinds of questions (DeAngelis and Mooij, 2005;Nabe-Nielsen et al., 2014;MacPherson and Gras, 2016;Van der Vaart et al., 2016;Udevitz et al., 2017;Van Beest et al., 2017).
An urgent case in which high quality model projections are needed is that of the Maui dolphin (Cephalorhynchus hectori maui). The Maui dolphin is, with its current estimated population size of 63 individuals (95% CI 57-75; Baker et al., 2016), classified as critically endangered by the IUCN (Reeves et al., 2008). The declining population (Martien et al., 1999;Slooten et al., 2000;Baker et al., 2013;Wade et al., 2012) that was once distributed along the west and south coasts of New Zealand's North Island (Du Fresne, 2010;Dawson et al., 2001) is now restricted to a much smaller area of approximately 300 km along the west coast, with most recent sightings within an area of about 140 km length (Oremus et al., 2012). Besides their small coastal geographic range and population size, Maui dolphins have a slow reproduction rate, endangering their survival (Davidson et al., 2012). Although the New Zealand government has increased restrictions to fishing activities in a zone along the west coast after the death of a Maui dolphin in 2012 (Hamner et al., 2014), current restrictions are not yet sufficient to result in population recovery.
Spatially implicit, population level models have successfully demonstrated the link between the Maui dolphin population decline and gillnet fishing (Martien et al., 1999;Slooten et al., 2000;Baker et al., T 2013;Wade et al., 2012), brought insight in the maximum by-catch that can be supported (Wade, 1998;Wade et al., 2012), and have prompted policy measures (Slooten and Lad, 1991;Martien et al., 1999;Slooten et al., 2000;Burkhart and Slooten, 2003;Slooten and Dawson, 2010;Slooten and Davies, 2012). However, population level models have lower spatial resolution as to where these policy measures should be in place. To be able to provide policy advice at high spatial resolutions, spatially explicit, individual-based models are needed (Van Beest et al., 2017).
We developed a tailor-made individual-based model of Maui dolphin movement, which can be used to project the impact of fishing activities on the Maui dolphin population under alternative policy scenarios. We model hourly movement of individual dolphins at a 1 km 2 resolution for a duration of 90 days, to account for the small spatial and temporal scales at which dolphin behaviour and human impacts occur in real life. This model will be useful for (i) spatially explicit modelling of interactions between fishery activities and dolphins and (ii) informing the spatially implicit models on values for parameters with high spatial and temporal sensitivity. Dolphin behaviour in our model is based on bathymetry, social drivers, and internal drivers. Because accurate data on the distribution and dynamics of prey species is insufficient, foraging is not explicitly included in the model.

The study system: Maui dolphins
The Maui dolphin population is distributed within a small range along the west coast of the North Island of New Zealand, ranging from Maunganui Bluff to Whanganui (Slooten et al., 2005(Slooten et al., , 2006Currey et al., 2012). The model was set up to represent the habitat of Maui dolphins (Fig. 1a). In the northern part of the Maui dolphin's range, protection implemented in 2008 prohibits gillnetting to 7 nautical miles (nmi) offshore and trawling to a maximum of 4 nmi offshore. Further south, between New Plymouth and Hawera, gillnets were banned in 2012 to 2 rather than 7 nmi offshore and trawling is permitted regardless of distance from shore. The southernmost part of the Maui dolphin range, from Hawera to Whanganui remains unprotected. Dolphins range throughout the area, from the shoreline to the 100 m depth contour, and occasionally enter harbours (Rayment et al., 2011). Maui dolphins were observed to prefer inshore waters in population surveys with equal survey effort with respect to distance from shore (Slooten et al., 2004(Slooten et al., , 2006Rayment and du Fresne, 2007). Dolphin movement behaviours depend on local environmental factors, including water depth, distance to the coast and the presence of other dolphins and fishing vessels . Maui dolphins are not attracted to gillnets, as the fish caught in these nets are too large for them to eat (Miller et al., 2013). But they are strongly attracted to trawlers, often aggregating in groups behind trawling vessels and diving down to the net ).

The model
A full ODD (Overview, Design concepts, Details; Grimm et al., 2006Grimm et al., , 2010 can be found in Supplementary Methods A. The objective of our model is to accurately simulate hourly movements of Maui dolphins as a function of dolphin swimming activity, water depth, distance to their home range centres, and the distribution of nearby conspecifics. We created an individual-based model to simulate the spatial distribution of dolphins over the coastal waters along the west coast of New Zealand's North Island. Fig. 1a represents the area of approximately 180 by 465 km that was used for the spatial configuration of the model. Based on a detailed bathymetry map, we included average ocean depth for each 1 by 1 km grid cell. During a simulation run, dolphins are initially placed in the landscape and iteratively move across by choosing the number of patch displacements per hour and choosing the directions of each of these displacements for a set number of hourly steps. First, 63 dolphins were initially placed within the simulated area. These dolphins are distributed along the coastline, with a doubled probability of initial placement in the centre area (see Fig. 1). Depth at which dolphins occur follows an exponential distribution with λ = -0.05 m −1 (based on field data; initial depth (in m) = ln(x) / -0.05 m −1 , where x is a random number between zero and one). After the initial distribution of dolphins along the coastline, hourly movements were simulated at a 1 km resolution. The number of moves per hour varies for each dolphin and every hour. We assumed that the distribution of the number of moves per hour can be described by a Poisson distribution, of which the scaling exponent μ is equal for all dolphins per simulation. Empirical data on the distance moved per hour by Maui dolphins fits extremely well to a Weibull distribution (see Suppl. Fig. 1). Using a Poisson distribution of the number of steps, with steps of approximately 1 km in 8 random directions, this Weibull distribution can be approximated. The movement direction per discrete move depends on a dolphin's preferences of shallow waters, staying close to its home range centre, and joining other dolphins. For the eight neighbouring grid cells that fall within the 1 km range, the weights of these three factors are multiplied: where W patch is the overall weight of a patch and W depth , W home , and W schooling are the weights assigned to a patch for its depth, closeness to the focal dolphin's home range centre, and the presence of dolphins, respectively. W depth follows a sigmoid function: where d is the depth (in m) and δ is the depth-preference exponent, which equals the depth where W depth = 0.5. At 0 m depth, W depth is kept at a constant value (0.99). W home also follows a sigmoid function: where h is the distance to the home range centre (in km) and η is the home-range exponent, which equals the distance where W home = 0.5. At the home range centre, W home is kept at a constant value (0.99). W schooling is a simple step function of whether other dolphins occur at a patch: where σ is the schooling preference and n is the number of dolphins present. When running a simulation, dolphin's home range centres were calculated at an hourly rate as the centre of all previously visited locations during a spin-up phase of 30 days, after which home range centres remained unchanged and dolphin movements were recorded for 90 days. Afterwards, a summary including home range characteristics and distributions of hourly displacements, group sizes, and depth was created. More detailed information on the model can be found in the appended ODD (Suppl. Methods A); a summary of parameters and their descriptions is given in Table 1.

Sensitivity analysis and calibration
The main model flexibility lies with four parameters that cannot be directly estimated from observations; the movement activity exponent (μ), depth-preference exponent (δ), home-range exponent (η), and schooling preference (σ). Because direct estimation of these parameter values is impossible -the observed patterns emerge through the interaction of these parameters -we determined the model sensitivity to these parameters and calibrated the sensitive parameters to observed patterns of dolphin behaviour. For this process, we (i) reduced the number of parameters to be calibrated based on an estimate of the global sensitivity of the goodness-of-fit from multiple one-at-a-time sensitivity analyses (OAT) through linear regression (Ten Broeke et al., 2016) and (ii) determined the optimal combination of parameter values through regression of the idealised goodness of fit landscape to the results of a full factorial assessment of the model over a wide range of parameter combinations (see Suppl. Methods B; Thiele et al., 2014).
Because we did not have observations that could be directly associated with model outputs, we compared model outputs with observed distributions of four key variables, as shown in Table 2. We measured the degree of mismatch between model outputs and observed distributions using a cost function. Following Thiele et al. (2014), we defined the cost as , , .
, , where i indicates the data category number, j the bin number, E i,j,lower the lower boundary of the expected range, E i,j,upper the upper boundary of the expected range, and O i.j the simulated results within that category and bin. Dolphin behaviour in our model is based on survey data for Maui dolphin (the North Island subspecies) and Hector's dolphin (the species as a whole). Data has been collected on i) depth preferences, ii) move length distributions, iii) grouping behaviour, and iv) home range sizes (Table 2; see Suppl. Methods C).
To compare values between the different cost measures, the values were rescaled, by which a cost of 1 indicates a deviation double the expected value. A cost of zero would indicate an optimal fit of the model to actual observations of Maui dolphin behaviour. Assuming that the cost function follows a quadratic form through the 3-parameter space, we calculated the parameter combination that minimizes the total cost (i.e. maximum of the four costs per parameter combination) using an analytical approach (Suppl. Methods D; Suppl. Fig. 2). By using the maximum of the four costs, we decrease the chance that the best performing parameter combination results in high deviations in one of the key behaviours as delineated in Table 2. Using summed up costs for model calibration resulted in estimated parameter values that Probability distribution of the occurrence of Maui dolphins along the west coast of New Zealand's North Island, according to our calibrated model (movement activity exponent μ = 5.1, depth preference δ = 50, home range exponent η = 27.5, and schooling preference σ = 322; depth cost = 0.028, home range cost = 0.076, hourly displacement cost = 0.105, and group size cost = 0.023). The map represents the total probability distribution of 100 simulations. Colours indicate estimates of K 50 (red), K 95 (orange), and K 99 (yellow). The dashed and blocked areas represent the current gillnet prohibition and trawling and gillnet prohibition zones, respectively. M. de Jager, et al. Ecological Modelling 402 (2019) 59-65 provided qualitatively similar simulated dolphin distributions as when the maximum of the four costs was used for calibration (Suppl. Fig. 3).

Interpretation
Once the model parameters were calibrated, simulations with parameter values that minimized cost were run to generate probability distribution maps of dolphin presence in the study area. To create these maps, dolphin locations were recorded during a single season (summer). To examine the robustness of our results, we repeated simulations 100 times for each likely parameter combination, as well as for parameter combinations at the edges of the simulated parameter space. We estimated the utilization distribution using the adehabitatHR package in R (R version 3.5.1 © 2018 The R Foundation for Statistical Computing) and subsequently calculated the kernel density estimates K 50 , K 95 , and K 99 , which are the areas in which all individuals have resided 50, 95, and 99% of the simulated time, respectively ).

Sensitivity analyses
Before initializing a full calibration of the model, we performed oneat-a-time sensitivity analyses to examine which parameters would influence model output. Analysing the set of one-at-a-time sensitivity analyses, we observed that one of the parameters, the depth preference exponent, did not significantly affect any of the four cost estimates (see Suppl. Results A). Exclusion of this parameter resulted in a more efficient search for the optimal parameter combination in three dimensions.

Calibration
After excluding variation in the depth preference exponent, we performed a full orthogonal sampling of different combinations of the three remaining parameters (i.e. movement activity exponent (μ), home-range exponent (η), and schooling preference (σ). We found the parameter values giving outputs which most closely match the observed distributions shown in Table 2 to be μ = 5.1, η = 27.5, and σ = 322, at Table 2 Observed distribution of four key variables, obtained from field observations as follows. Depth distribution data was obtained from population surveys (Slooten et al., 2004(Slooten et al., , 2006Rayment and du Fresne, 2007). We collected group size distribution data from Webster and Edwards (2008) and Oremus et al. (2012). The home range characteristics have been estimated by  and Stone et al. (2005) and the distribution of hourly displacements was collected from behavioural and photo-identification surveys (Slooten et al., 1992(Slooten et al., , 1993Slooten, 1994; Fig. 2. Maximum cost of the four cost measures averaged per parameter combination. In all panels, the movement activity exponent (μ) is displayed on the x-axis (ranging between 1 and 7) and schooling preference (σ) on the y-axes (ranging between 10 and 1000, shown here logarithmically). The panels differ in home range exponent (η), ranging from η = 10 in the bottom left panel to η = 50 in the upper right panel. A cost of 1 indicates a deviation from the expected value that is twice the expected value of that observed behavioural pattern. δ = 50 ( Fig. 2; Suppl. Results B).

Interpretation
Using these parameter values, we estimated the area used by Maui dolphins within the study site. Fig. 1b illustrates the K 50 (red), K 95 (orange), and K 99 (yellow) of dolphins simulated with 100 separate runs. Our estimates of the K 95 and K 99 lie well beyond the boundaries of the current marine mammal sanctuary (Fig. 1b).
To check the robustness of our model, we also ran 8 × 100 simulations with parameter combinations from the edges of the parameter space (μ = 1 or 7, η = 10 or 50, and σ = 10 or 500; Fig. 3). The badly parameterized model simulations (with μ = 1, top panels Fig. 3) result in a more spread out dolphin distribution compared to the calibrated model.

Discussion
To promote the survival of the critically endangered Maui dolphin population, conservation policies regarding restrictions on fishing and trawling activities need to be fine-tuned on the currently only coarsely known Maui dolphin distribution. We used an individual-based model to estimate the occurrence of Maui dolphins along the west coast of New Zealand's North Island and thereby predict local interactions between Maui dolphins and fishing activities. To provide high quality results, four model parameters had to be calibrated. Our model incorporates small-scale spatial detail in the estimates of the spatial distribution of Maui dolphins that can be used to estimate the effect of human (fishery) activity on population dynamics. By calibrating the model, we increased the credibility of the simulation results, as the model mimics the observed Maui dolphin behaviour best when using the calibrated parameter values. Furthermore, the high spatial resolution of the model can increase the legitimacy of applying models in policy design.
It should be noted that, in its current form, the model does not incorporate any interactions between dolphins and vessels. The dolphin distribution estimated with the calibrated model results in most of the dolphins being found inside the gillnet restricted zone; yet the dolphin distribution extends far beyond the current limits of this area (Fig. 1b). This is consistent with sightings outside the protected area, at least as far as Whanganui, with occasional sightings as far as Wellington and up the east coast of the North Island of New Zealand.
We focus on calibrating our individual-based model to observed real world patterns. It is important to note that taking values directly from field data as parameter input in an IBM will not result in model outcomes resembling observed behaviour patterns, because the different behaviours of individuals in an IBM interact (De Jager et al., 2014). For example, the distance moved per hour by Maui dolphins can be derived from field data (Slooten et al., 1992(Slooten et al., , 1993Slooten, 1994;). Yet, we cannot use these values as input and assume that the model is accurately parameterized. The observed movements are the consequence of interactions between intrinsic movement strategy, attraction by dolphins, depth preference, home range dependencies, and other behaviours. By using such an observed movement pattern as Fig. 3. Probability distribution of the occurrence of Maui dolphins along the west coast of New Zealand's North Island, according to our model using the parameter combinations from the edges of our simulated parameter space (movement activity exponent μ = 1 & 7, depth preference δ = 50, home range exponent η = 10 & 50, and schooling preference σ = 10 & 500). From top left to bottom right, depth cost = 0. 025, 0.030, 0.026, 0.026, 0.024, 0.018, 0.031, and 0.025; home range cost = 0.267, 0.101, 0.267, 0.157, 0.260, 0.184, 0.257, and 0.272; hourly displacement cost = 0.390, 0.387, 0.532, 0.553, 0.385, 0.424, 0.089, and 0.103; and group size cost = 1.236, 1.354, 0.033, 0.020, 1.220, 1.361, 0.021, and 0.050, respectively. The maps represents the total probability distributions of 100 simulations. Colours indicate estimates of K 50 (red), K 95 (orange), and K 99 (yellow). The dashed and blocked areas represent the current gillnet prohibition and trawling and gillnet prohibition zones, respectively.
the intrinsic movement strategy in a model in which interactions are still to come will certainly result in erroneous outcomes (De Jager et al., 2014). Hence, parameter values, which were given at first, but interact with other modelling assumptions, should be incorporated in the calibration process as parameters of unknown value.
Two of the downsides of spatially explicit individual-based modelling are that it is computationally and data intensive (Grimm and Railsback, 2005). Each run of a well thought-out model may require at least several minutes to weeks in extreme cases, while population-level models may take seconds or less. Due to these long running times, parameter estimation and model calibration become challenging. Furthermore, the data that is required to base the model on should be highly detailed at both the spatial and temporal scale (as these would be the resolutions at which behaviour and interactions are modelled in IBMs). Models that provide a good representation of a study system generally include several parameters whose value cannot be determined directly by empirical data. To calibrate a model, results from simulation runs with different combinations of parameter estimates need to be compared to real data (Grimm and Railsback, 2012). Yet, the number of simulation runs required for such an analysis increases exponentially with the number of parameters; hence, calibration remains a modelling process that is often neglected (Thiele et al., 2014), but which should be incorporated in individual-based modelling to generate more accurate predictions.
Different approaches have been used to optimize model calibration, such as manual calibration, statistical calibration techniques using maximum likelihood or Bayesian approaches, or inverse modelling (Sobol, 2001;Wiegand et al., 2003;Grimm and Railsback, 2005;Thiele et al., 2014;Van der Vaart et al., 2015;Ten Broeke et al., 2016). For instance, Wiegand et al. (2003) use a Monte Carlo Filtering method, in which the different cost measurements are successively applied as filters. We used a different approach: we used the highest cost out of the four costs for home range characteristics, hourly displacements, depth preference, and group size distribution to represent the maximum cost per parameter combination. Using this method, the worst fit of the model to the empirical data (indicated by a high cost score) represents a parameter combination's cost, which will result in similar results as the method proposed by Wiegand et al. (2003).
To improve the current model, more survey data is needed, especially on fish distributions. In our model, we assumed a homogeneous distribution of prey, as data on fish distributions is lacking. Yet, it is widely known that fish are not distributed evenly across the ocean, but are rather clustered in patches of high nutritional value. Dolphins subsequently search for such dense aggregations, a behaviour which is currently not incorporated in our model. Information on fish distributions would much improve the current Maui dolphin movement model.
Distribution maps generated with this model will improve the current spatial distribution estimates of the Maui dolphin population and provide key input for comparison of different conservation measures. By overlaying the estimated Maui dolphin distribution map with a map indicating current fishing activities, we can pinpoint the areas where human-dolphin interactions are most likely to cause casualties. Furthermore, for different future scenarios of protection zones, we will be able to estimate the impact on the Maui dolphin population. Calculation of the overlap between Maui dolphin habitat and fishing activity as well as predictions of future scenarios will be important future directions of research.

Declarations of interest
None.