Using Markov chains to quantitatively assess movement patterns of invasive fishes impacted by a carbon dioxide barrier in outdoor ponds

Natural resource managers use barriers to deter the movement of aquatic invasive species. Research and development of new invasive species barriers is often evaluated in pond and field scales using high‐resolution telemetry data. Telemetry data sets can be a rich source of data about fish movement and behavior but can be difficult to analyze due to the size of these data sets as well as their irregular sampling intervals. Due to the challenges, most barrier studies only use summary endpoints, such as barrier passage counts or average (e.g., mean or median) fish distance from the barrier, to describe the data. To examine more fine‐scale fish movement patterns, we developed a first‐order Markov chain. We then used this model to help understand the impacts of a barrier on fish behavior. For our study system, we used data from a previous study examining how bighead and silver carp (two invasive fish species in the United States) responded to a carbon dioxide (CO2) barrier in a pond.

assessments of fish responses. However, tracking of individual fish could allow for a more robust evaluation of other behavioral metrics that better describe how fish respond at finer scales during interaction with a deterrent stimulus.
High-resolution spatial and temporal telemetry data provide an excellent opportunity to characterize patterns of fish movement and behavior that are not possible to determine with conventional summary endpoints. Historically, the utility of telemetry data for evaluating fish behavior in captivity was limited by the size of telemetry units (Baras & Lagardre, 1995). Technological advances now allow high-resolution telemetry data with more precise fish locations that can be used to characterize fish behavior over fine spatial and temporal scales (e.g., Cupp et al., 2017 who had resolution <1 m). Unfortunately, high-resolution data sets can present analytical challenges (Edelhoff, Signer, & Balkenhol, 2016). In particular, the movement of study organisms may not be completely captured despite the large size of the data sets because of gaps due to missing points and erroneous points caused by the data collection process. Furthermore, signals can reflect off hard surfaces and result in two equally plausible points occurring nearly simultaneously. Consequently, data processing can not only be challenging due to the large size of the data sets themselves, but also because of the need to address the problem of multiple plausible locations.
Alternative methods that examine the path of the target organisms rather than the direct coordinate data can mitigate some of these challenges (Cushman, 2010). Edelhoff et al. (2016) described methods to analyze spatial data and trajectories of animals by studying movements and paths detected with tracking systems. The trajectories of animals are composed of several different parts, including direction, velocity, and location (Calenge, 2006;Calenge, Dray, & Royer-Carenzi, 2009). Traditionally, these approaches have been applied to large-scale animal movements (e.g., free-ranging wildlife moving across landscape rather than fish in ponds) and time series with greater time intervals (e.g., hours or days rather than seconds) and regular time intervals between data points (e.g., Langrock, King, & Matthiopoulos, 2012;Patterson, Basson, Bravington, & Gunn, 2009). Although models exist that can analyze continuously collected telemetry data (e.g., Blackwell, 2003), these models present challenges for natural resource managers. The models can be mathematically complex, which places the models beyond the reach of many applied ecologists and resource managers (cf., Ellison & Dennis, 2010) and difficult to implement computationally. In contrast to these previously described models, firstorder Markov chains use linear algebra (Allen, 2011), a topic more likely to be accessible and useable by applied ecologists and resource managers. Furthermore, most spatially and temporally continuous movement models focus on habitat selection for foraging, as opposed to avoidance of adverse conditions (Hussey, Kessel, & Aarestrup, 2015;Kays, Crofoot, Jetz, & Wikelski, 2015;Owen-Smith, Fryxell, & Merrill, 2010).
Markov chains have a few key assumptions as part of their definition that are worth noting for our example, which is a discrete-time, finite state-space Markov chain (Allen, 2011). First, Markov chains assume discrete events (e.g., spatial grid cells) in contrast to modeling continuous events (e.g., exact locations with coordinates). In our example, we modeled fish moving between discrete cells (e.g., cell 1 to cell 2) over a countable number of cells (mathematically, this would be the set of cells, specifically for our example {1, 2, …, 35, 36}) rather than specific, continuous coordinates. Second, Markov chains assume a "lack of memory" between events. That is to say, outcomes at one time step only depend upon the current state, but are independent from previous time steps. In our example, a fish movement at time t does not impact the fish's movement at t + 1, but the fish's state at time t is important. We explored this assumption during our modeling efforts and found a second-order Markov chain (where a fish's movement depended upon its previous movement) produced similar results as a first-order Markov chain. Third, Markov chains assume transition probabilities are constant (homogeneous) through time. In our example, we assumed each fish had its own transition probabilities for each CO 2 -exposure time period, which are described in Section 2. Fourth, Markov chains assume the transition matrix between cells includes every possible outcome. In our model, this is a reasonable assumption for fish movement. However, two spatial components are worth noting. First, fish can have a nontrivial (nonzero) movement probability between two nonadjacent cells. This could happen because fish move farther than one cell during the time interval recorded by the hydroacoustic telemetry gear; the telemetry system and data processing erroneously placed fish, or the telemetry did not capture a fish's use of a cell. Second, adjacent cells might have zero (trivial) movement probabilities. This could happen if a physical barrier separated the two cells or a fish did not move between two cells during the study. The fourth assumption that movement captures all possible outcomes also applies to our model because fish mortality did not occur during the studies generating the data.
We explored this approach by developing a Markov chain model to describe the movement of fish during CO 2 exposure using a high-resolution telemetry data set (Cupp, 2020). As part of this process, we first described fish movement metrics (i.e., distance traveled, directional change, velocity, and acceleration). Second, we developed first-order Markov chain models and applied them to spatially discretized data from the ponds. Third, we used the results of the Markov chain with Moran's I to examine spatial autocorrelation during exposure. Fourth, we evaluated the outputs from the previous three endpoints using regression analysis. Specifically, we asked: (a) Do fish alter their direction of movement (i.e., absolute angle of travel, such as heading north), relative angle of travel (e.g., turning right), distance, velocity, and acceleration in response to changing CO 2 concentrations? (b) Does the likelihood that a fish remains in a location depend on CO 2 exposure? and (c) Does the spatial distribution of movement metrics vary throughout the duration of an experiment? 2 | METHODS

| Telemetry data set
We used fish location data from outdoor pond trials developed to test the effects of a CO 2 barrier on bighead and silver carp described in Cupp et al. (2017). Briefly, Cupp et al. (2017) constructed CO 2 barriers in an outdoor concrete pond (10.0 m long, 4.9 m wide, and 1.2 m deep) at the Upper Midwest Environmental Sciences Center in La Crosse, WI. The pond included a shaded area to attract carp (carp prefer shaded areas) and a CO 2 barrier to keep carp away from the shaded area. Five bighead and five silver carp were placed in each pond during each of the three trials and each fish was monitored using acoustic telemetry (Hydroacoustic Technology Inc. (HTI) model 290 acoustic tag receiver, HTI, Seattle, WA). Fish were individually monitored in space using an acoustic telemetry array of 16 individual hydrophones (HTI model 590 series) for a period of 72 hr with ping rates ranging from 2059 to 2675 ms; the temporal resolution of the data was thus irregular. During each trial, fish were allowed to acclimate to the pond for 24-hr before the CO 2 barrier was activated for 24-hr. Fish were monitored for an additional 24-hr following barrier deactivation. Carbon dioxide concentrations were monitored in each trial using Hatch digital titration (Method 8205; Hatch Instruments Loveland, Colorado). See Cupp et al. (2017) for a detailed description of data collection methodology.
We used the telemetry data set from these pond trials as the basis for understanding behavioral responses to the CO 2 barrier using movement metrics and Markov chains. The telemetry data set consisted of 2,094,720 Universal Transverse Mercator coordinates from 30 individual fish (15 bighead carp and 15 silver carp) with an accuracy of <0.2 m . In some cases, multiple location observations of an individual fish were reported with the same time stamp. We retained all observations because true fish location was undeterminable. For example, sometimes two locations are estimated because reflections cause two equally plausible points to be estimated. In general, the undeterminable points occurred because of both physical limitations of the study system (e.g., the reflections described above) and software limitations of the algorithms used to estimate points.
We classified observations from each of the 72-hr trials into five unique time periods before the analysis to test how fish behavior changed in response to the CO 2 barrier. In their final analysis, Cupp et al. (2017) classified the first 24-hr period as the pre-CO 2 period, the middle 24-hr as the CO 2 -exposure period, and the last 24-hr as the post-CO 2 period. However, a concentration ramp-up period occurred between the activation of CO 2 barrier and reaching equilibrium CO 2 concentrations. Similarly, a dissipation period occurred after barrier deactivation. We included a 6-hr ramp-up and a 6-hr dissipation period which were based upon measured concentrations of CO 2 in pond water. Thus, we used five time periods for analyses: the 24-hr period before CO 2 additions ("Pre-CO 2 "), a 6-hr period of increasing CO 2 concentrations ("Increasing CO 2 "), the 18-hr exposure period where CO 2 concentrations were maintained at a constant level ("During CO 2 "), the 6-hr period of decreasing CO 2 concentrations ("Decreasing CO 2 "), and the 188-hr period following the CO 2 additions ("Post-CO 2 ").

| Fish movement metrics
We calculated distance traveled, directional change, velocity, and acceleration using pairs of sequential location observations to describe the movement patterns of each fish. We calculated distance traveled in meters (m) between each sequential pair of location observations using the adehabitat package (Calenge, 2006; see our Statistical analysis and numerical methods section for version details). We characterized directional movement using absolute and relative angles. The absolute angle was defined between the next step and the baseline direction (Marsh & Jones, 1988) and the relative angle was defined as the angle between the next step and current direction, which is often called the "turning angle" (Marsh & Jones, 1988;Root & Kareiva, 1984). We used the relation between distance traveled and change in time between observations to calculate velocity (m s −1 ) and acceleration (m s −2 ).

| Markov chain model
We modeled fish movement using first-order Markov chain models applied to a spatially discretized representation of a study pond to overcome potential location errors related to multiple detections at a given time stamp. We discretized the study pond into 36 cells. Each cell was 11.2-m wide × 1.11-m long and distributed in a 4 × 9 grid. We used numerical experimentation to evaluate the number of grid cells and sought to balance the number of observations per cell versus the total number of cells and led to our choice of 36 grid cells. Furthermore, this created grid cells that were comparably sized to the 0.2-m resolution of the telemetry system . The 36 grid cells had a mean of 1,910 observations per fish, a median of 1,423, and a minimum of 36. Probabilities of fish movement into and out of cells were calculated using paired sequential location observations for each fish. We estimated the probability of fish movement between cells using the empirical movement data for fish (i.e., counting how many times a fish transitioned from one cell to another cell and then dividing by the time interval for that period). This was done for each fish for each exposure period. This allowed us to compare not only how likely a fish was to change cells during each period, but also identify the cells to which a fish was most likely to change during each period. We used violin plots and summary statistics to characterize the distributions of fish movement probabilities from one cell to another during each time period.
Empirical movement probabilities were then used to parameterize first-order Markov chain models (Allen, 2011) for each fish. Specifically, the parameters for the first-order Markov chain were the transition probabilities between cells. We ran each Markov chain model using 6,000 simulation steps. We qualitatively evaluated the agreement between empirical movement measures and those predicted by the Markov chain model. We did this because agreement between the model and the data used to parameterize a model does not necessarily provide strong support for a model, but a lack of agreement would raise concern about a model's structure or assumptions.

| Moran's I
We calculated Moran's I (Moran, 1950) for each mean movement metric to assess the degree of spatial autocorrelation across the five time periods. Moran's I describes the degree of spatial autocovariance in relation to the variance of the data. Moran's I thus has utility in evaluating whether fish exhibit a behavioral response in only certain areas of the pond (e.g., nearest the CO 2 barrier) during any given time period versus little or no spatial patterning of behavioral responses. Moran's I values vary from negative one to one. Negative Moran's I values indicate more spatial dispersion than expected by random chance. Positive Moran's I values indicate more spatial autocorrelation than expected by random chance and thus more "clumpiness" in space. Values approaching zero indicate no spatial correlation. The layout of squares on a chessboard can be used to describe Moran's I. A chessboard of a Moran's I of negative one would be perfectly dispersed, like the black and white squares on a normal chessboard (for an illustrate example of different boards, see fig. 1 In this equation, n is the number of samples, y an observed value, w ij the distance weight, and S w = i n j n ij 0 =1… =1… ∑ ∑ . In this equation, weights for each cell were the cell's Markov chain value. In doing so, we assumed that our empirically driven models of fish movement among grid cells better represented the spatial relations among movement metrics than other weighting schemes (e.g., inverse distance). Our final result was a set of Moran's I values describing the spatial autocorrelation of each movement metric in the discretized pond during a given time period within a trial for each fish.

| Statistical analyses and numerical methods
We developed a series of mixed-effects models using (a) movement metrics, (b) movement probabilities, and (c) Moran's I value to assess whether fish exhibited behavioral response to the CO 2 barrier. First, we tested whether each movement metric varied in response to the CO 2 barrier treatment using a mixed-effects model with "period" and "species" being fixed effects and "trial" being a random effect. Our models used the 75th quantile of each movement metric for each grid cell for each individual as response variables, allowing a comparison of the tail of the distribution rather than the central tendency due to highly skewed distributions of movement metrics. For example, distance traveled had skewness values that ranged from 3.8 to 4.8 (see Figure 1 for an example of distance traveled). This skew occurred because fish movement has a minimum distance traveled of zero, whereas the observed maximum distance traveled was much greater than the mean distance traveled. The 75th quantile was chosen because the distributions had similar medians, but high right skews in their distributions. This approach is similar to quantile F I G U R E 1 Example distributions of data showing distance traveled (m) for a fish during different periods of CO 2 treatments. Boxplots have the median line in the middle of the box. The box is the middle 50% of the data and the stems are the upper and lower quarters of the data. Points are outliers. The purpose of this plot is to demonstrate how data have a positive skew. BHC are bighead carp and SVC are silver carp regression (e.g., Cade & Noon, 2003), but precalculates the quantiles before fitting the statistical model. This precalculation solves three problems. First, traditional quantile regression approaches would likely have produced statistically significant results regardless of biological differences due to large sample sizes. Second, traditional quantile regression can be computationally expensive on small data sets and difficult if not impossible with our data set due to the large sample size observations. Third, using the 75th quantile avoids pseudo-replication at the fish level. For our second mixed-effects model, we predicted the empirical probability of a fish not leaving a cell based on the results of the Markov chains with "trial" as a random effect and "period" and "species" as fixed effects. Finally, we developed a mixed-effect model to predict Moran's I values with "trial" being a random effect and "period" and "species" being fixed effects to formally test for differences in the spatial autocorrelation of movement metrics among time periods.

| Fish movement metrics
For all movement metrics, mixed-effects model coefficients for pre-and postexposure movement metrics were not different across trial time periods, indicating fish movement patterns were similar before and after CO 2 barrier activation ( Figure 2). Specifically, the postexposure coefficients were all near zero and had 95% CIs that included zero. However, fish movement varied in response to CO 2 exposure during the experimental treatment: carp traveled farther during increasing CO 2 -exposure periods (0.21 m, 95% CI 0.177-0.277), full CO 2 -exposure periods (0.133 m, 95% CI 0.08-0.188), and decreasing CO 2 -exposure periods (0.067 m to 0.013-0.110). The absolute angle varied in response to CO 2 exposure during increasing CO 2 -exposure periods (−1.64, 95% CI −0.21 to −0.12) and during decreasing CO 2 exposure (−1.09, 95% CI −0.142 to −0.077). Likewise, the relative angle varied in response to CO 2 exposure during increasing CO 2 -exposure periods (−0.244, 95% CI −0.28 to −0.199) and decreasing CO 2 exposure (−0.144, 95% CI −0.199 to −0.099). No upper 75th quantiles for acceleration differed from the pre-exposure period (all estimates were near zero and the 95% CI coefficients include 0; Figure 2).

| Markov chain model
Simulations of fish movement probabilities using the Markov chain qualitatively matched empirical data used to parameterize the model (Figure 3). The transition probabilities for staying in the same cell (i.e., the mean of the diagonal of the Markov chain) differed across time periods for both carp species (Figures 3 and 4). Fish stayed in their same cell during 54.2% (95% CI 50.9%-57.5%) of the time during pre-exposure and 55% (95% CI 51.9%-58.5%) of the time during postexposure. In contrast, fish only stayed in their cell 44.9% (95% CI 41.5%-48.2%) of the time during increasing CO 2 exposure. When CO 2 reached equilibrium, fish stayed in their cell 51.2% (95% CI 48.5%-55.2%) of the time. During the decreasing CO 2 period, fish stayed in their cell 49.3% (95% CI 46.0%-52.7%) of the time.
F I G U R E 2 Regression coefficient estimates (solid black points) and 95% confidence intervals (horizontal black lines) for four mixed-effects models predicting absolute relative angle, acceleration, distance, and relative angle. Coefficient estimates include a reference intercept (pre-CO 2 exposure for bighead carp, "Intercept"), silver carp coefficient ("silver carp"), and four other CO 2 -exposure periods ("Increasing CO 2 ," "During CO 2 ," "Decreasing CO 2 ," and "Post-CO 2 "). The vertical red line indicates zero F I G U R E 3 The spatial distributions of observed and simulated occurrence probabilities of bighead carp and silver carp during five CO 2 -exposure periods. Occurrence probability is the probability a fish used a cell and was calculated across all individual fish. For example, if a cell had an occurrence probability of 0.10, fish spent a mean of 10% of the time that specific cell. Within each rectangular pond (size of 4 × 9 cells), the white line indicates a physical barrier that creates a separation in the upper portion of the pond, and the orange horizontal bar indicates the position of the CO 2 barrier F I G U R E 4 Violin plots for the mean probability that a fish stayed in the same cell during each exposure period. The widths of violin plots reflect a kernel density estimate of observations (mirrored for symmetry), with median values of the distributions indicated by horizontal bars. Points within each violin plots are the means of individual fish. Changing CO 2 concentrations (i.e., time periods of increasing and decreasing CO 2 ) and elevated CO 2 during exposure increased the variability of individual fish staying in the cells. This demonstrates individual fish respond differently to CO 2 in the water 3.3 | Spatial autocorrelation of fish movement metrics All movement metrics became less spatially correlated during periods of CO 2 exposure as indicated by lower Moran's I values (Figures 5 and 6) and substantiated through regression analysis. The three distance-based metrics (distance, velocity, and acceleration) and both angle-based metrics (relative and absolute angles) had similar patterns of decreasing spatial correlations. The three distance-based metrics all had similar baseline Moran's I values with estimates between 0.55 and 0.68 across species and endpoints (lowest 95% CI 0.45 and highest 95% CI 0.75). Likewise, the exposure values for these endpoints were lower and had similar values to each other, ranging from 0.37 to 0.44 (lowest 95% CI 0.26 and highest 95% CI 0.52). Relative angles had a baseline Moran's I of 0.5 for both species (silver carp 95% CI 0.6-0.54; bighead carp 95% CI 0.39-0.54) and exposure values of 0.24 for silver carp (95% CI 0.15-0.33) and 0.37 for bighead carp (95% CI 0.29-0.45); absolute angles had a baseline Moran's I of 0.6 (silver carp 95% CI 0.54-0.68; bighead carp 95% 0.59-0.70) and exposure values of 0.46 for silver carp (95% CI 0.39-0.53) and 0.56 for bighead carp (95% CI 0.50-0.62).

| DISCUSSION
We documented how two invasive carp species changed their behavior in complex ways in response to a CO 2 barrier. This supports previous work showing that CO 2 can be an effective deterrent to the spread of invasive fish populations (e.g., Cupp et al., 2017). For example, our Figure 3 shows the carp moved as far away from the CO 2 barrier as possible. Our results are also in concordance with work by others that found molecular-level responses to CO 2 in Asian carp F I G U R E 6 Regression coefficients estimated for Moran's I mixed-effects models. Points indicate the coefficient estimate and horizontal lines correspond to the 95% confidence intervals. The vertical red line is zero. The intercept corresponds to the pre-CO 2 -exposure period. For all movement metrics, the average fish had lower Moran's I estimates during the increasing, elevated, and decreasing CO 2 -exposure periods. This change was similar for all three CO 2 -exposure periods. This demonstrates CO 2 all three time periods had similar disruptions on fish movement behavior and effects on Asian carp behavior (Dennis, Adhikari, & Suski, 2015;Dennis, Kates, Noatch, & Suski, 2015). Furthermore, our approach captured more complex responses than typically assessed in fish movement studies and has the potential to provide insights into fish behavior and response to stressors. Both bighead and silver carp altered their distance traveled, direction, velocity, and direction angles when the CO 2 barrier was initiated and returned to original movement patterns once CO 2 exposure ended, suggesting that the exposure did not have residual effects on behavior. In addition, fish were more likely to move during periods of increasing CO 2 compared to pre-and post-CO 2-exposure periods. Specifically, the fish moved as far away from the CO 2 injection area as possible during exposure (Figure 3), which was also captured by the change in absolute angle of travel (i.e., fish changed their movement direction to avoid CO 2 ). The spatial patterns of the movement metrics changed throughout the trials: Movement metrics were highly clumped in the pond during pre-and post-CO 2 exposure, but were more spatially disaggregated during periods of CO 2 barrier activation. These patterns may be influenced by the spatial distributions of fish themselves within the pond area but also may reflect underlying spatial gradients of altered fish behavior. That is, fish may exhibit a behavioral response during CO 2 exposure only in certain areas of the pond-nearest the CO 2 barrier-while other areas may not elicit a response, resulting in less spatial autocorrelation in distance and angle metrics. Spatial analyses designed to discriminate among such possibilities were beyond the scope of this study but could be fruitful endeavors for future work. Together, our results agreed with previous findings that CO 2 barriers are avoided by Asian carp (Dennis, Adhikari et al., 2015;Donaldson, Amberg, & Adhikari, 2016) but provide much more detailed spatial and temporal information on their avoidance behavior.
The quantitative movement metrics and numerical models we used in this study provide greater insight into the mechanics of CO 2 barrier responses than location as an endpoint for behavior and avoidance alone. The variation in both the absolute and relative angles of carp movement, as well as decreased velocities and shorter distances traveled during exposure periods, reflect how CO 2 exposure changes the basic movements of carp. The smaller angles and shorter distances of fish movement may be because carp are limited to a smaller area in the pond as they avoid exposure to CO 2 , which can cause physiological stress (Eddy, Lomholt, Weber, & Johansen, 1977). In natural settings, carp would likely be able to simply swim away from the barrier, but the more variable swimming behavior immediately following exposure indicate that CO 2 may temporarily stress and disorient the fish. However, our analysis did not allow us to tease apart if fish were seeking out areas with lower CO 2 or had their behavior affected by CO 2 . Although beyond the scope of our work, multiple possible approaches could be used to compare these possible reasons. First, the movement of individual fish could be modeled and examined through time using Markov chain state-space models on the level of the individual fish. Second, studies similar to Cupp et al. (2017) have been conducted examining other stressors (e.g., sound by Murchy et al., 2017). Data sets, such as Murchy et al. (2017), could be re-analyzed that would allow comparisons between stressors. The responses to these other stressors would provide empirical support for why fish movement changes. We also observed that fish were less likely to stay in the same grid cell during periods of CO 2 exposure, indicating the fish were moving around more within the ponds, likely to find areas with lower CO 2 concentrations. Our findings that movement metrics became less spatially clustered during exposure periods support these interpretations as well. Whether these types of movements occur in response to barriers in the field is unknown, but these types of metrics could be used to quantify effective deterrence or in studies evaluating stress responses.