A logistic function to track time-dependent fish population dynamics

This paper uses a two-parameter logistic function to model the dynamics of length-at-maturation for the Barents Sea capelin over the past 47 years. We estimate the function parameters using a combination of length-age data from scientific surveys, and commercial catch statistics. Using temporal variability in the function parameters, we demonstrate that the time series of stock biomass defines a three-state Markov process, that qualitatively represent high, moderate, and collapse states of the stock biomass. We make inference about transition times between the states by calculating the mean passage times for the Markov process. Our analyses also show that maturation intensity is higher at low stock size (leading to shorter lengths at maturation), compared to when biomass levels are either high or moderately high. Our results are central to management of this stock, as uncertainty in estimating the proportion of maturing biomass affects harvest decisions and ultimately, the sustainability of the stock.


Introduction
Maturation rates are key drivers of fish stock population dynamics, as they are intrinsically linked to the reproduction potential of the stock (see e.g., Trippel, 1995). This link is also important from a management perspective, as most biological reference points or harvest rules (e.g., spawner escapement) are defined in terms of the amount of potential spawners (see Clark, 1991;Hannah et al., 2009;Murawski et al., 2001;Vitale et al., 2006). Age and mean body weight at maturity can also be used to predict the risk of overexploitation (Reynolds et al., 2005). Hence determining how maturation may depend on age and/or size (e. g., length or weight) is of ecological importance, and central to the sustainable management of most marine species.
In principle, the maturation stage of individuals may be determined (for most species) by direct inspection of gonads during periods preceding the spawning period (Flores et al., 2015;Williams and Babcock, 2005;Smith and Walker, 2004;Berglund, 1992). This is however, not always feasible, especially for species where there is a large time-span between when scientific surveys are conducted, and the onset of maturation. For such cases, projections of maturation scheduling may be estimated by linking the maturation process to some physiological species property obtained during the scientific survey. The ability to link growth, especially in length, to maturation is attractive for several reasons. For most species, length measurements are most reliable, as body length increases monotonically with age, in contrast to e.g., weight. Furthermore, length is perhaps the easiest physiological measurements that can be obtained at reasonably low operational cost, even for large population samples. However, maturation rates may change in response to variability in the marine environmental, such as, the availability of food (Rilling and Houde, 1999), and abiotic factors, such as, ambient temperature (Sumpter, 1992), salinity and ocean acidity (Brett, 1979). Hence establishing a link between maturation and a physiological property as length may be challenging. This challenge may be addressed by the use of mathematical/statistical growth models.
Mathematical/statistical growth models quantify the increase in body dimensions over time. The literature catalogs several models, of which the four most common are the Logistic (Verhulst, 1845;Schnute, 1981), von Bertallanfy (von Bertalanffy, 1938), Gompertz (Gompertz, 1825) and Richards (Ricker, 1979) functions. Though structurally different, the models share a common characteristic of being monotone increasing functions of age. In general, growth models also provide parsimonious tools for linking both biotic and abiotic factors to variability in the growth process. The logistic function, for instance, has been employed in modeling growth in bacterial cultures (McKendrick and Pai, 1912), birds (Tjørve and Tjørve, 2017), and trees (Zeide, 1993), and demonstrated to give good fit to fish weight and length data (Katsanevakis and Maravelias, 2008;Lugert et al., 2016;Bennett et al., 2007;Llamas and Ratkowsky, 2004;Zhu et al., 2016).
Since it is only mature fish that spawn, the assumption that fecundity in fish is size dependent (Tsoukali et al., 2016) establishes a link between maturation and size. Hence, changes in growth rates are assumed to affect maturation schedules (Hunter et al., 2015). For species where maturation is considered a function of length, the two-parameter Logistic function has proven to be a useful model for fish maturation. For instance, the two-parameter function was used by Freitas et al. (2016) to identify the attainment of maturity in catfish, and by Chen and Paloheimo (1994) to estimate the fish length and age at 50% maturity, while Bogstad and Tjelmeland (1990) (see also Gjøsaeter et al., 2002a) used it to estimate the annual proportion of the stock that is mature, based on scientific survey data.
For practical and operational reasons the capelin stock size is estimated by an annual scientific survey (conducted since 1972), that is conducted six month ahead of the spawning time. This survey gives estimates of stock size (numbers and biomass) at October 1. An anticipatory Total Allowable Catch (TAC) is set so that there is a 95% probability of the mature stock biomass is above 200 kilotonnes by April 1 (assumed spawning time), the following year. Hence determining the proportion of the stock that is maturing, is at the heart of the management decision.
Given the early timing of the scientific survey, gonad inspection to determine the proportion of potential spawners, will be highly uncertain. On the other hand it has been found that dependence of maturation on length, rather than age, is a more stable relationship (Tjelmeland and Bogstad, 1993;Forberg and Tjelmeland, 1983), with faster growing capelin having higher propensity for early maturation than those with slower length growth rates (Hopkins and Nilssen, 1991;Forberg and Tjelmeland, 1983). Hence the survey data is used to split the Barents Sea capelin stock into maturing and immature proportions, based on body length (Tjelmeland and Bogstad, 1993). A two-parameter logistic function (which is described shortly) is used to describe the maturation schedule for different age-classes.
We hypothesize that temporal variability in the logistic function parameters can be used to track time-dependent changes in the population dynamics of the Barents Sea capelin. We investigate this by analyzing parameters derived from fitting the two-parameter logistic function to observation data from scientific surveys. The capelin stock has also been prone to several episodic collapses, occurring in 1985-1989, 1994-1997-2007and recently in 2017(Gjøsaeter et al., 2002aWGIBAR, ICES, 2019). We investigate whether the logistic model parameters may be used to detect changes in population dynamics e.g., transitions from collapse to stable states, or vice-versa.

The data
Since 1972, annual joint scientific surveys by the Institute of Marine Research (IMR, Norway) and the Polar Research Institute of Marine Fisheries and Oceanography (PINRO, Russia) have collected data on capelin in the Barents Sea. In this paper, we use data consisting of estimated population numbers (specified by age and length) from the joint scientific surveys, and commercial catch data reported for fish of length less than 14cm in different seasons. Fig. 1a-b shows (on log-scale) survey estimates of capelin population size, and reported annual commercial catch of ages 2-4 capelin with length less than 14 cm, during 1972-2018. All data used in this manuscript have been obtained from the database of the ICES Working Group on the Integrated Assessments of the Barents Sea (WGIBAR, ICES, 2019).

The maturation model
We follow (Gjøsaeter et al., 2002a) in modeling r as a logistic function of two parameters p 1 and p 2 , as: where p 1 (positive number) is the maturation intensity, p 2 is the median (50%) length at maturity, and μ l is the mid-point of length-class l (Bogstad and Tjelmeland, 1990). If we define the maturation function more generally as in (2), such that Hence, by setting α = 4 in (3), we can more easily interpret p 1 as the change in maturation proportion when l = p 2 . Fig. 2 is an illustration of how r varies with p 1 , and p 2 . Let a, l and y represent age, length-class, and year respectively. Define N y,l,a as the number of fish estimated by the acoustic survey, and r as the proportion of maturing fish in a given length-class. Given r(l, p 1 , p 2 ), (4) gives the number of immature fish, The goal is to determine the parameters p 1 , p 2 , and a constant mortality term p 3 , using N y,l,a data, over a number of years.

Prediction and parameter estimation
The estimation of p 1 and p 2 is based on a likelihood function built on comparing the immature capelin projected one year ahead (prediction) to the total stock measured in that year (data). The prediction process is started at the first of October of the year y, taking N y,a (p 1 , p 2 ) as input by considering monthly catch and natural mortality (p 3 ), and finishes at the first of October of the year y + 1 when the total stock size is re-estimated by the acoustic survey, N y+1,a . The prediction procedure only considers the immature component of the stock, which can be compared to the stock statistic in the following year. The mature component dies after spawning during the Spring. As will be shown later in (5), the estimation of p 3 is confounded with the estimation of p 1 and p 2 .
Define C y,a,k as reported catches per month k of fishes with less than 14 cm body length for a given year y and age-class a. Then (5) defines Ñ y,a (p 1 , p 2 , p 3 ), the total population stock size resulting from one year ahead projection of the immature capelin calculated from (4).
Eq. (5) can be obtained from (6) (see Pope, 1972) where (6) is applied for monthly projections, We adopt the Pope's approximation to the Baranov equations, (Pope, 1972) in the equations that model these processes. In particular, the Pope's approximation assumes pulsed harvesting, where fishing takes place as a single event in the middle of each month.
The predicted stock size Ñ y,a in Eq. (5), represents the expectation of the true population stock size N y,a = ∑ l N y,l,a . Following Gjøsaeter et al. (2002a), we let N y,a follow a Gamma-distribution with unknown variance Ñ 2 y,a ν , i.e., using a mean-variance parameterization, Finally we determine the parameter set Π = {p 1 , p 2 , p 3 , ν} by minimizing the negative log-likelihood function We estimate the parameter set Π for two different sets of numerical predictions. In the first set, we consider year-to-year variation in the maturation process, which leads to the estimation of 47 parameters based on (8). The second set of numerical predictions involves the estimation of a single parameter set, Π, which represents average maturation over a period of several years.

Constraints
A combination of stock size values lower than catches over a given period could result in a negative prediction of stock size in equation (5). We follow the approach in (Gjøsaeter et al., 2002a), by setting , ifÑ y,a < 0.5.
so that the stock size is assigned a value in the interval [0.25 0.5], when the predicted stock abundance is negative. We also impose the following constraints on the parameters, p 1 and p 2 , in order to avoid computationally plausible, but biologically unrealistic estimates. For instance, for p 1 ≤ 0.06 results in a non-sigmoidal function (see Fig. 2a), while the constraint on p 2 is derived from empirical data on capelin.
We use the TMB (Kristensen et al., 2016) package (version 1.7.15) in R (version 3.6.0) for minimization of the objective function in Eq. (8).

Parameter-based inference
Firstly, we investigate whether the parameter tuple (p 1 , p 2 ) define distinct clusters in the p 1 : p 2 plane, and whether the clusters can be linked to specific states (e.g., moderate, high, collapsed stock) of the population dynamics. Next, we investigate whether these states are Markovian, i.e., a current population is predetermined only by the state preceding it. We use the K-Mean clustering algorithm to investigate clustering in the p 1 : p 2 plane, and Markov process analysis to address the transitions between states. A brief discussion of the methods K-Mean clustering and Markov analysis are presented in the supplementary material.  Variation of r(l, p 1 , p 2 ), with (a) fixed p 2 , and (b) p 1 , where μ l is the length in cm, p 2 is the length at 50% maturity and p 1 describes the increase in maturity by length at p 2 .

Results
Though our analysis covers capelin of ages 2 and 3, here, we present result for age-3 capelin. This is for the sake of brevity and because this age-group gave the most consistent results. Our results for age 2 capelin are presented in the supplementary material.
We present the results from two sets of numerical predictions. The first set of predictions (annual predictions) consider year-to-year variations in the maturation process, and therefore, in the maturation function parameters. In the second set of experiments (periods of fishing moratorium) we considered the average maturation dynamics over time interval of several years. A detailed discussion and comparison of the results is presented later.

Year-to-year predictions
Time trends of the maturity parameters p 1 , p 2 and the mortality parameter p 3 derived from annual predictions are presented in Table 1 and Fig. 3, which show the high variability in maturation and natural mortality typically associated with the Barents Sea capelin dynamics. We identify four spikes in Fig. 3a, for which p 1 > 0.6. We investigate the link between these spikes and the population dynamics by plotting moving average p 1 and the capelin biomass in the same graph (see Fig. 4a).
The results show that the spikes are coincidental with years of relatively low Capelin stock biomass, and p 1 is inversely related to the capelin biomass. Fig. 4b shows the 3-year moving average for p 2 , whose trend appears to be synchronous with the capelin biomass, for all the years considered. For the natural mortality, we observe synchrony between p 3 and stock size (see Fig. 4c).

K-mean clustering
Most of the maturation curves resulting from substituting estimated parameters p 1 and p 2 in the logistic model (1) have smooth monotone increasing (from zero to one) trajectories. However, variations exist in terms of steepness of the maturation functions especially at lengths close to p 2 . We investigated the variability in the maturation functions by determining whether the tuple (p 1 , p 2 ) define distinct classes of maturation ogives. We use the K-Mean (unsupervised) clustering algorithm to partition the maturation curves (based on (p 1 , p 2 )) into k distinct number of groups. A central step for such clustering algorithms is to determine the optimal number of k clusters into which the data may be clustered. We use the Elbow Method (Thorndike, 1953) (a popular heuristic) to determine this optimal value of k. The method is akin to the mathematical heuristic in optimization, where the search algorithm is terminated when no further improvement can be observed in the objective function. For the Elbow method, starting from a very low value of k (e.g., k = 1) one calculates the within-cluster sum of squares (wss) for increasing values of k. If no improvement is observed in the wss for k≥ k *, then k* may be chosen as the optimal number of clusters. A graphical approach (plotting k against wss) offers an easy way of determining k*, which is located at the bend (elbow) in the plot. Fig. 5 shows the result obtained using the Elbow method, which suggests k = 3 as optimal number of clusters. Fig. 6a-d show results from allocating the tuples (p 1 , p 2 ) to k = 3 clusters. Notably, group 1 consists of maturation functions with large p 2 values, while group 2 functions have steep gradients that correspond to larger maturation intensities.

Markov properties
We use the clustered groups of maturation functions to assign biomass states to the time series of capelin biomass. For each year, we identify the group to which the maturation function belongs, and derive the state from the cluster group number, see Fig. 7. We deduce that state 1 is a relatively high biomass indicator, while state 2 is mainly associated with years of stock collapse. Moderate biomass levels are associated with state 3, which appears to persist, when it first occurs. We analyze the sequence of identified biomass states to determine whether they are Markovian, and if so, derive relevant Markovian properties that characterize the biomass dynamics. We tested for Markovianness of the sequence using the method presented in Hart and Martínez (2014), for which a brief description is given in the supplementary material. Our result shows that we cannot reject the null hypothesis (the sequence is Markovian) for 99 of 100 replications. Next, we derived the transition  and stationary probabilities, as well as the mean passage times, see Fig. 8. The relevance of quantifying mean passage times, is to determine e.g., how long it will take for the stock in a collapse state, to reach recovery.

Periods of fishing moratorium
The capelin stock has been through multiple episodes of fishing moratorium, in response to stock collapse (biomass <100kt, (see Gjøsaeter et al., 2015)). We used all data from 1972-2019, to investigate   the maturation dynamics during periods of fishing moratorium. We do so by estimating the maturation curve per scenario, for two different data scenarios, namely, for years with (a) moratorium on commercial fisheries (bold years in Table 1), and (b) commercial fisheries (black color years in Table 1). Table 2 summarizes the parameter (p 1 , p 2 and p 3 ) estimation results. The table shows the average of individual annual estimates (Avg. mean), which is presented in Section 3.1, as well as the case where all data  were used in estimating the set of parameters (Estimated). The results show that the two approaches in estimating the maturation parameters give identical results. Fig. 9a shows the derived maturation functions for the two scenarios. This figure shows that there is a higher propensity for early length at maturation during years of commercial fisheries (r b , red), than for years when there is moratorium on fisheries (r a , black). We also investigated whether the two maturation functions scale, i.e., whether at any given length, there exists a scalar or functional relationship between the maturation rates for the two scenarios. We do this by calculating the ratio r ab = r a /r b . We estimated the parameters α and β such that, where (p 1 ,p 2 ), and (p ′ 1 , p ′ 2 ) are estimated parameters, in the second and fourth columns respectively, of Table 2. Fig. 9b shows results for the scaling function, using derived parameters α = 1.5313 and β = 13.4701.

Discussion and conclusions
In general, ages 2 and 3 capelin dominate the stock biomass for any given year and while data on age-4 is highly variable, data on age-5 is mostly zero's and not informative. Our analysis found comparable results for ages 2 and 3, though those for age-3 capelin were more consistent and therefore, easier to interpret. This is consistent with our expectations, as age-3 capelin is dominated by maturing capelin during the autumn, when the scientific survey is conducted.
The result in Fig. 4a-b appear to show consistent patterns (inverse, and synchronous relationships, respectively for p 1 and p 2 ) from 1972-2018. However, in Fig. 4c, the relation between stock size and natural mortality (p 3 ) may be divided into two different time periods. While stock size and natural mortality appear to be synchronous for 1990-2018, and inverse relationship is observed for 1972-1990. The 1990-2018 results meet our expectation of natural mortality (which, here includes predation mortality) being proportional to stock size. The reasons for the discrepancy observed for 1972-1990 may be multiple, confounding, and hence, difficult to disentangle. Underlying these are uncertainties associated with data collection protocols, assessment of stock size, and even management. We find the discussion to be beyond the scope of this manuscript, hence refer the reader to (Gjøsaeter et al., 2002a;Gjøsaeter et al., 2015;Tjelmeland, 2005).
Results from the clustering procedure (for p 1 and p 2 in Fig. 6, and Table 2) show that fish reached maturation at shorter median lengths (lower p 2 ) and lower maturation intensity (lower p 1 ) during nonmoratorium years, than otherwise. For the three states identified, the largest range for median length at maturation is associated with years of moratorium (stock collapse). It is plausible that the amount of prey per fish is higher when the stock is low. This may be a driver for growth spurts (intensity spikes), that result in early maturation at length. This explanation is supported by (Ingvaldsen and Gjøsaeter, 2013), which documents density-dependent (stock abundance) effects on distribution, feeding, and growth of capelin. There are however, other factors (e.g., degree of overlap with prey, temperature, ice-cover, see (Ingvaldsen and Gjøsaeter, 2013)) that may drive variability in the growth rates and maturation schedules during periods of fisheries closure. In a semelparous population, a non-selective terminal fishery should not generate strong selection for changed age and size in maturation. Baulier et al. (2012) found this theory to be valid by applying the probabilistic maturation reaction norm method to empirical data.
Of the results from the Markov process analysis, Fig. 8b-c are especially interesting, as they show the long-term propensity of stock transitions from/to different states (Fig. 8b), and the minimum time it takes to re-encounter a current state, or transition to another (Fig. 8c). The highest propensity (0.48) is observed for transition from collapse-or high-, to an intermediate biomass state. In spite of the observed episodic collapses, our analysis shows that there is a low propensity (0.17) to move from a high/moderate stock-to a collapse state, and an equal propensity to remain in this state, once this is encountered. We deduce from Fig. 8c that the minimum time it takes for the stock to transition into collapse, stating from a moderate or high biomass state is approximately 7 years. The time lapse between two consecutive episodes of collapse is approximately 6 years, consistent with observed episodes capelin collapse being reported during the periods 1985-1989, 1993-1997, 2003-2007(Gjøsaeter et al., 2015. The periodicity of collapse (6 years) is also consistent with (Solvang et al., 2018), and the 6.2 year cycle reported in (Yndestad and Stene, 2002) which, according to the authors, is an optimal cycle for natural environmental adaptation, and optimal strategy for growth and survival of the capelin stock. Despite their consistency, these results must be interpreted with a caveat, since this is a population that is subject to commercial fisheries. We can therefore, not exclude the possibility that our inference may be influenced by a combination of intrinsic stock dynamics and management harvest decisions.
The very low estimates of p 3 in Table 2 during the moratorium period may lead to the conclusion that absence of fishing mortality (see (5)) will result in unrealistic longevity of each cohort. This conclusion however, is erroneous. Recall that p 3 is the mortality of the immature stock from October in year y to January in year (y + 1), prior to commencement of the fisheries. Each cohort, even in the absence of fisheries, will experience 100% natural mortality after 4 years since the capelin stock is semelparous, and the longevity of a single fish seldom exceeds 4 years. Another plausible reason is that during moratorium periods, stock sizes are very low, which may imply less competition for food and therefore, high probabilities for survival.
There is need to understand the causes of large variability in the median length at maturation, and drivers of high maturation intensity during years when there is a moratorium on fisheries. Forberg and Tjelmeland (1983) used the two-parameter logistic model (same as in this manuscript) to study the median length at maturation (p 2 ) for 3-year old capelin in different areas of the Barents sea. Significant variations were found, which were linked to variability in physical and feeding conditions. We are unable to draw inference on spatial scales, as this is beyond the scope of this manuscript. A potential area of research is to investigate whether maturation functions from different spatial locations may scale with environmental/physical variables, i.e., whether for two spatial regions A and B, we can write where the "~" symbolizes the existence of a relationship (most plausibly non-linear) between the two functions. In (12), Ω A (similarly for B) denotes any spatial co-ordinate in A, and f is a function (discrete or continuous) that represents the environmental variable(s) (e.g., food density, temperature) at Ω A . This will be a natural extension of our results in Fig. 9. Note, however, that the introduction of a spatial component implies that the functional forms in (12), will be essentially different from those in (11). The same analysis may be extended to linking the environment to p 2 , and in investigating whether the links scale (linearly/non-linearly) during fisheries, and moratorium years. The established connection between maturation parameters (p 1 and p 2 ) and stock status has relevance for how this stock is managed. A practical use of this information is in updating the stock assessment model with a metric that quantifies stock statushigh, moderate, collapse. The status metric should determine the values assigned to the maturation parameters. Such an approach will reduce the uncertainty in modeling maturation from the autumn (survey time) to spawning time. The impact of such uncertainty reduction on the harvest rule can be demonstrated through a retrospective analysis of historical management decisions using fixed versus dynamic (based on stock status) maturation parameters.
This manuscript has used a sequential approach, where data clusters were derived from estimated (p 1 , p 2 ) tuples, which were then used as input in multiple steps of Markov chain analysis. While this step-wise approach may work for short time series (as in this paper) the approach may be prone to bias if the uncertainties at each step of the sequential process are not insignificant. This may, for example, be the case when dealing with long time series of observations that yield a large set of (p 1 , p 2 ) tuples. A viable alternative will be to use the estimated (p 1 , p 2 ) tuples as input to a Hidden Markov Model (HMM) (Rabiner, 1989;McClintock et al., 2020) and thus, derive biomass states and their associated Markov state properties in a manner that is less prone to bias and uncertainty.

Declaration of Competing Interest
The authors report no declarations of interest.