Improving the gnomonic approach with the gnomonicM R-package to estimate natural mortality throughout different life stages

Gnomonic interval model and new features

According to Caddy (1996) and Martínez-Aguilar, Arreguín-Sánchez & Morales-Bojórquez (2005), the gnomonic method is supported by a negative exponential function in which the independent variable is Δ_i, representing the temporal duration of the ithgnomonic interval; for i = 1, 2, 3, …, n, the equation is expressed as follows: (1)${\displaystyle N_i=\left\{\begin{array}{c}MLF\ast e^{-\left(M_i\ast {\Delta }_i\right)};fori=1\hfill \\ N_{i-1}\ast e^{-\left(M_i\ast {\Delta }_i\right)};fori>1\hfill \end{array}\right.}$ where M_i is the average value of the natural mortality rate, which integrates the declining death rate through Δ_i and N_i are the survivors from the previous Δ_i. The initial population for the first gnomonic interval could be assumed as: (i) the number of hatching eggs (Caddy, 1996), (ii) the mean lifetime fecundity (MLF) (Caddy, 1996), or (iii) the number of offspring per mating event (Lambert, 2008).

In the gnomonic model, the estimation of M_i for each Δ_i requires biological information about (i) the number of developmental stages throughout the life cycle i ∈ {1, …, n}, (ii) the duration of the first life stage corresponding to the first gnomonic interval ( Δ₁ = egg stage duration), (iii) the MLF, and (iv) the lifespan of the species. In biological terms, the MLF is inversely correlated with lifespan, except for specific groups such as semelparous populations (e.g., squids) (Caddy, 1996). As additional information, the duration of the other developmental stages (larvae, juvenile, adults) can be provided. The duration of the second gnomonic interval is defined as Δ₂ = α∗t₂₋₁, where α is a proportionality constant, while the successive gnomonic intervals are calculated as ${\Delta }_i=\left(\alpha \ast t_{i-1}\right)+t_{i-1}$, where i ≥ 3 up to the ith gnomonic interval. M_i is proportional to the life stage duration since G is constant, which was expressed as follows: ${\displaystyle M_i=\frac{G}{{\theta }_i-{\theta }_{i-1}}}$ where G is the proportion of the overall natural death rate, being constant for all gnomonic intervals (Caddy, 1996), which is the product of M_i∗Δ_i and ${\theta }_i=\left({\Delta }_i∕t_n\right)∕365$, representing the annual proportional duration of each interval, where t_n is the longevity of the species in days (Martínez-Aguilar, Arreguín-Sánchez & Morales-Bojórquez, 2005).

The gnomonic approach is not applicable to marine mammals or sharks; the primary assumption of the gnomonic method is based on the estimation of natural mortality from a negative exponential function, similar to those used to estimate Z (total mortality). According to this mathematical solution, there must be abundant individuals at time 0 (N₀) to enable the decay in the number of individuals in the population to be estimated. Conversely, whether only one pup or individual is born, the negative exponential function cannot be solved.

Mathematical simplification, uncertainty, and sensitivity

In this study, some equations have been modified from Caddy (1996) related to (i) the calculation of the duration of each subsequent gnomonic interval after the egg stage and (ii) the estimation of the constant proportion of the overall natural death rate (G) to improve the model performance during the computational process.

The duration of the first gnomonic interval $\left({\Delta }_1\right)$ is equal to the time elapsed after the moment of hatching (t₁). The durations of the subsequent gnomonic intervals (i ≥ 2) are estimated as follows: ${\displaystyle {\Delta }_i={\Delta }_1\ast \alpha \ast {\left(\alpha +1\right)}^{\left(i-2\right)}}$ where

Δ_i is the duration of the gnomonic interval when i ≥ 2;

Δ₁ is the duration of first gnomonic interval t₁;

α is the proportionality constant; and

i is the i th gnomonic interval.

According to Caddy (1996) and Martínez-Aguilar, Arreguín-Sánchez & Morales-Bojórquez (2005), the α and G parameters are estimated by numerical iterations using Newton’s algorithm (Neter et al., 1996). The α value is estimated by minimizing the residuals until the sum of the durations of the gnomonic time intervals is equal to t_n. Unlike the iterative procedure previously used to estimate G, we propose an analytical solution as follows:

N_i+1 = N_i∗e^−G.

This equation could be expressed as follows:

$N_n=N_0\ast {\left(e^{-G}\right)}^n$.

According to this equation, G can be estimated based on the total number of gnomonic intervals $i\in \left\{1,\dots ,n\right\}$ and the value N₀ = MLF. Furthermore, G is chosen so that the number of individuals surviving to the last gnomonic time interval is N_n = 2, following the assumption of stable population replacement with a 1:1 sex ratio; therefore, one female fulfils the requirement for population replacement if the eggs are fertilized (Caddy, 1996). The new equation for G is expressed as follows: ${\displaystyle G=-ln\left[{\left(\frac{2}{MLF}\right)}^{\frac{1}{n}}\right].}$ Therefore, only the constant proportionality (α) parameter needs to be estimated. For this purpose, the NEWUOA numerical optimization algorithm (Powell, 2006; Powell, 2008) was used, available from the “minqa” package version 1.2.4 of R software (Bates et al., 2014; R Core Team, 2020).

According to Martínez-Aguilar, Arreguín-Sánchez & Morales-Bojórquez (2005), the variability in M was assessed assuming that MLF was the main source of uncertainty, therefore simulating a total of j samples of MLF with a uniform distribution to determine the uncertainty of M_i. Then, M_i estimates were obtained from n simulations per gnomonic interval, obtaining the mean natural mortality rate $\left({\overline{M}}_i\right)$ and the standard deviation $\left({\sigma }_{\overline{M}i}\right)$. Another modification in the gnomonicM package is the assessment of the uncertainty via a Monte Carlo simulation, which was improved via the inclusion of three assumed probabilistic density functions for MLF defined as (i) uniform MLF ∼ U(MLF_min, MLF_max), (ii) normal $MLF\sim N\left({\overline{\mu }}_{MLF},{\sigma }_{MLF}\right)$, and (iii) triangle MLF ∼ Triangular(MLF_min, MLF_max, c_MLF), where MLF_min and MLF_max represent the minimum and maximum of the observed MLF, respectively; ${\overline{\mu }}_{MLF}$ and σ_MLF are the mean and standard deviation of the observed MLF, respectively; and c_MLF is the mode of the MLF in the triangular distribution. A Monte Carlo simulation must provide a correct stochastic orientation for estimating confidence intervals because the procedure (a) quantifies uncertainty based on statistical distributions derived from data rather than arbitrarily chosen distributions, (b) is unbiased, (c) is accurate, and (d) uses few distributional assumptions (Haddon, 2011; Magnusson, Punt & Hilborn, 2013). Additionally, in this study, sensitivity can be tested assuming different values in the number of gnomonic intervals, longevity, or egg stage duration in the gnomonicM package. The sensitivity analysis does not require assumptions regarding statistical distributions, such that the user can choose, even arbitrarily, the values of the input parameters (Blackhart, Stanton & Shimada, 2006; Magnusson, Punt & Hilborn, 2013).

The gnomonicM package

The source code of the gnomonicM package is freely available from CRAN (https://cran.r-project.org/web/packages/gnomonicM/index.html) or on GitHub at https://github.com/ejosymart/gnomonicM. This package has been built to provide a user-friendly method for estimating the natural mortality (M_i) and temporal duration of each gnomonic interval, the α parameter of the gnomonic model, and the proportion of the overall natural death rate (G). The main arguments from the gnomonicM package are described in Table 1, and a detailed description is presented in the package manual (https://cran.r-project.org/web/packages/gnomonicM/gnomonicM.pdf). Additional print and plot methods for the deterministic and stochastic approaches are provided to show the results.

Testing the gnomonicM package

For the application of the improved gnomonic approach and to show the functionality of the gnomonicM package, gnomonicM was tested via the deterministic method and by comparing the estimates with the results of two species reported by Caddy (1996); the species had (i) seven gnomonic intervals, (ii) longevity of one year (365 days), (iii) egg stage durations of 2 days, and (iv) MLF values of 200,000 and 135 eggs. Additionally, the methodology was applied to published data that used the gnomonic model, and the estimates were compared with the results provided by the cited authors (see Ramírez-Rodríguez & Arreguín-Sánchez, 2003; Martínez-Aguilar, Arreguín-Sánchez & Morales-Bojórquez, 2005; Giménez-Hurtado, Arreguín-Sánchez & Lluch-Cota, 2009; Martínez-Aguilar et al., 2010; Aranceta-Garza et al., 2016; Romero-Gallardo et al., 2018 for details). This approach allowed the assessment of the application of gnomonicM for different taxa (fish and invertebrates) and life histories (demersal, pelagic, benthic, and short and long-life spans). The gnomonicM package was supported by the reproducibility and verification of the results obtained from different reports, thus guaranteeing its functionality, applicability, and performance in estimating M for different ontogenetic developmental stages.

The case of Pacific chub mackerel

For illustrative purposes with the entire estimation process, Pacific chub mackerel was used as an example. Thus, this section indicates the data requirements and the steps taken to apply the gnomonicM package, highlighting its flexibility in parameter estimations, variability in fecundity, the selected probabilistic density function, and the inclusion of auxiliary information on known gnomonic intervals different from the egg stage.

(a) Choosing the number of gnomonic intervals. These intervals are defined a priori and should be well represented since they are associated with the ontogenetic development stages, exhibiting realistic subunits of biological time. For Pacific chub mackerel, the biological-ecological criteria indicated the presence of eight gnomonic intervals, which are shown in Table 2. This means that the intervals should be defined following the gnomonic time framework since they are a key input in the model.

(b) Defining the egg stage duration. For Pacific chub mackerel, Hunter & Kimbrell (1980) reported that eggs of this species hatched in 56 h at 19 °C; the author also provided information about the fluctuation of the elapsed hatching time as a function of temperature under controlled conditions, ranging from 33 to 117 h at 23 and 14 °C, respectively.

(c) Including additional information if available (such as the duration of a specific cycle life stage different from the egg stage). Occasionally, some biological information about the duration of a specific life cycle stage is reported, mainly including knowledge of early stages (e.g., the larvae stage) obtained from rearing conditions. If additional information on a known gnomonic interval (life development stage) is provided, then this time duration will not be estimated. For Pacific chub mackerel, no additional information was used in the estimation process.

(d) Including the longevity of the species. The concept of longevity in this study is related to the maximum age that a species could reach, and it is based on growth studies supported by reading otoliths, statoliths and other hard structures, such that the age structure of a population can be known (Beverton, 1987). For Pacific chub mackerel, according to Mendo (1984) and Caramantin-Soriano, Vega-Pérez & Ñiquen (2008), a longevity of 8 years (2,920 days) was used.

(e) Assigning fecundity values. Fecundity is defined as the number of offspring per mating event (Lambert, 2008), and for the gnomonicM package, fecundity represents the population at time 0; this input can be highly variable depending on several biological, physiological, and environmental conditions (Kjesbu et al., 1998; Zwolinski, Stratoudakis & Sares, 2001; Lambert, 2008). Therefore, the gnomonicM package includes two options: the first is a deterministic approach in which the uncertainty in fecundity is ignored; the second approach involves stochasticity and includes the uncertainty of the fecundity value using a Monte Carlo simulation. For the latter option, the user must assume a probabilistic density function (uniform, normal, or triangular) linked to the fecundity. The procedure provides as outputs the precision of the natural mortality value for each gnomonic interval selected in step a). For Pacific chub mackerel, a uniform probabilistic density function was used based on the fecundity values reported by Peña, Alheit & Nakama (1986) and Buitrón & Perea (1998) (Table 3).

(f) Assigning an initial value to the α parameter. This step allows a statistical solution when the gnomonic model is optimized. The α parameter has a default initial value equal to 2. In cases where the information on this parameter is limited, the user should provide an acceptable α parameter according to the taxonomic group studied, using references previously reported for fishes (demersal and pelagic fishes), crustaceans (shrimp), molluscs (cephalopod and clams), or holothurians (sea cucumbers) (see Table S1). In this way, the α parameter represents an initial value that is able to prevent the algorithm from becoming stuck in local minima.

The gnomonicM package and other methods

Compared to the method employed by Caddy (1996), the gnomonicM package uses simplified algebra, increasing its parsimony. The original gnomonic model proposed estimations of the α and G parameters in two independent numerical optimizations, while the gnomonicM package used an algebraic estimator for G and an analytic estimator for α supported by the NEWUOA numerical optimization algorithm (Powell, 2006; Powell, 2008). Another feature of the gnomonicM package is the biological and ecological sense linking the gnomonic concept to the life cycle of any given species, expressed as the development of the ontogenetic stages (Martínez-Aguilar, Arreguín-Sánchez & Morales-Bojórquez, 2005). Moreover, the gnomonicM package allows the incorporation of auxiliary information for specific gnomonic intervals related to their observed time durations; this contribution is useful for calibrating estimates of M vectors. Furthermore, the gnomonicM package provides confidence intervals (CI 95%) of natural mortality for each estimated gnomonic interval, assuming that the main source of uncertainty is the mean lifetime fecundity. The latter can be assumed to be non-informative (i.e., with a uniform distribution) or informative (i.e., with normal and triangular distributions). The choice of the probabilistic density function will affect the natural mortality estimates; however, the uncertainty regarding fecundity must not be limited to only non-informative distribution (Martínez-Aguilar, Arreguín-Sánchez & Morales-Bojórquez, 2005).

According to Quinn & Deriso (1999) and Kenchington (2014), the natural mortality estimation methods require different input data and variables, and some of these data could be difficult to obtain because they are usually associated with extensive time series (i.e., length-frequency analysis and some stock assessment models); others, such as tag-recapture and telemetry in oceans, are usually expensive and sometimes provide imprecise estimates (Hearn, Pollock & Brooks, 1998; Pine et al., 2003; Pollock, Jiang & Hightower, 2004), while some are supported by knowledge of the age structure of a population and its growth parameters (in methods based on the life history and maximum observed age of a species). These meta-analysis estimators are biological generalizations that use multiple regression analysis with independent variables, such as growth and environmental variables, whose application is limited to certain marine taxa. Regarding the uncertainty, a common feature of the above methods is the absence of an uncertainty calculation, expressed as a confidence interval (Kenchington, 2014). Although the gnomonicM package requires specific biological data, it provides estimates of M for the entire life cycles of marine organisms; estimations are commonly scarce for egg and larval stages, and the gnomonicM package also includes an analysis of uncertainty represented by the confidence intervals of M. Furthermore, the sensitivity of the estimated M values can be tested by varying the duration of the egg stage, since this stage is a critical input data for the gnomonic model, affecting the estimated M values mainly at the early stages of development. Finally, the gnomonicM package provides estimations of M via numerical optimization in comparison to deterministic estimators (see Kenchington, 2014).

The Pacific chub mackerel case and details about gnomonicM applications

The aim of presenting a new case study was to provide a comprehensive guide for data collection and obtaining results and to explain the details of the application of the gnomonicM package to avoid its misuse. Based on the gnomonic approach, biological information regarding the life cycle of a given species must include the temporal duration of the egg stage, longevity, and fecundity. For Pacific chub mackerel, the biological data were collected from several sources (Table 1). At this step, the main challenge was the different criteria among authors focusing on the definitions of biological stages and their durations. For example, in the larval stage, it is common to find subclassifications, referred to as pre-larvae, early larvae, larvae, yolk-sac larvae, and post-larvae stages; for some species, there are even subcategories defined as I or II (e.g., Martínez-Aguilar, Arreguín-Sánchez & Morales-Bojórquez, 2005; González-Peláez et al., 2015). Moreover, contrasting results regarding the development stage can be documented from rearing conditions or field observations. The solution to this problem can be addressed by the well-defined concept of the gnomonic interval (Caddy, 1991; Caddy, 1996). This means that some stages could be considered individually or grouped to make up a particular gnomonic interval considering the biological and ecological characteristics of the interval, with the purpose that each successive gnomonic interval is greater than the previous interval in its duration. If this condition is not satisfied, then the results will be biased due to the misspecification of data. Another cause can be related to the lack of specific biological data for the species studied; therefore, some generalizations could be assumed, such as, for example, using a characteristic (e.g., duration) of the development stages of other members of the same taxonomic genus, even from different geographical regions.

Once the data have been collected, they are introduced directly as input arguments (see Table 1). These inputs are simply composed of numbers, vectors, or characters depending on the approach used, either deterministic or stochastic, providing a quick, flexible, and user-friendly tool with characteristics that a software package should have (Wilson et al., 2017). This approach allows users to focus on obtaining and interpreting results rather than the calculation process. Therefore, when the deterministic approach is selected via the gnomonic function, the input data must be organized as follows:

x <- gnomonic(nInterval = 8,

eggDuration = 2.33,

longevity = 2920,

fecundity = 78174,

a_init = 2)

In the case of the stochastic approach via the gnomonicStochastic function, the input data must be organized as follows:

x <- gnomonicStochastic(nInterval = 8,

eggDuration = 2.33,

longevity = 2920,

min_fecundity = 11804,

max_fecundity = 144543,

distr = ”uniform”,

niter = 1000,

a_init = 2).

For any option selected previously, the numerical outputs can be obtained from print(x), and the graphic representation is available for the user from plot(x).

The scenarios described above represent a sensitivity analysis available from the gnomonicM package; the rationale is supported by the role of the environmental variability influencing the fecundity and the egg stage duration, both of which have impacts on the natural mortality of S. japonicus. The sensitivity analysis must be based on biological and ecological knowledge of the analysed species, such that the scenarios estimated from this approach can be plausible. In this way, the use of the gnomonicM package requires clear criteria for estimating the number of scenarios useful for each case. In this study, we tested scenarios based on the duration of the egg stage and MLF. We chose the scenario of 56 h (2.33 days) as the “best” scenario because the estimations of the durations of the gnomonic intervals obtained under this scenario were very similar to the observed values (Tables 2 and 5). Additionally, these durations increased over the life cycle following the gnomonic time concept (Caddy, 1996). The M values for different stages were similar to the reported values for the egg and larvae stages (Watanabe, 1970; Ware & Lambert, 1985), and the differences were presented for the adult phases due to the method used, such as catch curve analysis, stock assessment and empirical equations (Tables 5 and 6).

The natural mortality values estimated from the gnomonicM package for Pacific chub mackerel showed an adequate biological trajectory through the different ontogenetic developmental stages, with the highest values observed during the early development stages, indicating that the gnomonic times selected for this species were adequately grouped. Regarding ecological theory, the early stages of Pacific chub mackerel are exposed to high rates of predation by planktonic organisms and fish. This species has a variety of predators depending on the phase of development. The main predators at different sizes of Pacific chub mackerel, including juveniles, are hake Merluccius gayi and the eastern bonito Sarda chilensis (Ojeda & Jaksic, 1979; Fuentes, Antonietti & Muck, 1989). Cannibalism is another source that increases natural mortality, and several studies on this species have shown evidence of cannibalism in eggs and individuals 8-mm and larger, particularly in spawning grounds. It is important to mention that if the egg stage duration is long or if the transition from egg to larvae is slow, the organism may be subject to predation mortality for longer periods, resulting in high natural mortality (Pitcher & Hart, 1982).

Finally, the improved approach and the use of the gnomonicM package adequately represented the biological developmental stages of the species assessed in this study over its life cycle. Additionally, estimates of M produced by the gnomonicM package may provide basic input data to ecological models, enabling them to use estimates of variable M across different ages groups, mainly for the life stages affected by fishing.