Estimation of the present status of the species based on the theoretical bounds of environmental noise intensity: An illustration through a big abundance data and simulation

Sibly et al. (2005) described that most species have a fundamental characteristic to follow the theta-logistic growth trait with the convex downward trend. The fundamental yardstick of this research work builds under the deterministic setup, whereas the involvement of the external noise in any growth system is inevitable. But, the involvement of external affairs in any species growth can not be well judged only through its density dependence; it requires a further assessment. So, we frame the theta-logistic model with the stochastic analog in two directions, i.e., the discrete and continuous setup. The analysis of the discrete stochastic model is manifested by the bifurcation analysis, which shows that the attainment of the chaotic regime enhances with the increase in noise intensity level. Although the role of chaos in species extinction is debatable, a literature survey suggests that chaos with stochasticity accelerates the extinction of species. Similarly, in the case of the continuous version, we performed a theoretical study on the stochastic theta-logistic model to provide a critical value of the noise intensity parameter. This threshold magnitude acts as the sustainability criterion of any species environmental tolerability. In this connection, we use the data of four major taxonomic groups, i.e., Bird, Insect, Mammal, and Fish, from the GPDD database and classifies the species based on environmental sensitivity. The highly sensitive species have a low tolerance level, associated with the small magnitude of environmental noise intensity parameter. Moreover, the simulation prediction model on these four taxonomic classes also shows that the overall extinction probability of the considered birds in our research work is almost negligible for the current time window.


Introduction
A fundamental question in population ecology concerns the relationship between the population size and its per capita growth rate (pgr). This relationship has vital implications for understanding the animal lives, sustainability, and risk status of the species. A seminal contribution by Sibly et al. (2005) established this relationship using density-regulated population growth model based on the data of four taxonomic groups such as Birds, Insects, Mammals, and Bony fishes available in the Global Population Dynamics Database (GPDD) (NERC Centre for Population Biology 2010). For some time series of species abundance in this database, the density-regulated parameters have been allowed to take negative values, which ultimately lead to a negative estimate of instantaneous growth rate, an unrealistic phenomenon of the population ecology. But, Bhowmick et al. (2015) addressed this issue and provided a settlement by proposing a more general model incorporating cooperation and competition.
Although Sibly et al. (2005)'s proposal of using a thetalogistic model for understanding the pgr-density relationship suffers from some mathematical and model-specific limitations, it provides some fundamental insights into the population dynamics. The study of Sibly et al. (2005) and a series 1 3 of sequel studies (Reynolds and Freckleton 2005;Bhowmick et al. 2015) reached the same conclusion that most of the pgr-density relationship is concave upward (convex downward) in nature except for large mammals. This phenomenon can be well described by the theta-logistic model having the density regulation parameter less than unity. This phenomenon is probably because many species would like to spend most of their lifespan in the neighborhood of the carrying capacity because of its better survival prospect and optimizing habitat rationing. We know that carrying capacity is the only stable equilibrium point, and so naturally, species would like to spend most of their life span around the carrying capacity. There is some evidence of non-monotonic Allee type relationships between density and pgr, which is negligible in numbers compared to a concave relationship. Note that Sibly et al. (2005) and Bhowmick et al. (2015) found that in the percentage of species following the Allee, cooperative growth mechanism is very much negligible in respect to the theta-logistic growth trait. The question persists whether the species are safe from extinction while adopting the theta-logistic pattern for their density-pgr relationship under the realistic stochastic environment.
The fundamental yardstick of Sibly et al. (2005)'s work builds under the deterministic setup, which is well described through the first-order differential equation. It is worth mentioning that the theta-logistic model's discrete analog is well represented through the extended Ricker setup. A complete picture of the stable state can be obtained through the bifurcation diagram by varying the magnitudes of density-regulated ( ) and carrying capacity (K) parameters. Moreover, the model provides clear concepts of the estimated population sizes for which the species may enter sequentially into the limit cycle and chaotic dynamics.
Although it is a matter of debate whether the chaos is a precursor of extinction or not, Bhowmick et al. (2015) are informative in this connection. In a conclusive statement, Bhowmick et al. (2015) suggest that the probability of local extinction increases when the population enters the regime of chaos. Ritchie (1992) demonstrates that if chaotic dynamics occur more often, extinction of the local population is almost inevitable. Thomas et al. (1980) pointed out that deterministic chaos sometimes enhances the probability of extinction. For example, nonlinear dynamics increase the likelihood of species local and global extinction in experimental populations of brine shrimp. But chaos is not always detrimental on the ecosystem (Ruxton 1994); Huisman and Weissing (1999) describe the maintenance of chaos to be desirable for reducing the probability of global extinction to maintain species biodiversity and Bhowmick et al. (2015) suggest that system managers can regulate certain unwanted fluctuations. They will have to be careful whether the population size declines the viable size or not. In nonlinear dynamics, chaos is the measure of this system's entropy (Waliszewski and Konarski 2005). Hence, the beginning of chaotic dynamics may not always lead to extinction. Still, in the transition of chaotic dynamics, population growth is very irregular and unpredictable, and even a tiny external perturbation may lead the population to extinction, leaving a longterm impact on the population dynamics (Schreiber 2003). In summary, chaos with environmental stochasticity increases extinction probability (Berryman and Millstein 1989). Still, the coupled population can avoid global extinction utilizing chaotic dynamics (Allen et al. 1993). Lande (1993) point out that the incorporation of stochasticity will profoundly impact species growth as most of the population dynamics are driven by external fluctuations. Henceforth, IGR is the most vital parameter, which has been affected by demographic stochasticity (Loeschcke and Seitz 1991). There are also many instances where IGR is affected by both environmental and demographic stochasticity (Méndez et al. 2010). Mathematically, suppose we incorporate stochasticity in IGR of the theta-logistic model by introducing the Wiener process. In that case, the growth equations take the form of a stochastic differential equation complete with the multiplicative noise structure.
A key question for the population biologists is whether there is any theoretical bound for quantifying the intensity of stochastic fluctuation based on the theta-logistic model's parameters. We selected the theta-logistic model due to its applicability in explaining the relationship of pgr-density for most species. The empirical estimate of this bound, estimated using real data, will help provide ecologists with information on the actual status of the species in relation to its sustainability and extinction potential. A species can be threatened due to the multiplicative noise variation even if they exhibit stable dynamics under deterministic theta logistic setup. In this paper, we have threefold objectives: (i) estimate the theoretical bounds of the theta-logistic model by incorporating the multiplicative noise through IGR, (ii) comment on the species sustainability and risk status by evaluating the empirical estimate of the above bounds based on the GPDD data sets, (iii) the abundance data for species of several taxonomic groups is available in different past time windows. Based on this historical data, we like to predict the species present status through the probability of extinction.
1 3 the deterministic trait can reach a steady state. The species growth regulation can be viewed in two ways: continuous and discrete. Continuous setup is appropriate for the species with overlapping generation. Most interestingly, some species have life cycles maintained through non-overlapping generations (Rana et al. 2014). In this case, the discrete analog of the theta-logistic model is functional. Another critical question will arise how an experimental biologist can understand the extinction status of some species under a discrete growth rate. The advantage of a discrete system is that a single species model can exhibit the different dynamical phases such as stable, limit cycle, and chaos under varying parameters.
Does chaos act as an onset of species extinction? It is debatable whether chaos can be viewed as a precursor of extinction. Chaos is a non-repeating, erratic fluctuation trajectory that is bounded and sensitive to the initial condition. Under chaotic conditions, the system managers will have to be careful whether the population size declines the viable population size or not. Population size below this critical threshold may enhance the risk of extinction if stochasticity is present in the system. If we consider any simple environment, small vegetation-eating animals with high reproductive rates eat almost all the vegetation in one year. The following year, the vegetation will not have recovered, but the species richness will still be very high. Thus, the high growth rate causes a disconnect between the actual population size and the adverse effects of those individuals comprising the population. The negative effects of the actions of individuals (e.g., resource consumption) are felt by the offspring of those individuals rather than the individuals themselves. Hence, the birth of chaos occurs in the ecosystem. However, chaos is not always detrimental. Systematic chaotic intervention is desirable for reducing the probability of global extinction (Huisman and Weissing 1999).
Lyapunov exponent as an indicator of chaos: Several authors (Berryman and Millstein 1989;Thunberg 2001;Dennis et al. 2003) discuss how the Lyapunov exponent (henceforth, LE) can illustrate the chaotic dynamics in any species. The positivity of the LE signifies the chaotic regime into the dynamical system (see the diagram ref Deterministic discrete theta-logistic bifurcation diagram). Berryman and Millstein (1989) state that this positivity depends on two factors. The factors include (1) production of the offspring should be regular, and (2) positive feedback effects must dominate that dynamical system. Moreover, the species richness is displaced far from their equilibrium points due to the negative overcompensating effect. In that situation, positive feedback behavior must dominate.
Chaos in discrete theta-logistic model: an illustration: One of the primary objectives of the present study is to predict the species extinction status and estimate the amount of environmental perturbation a species can tolerate. Most of the species under consideration follow the law of overlapping generations. So, we focused our study on the continuous system for identifying the species status and the theoretical bound of environmental fluctuations. Before exploring the continuous version of the theta-logistic model, we feel that it would be appropriate to monitor the discrete system through bifurcation analysis. We need to identify when the chaotic regime starts by varying the values of the reproduction parameter. A bifurcation diagram can determine parameter values for which the dynamical system moves to period-doubling, limit cycle, and extinction corridor. This analysis will be helpful to predict the local extinction status of the species. We focus our study on a non-overlapping species, Oncorhyncus gorbusha. We compare the estimated values of the reproduction parameter with its theoretical values for which the system will move to a chaotic regime. In order to understand the nature of the bifurcation diagram, one should follow the protocol of Thunberg (2001). In the spirit of Thunberg (2001), we propose a theorem to justify whether the diagram for the species possesses the one-humped character or not. The area of the hump is an indicator of the onset of the chaotic region: more the area, less the risk of entering the species into the chaotic zone.

Species extinction status under continuous stochastic thetalogistic model:
We will express the theta-logistic growth model as a stochastic differential equation by incorporating the multiplicative noise in the reproduction term. We can theoretically estimate the upper bound of the environmental fluctuation using Ito calculus (Golec and Sathananthan 2003) and the Lyapunov function Kot (2001). The risk status of the species can be identified by comparing the theoretical estimates with that of the fitted model estimate of the noise intensity parameter. This comparison is made based on a systematic and thorough study on 1780 time series of species from different taxonomic groups. A natural and fundamental question will automatically arise: What is the current status of the species concerning the calendar year 2021? We have planned to propose a simulation-based analysis for identifying the extinction status of the species, which is contemporary to the PVA analysis in IUCN (Stephens 2016).

Theta-logistic model and density regulation
The growth dynamics of any species is characterized by the convolution of two opposite growth pulses, i.e., (1) natural inherent affinity of the species for unlimited growth and (2) the nature's interference as the negative feedback. The overcompensatory dynamics of the density-dependent phenomena is one of the significant components behind these two characteristics (Nicholson 1954;Kot 2001;Bhowmick et al. 2015). The negative feedback rule of any ecosystem hinders the species unlimited growth as a consequence based on the available resources and favorable conditions species have a chance to change their habitat through the adaptation on the stage of evolutionary theater due to the long time survival prospect (Johnson and Agrawal 2003). Although several controversial statements exist (Andrewartha and Birch 1984;Foley 1994) regarding the influence of the density dependence, recently many authors (Fagan et al. 2001;Henle et al. 2004;Heering and Reed 2005;Lande et al. 2009) acknowledge the importance of it upon the species growth regulation.
The classical model describing the density-dependent phenomena is the logistic growth law (Kundu et al. 2021). The linear decrement of the relative growth rate with increasing population abundance will put a restriction on the explanation of density dependence through the logistic growth law (Ross 2009). The study of Sibly et al. (2005) reveals that the density dependence is not well captured through the linear association between the species fitness and its abundance. The author demonstrates this relevance by showing that the fitness profile of most of the species follows a non-linear trend, i.e., either convex or concave. Several authors (Lande et al. 2009;Ross 2009) also admit the non-linear associations between the species fitness and its abundance to illustrate the density dependence issue. In this connection, Sibly et al. (2005) uses the theta-logistic growth curve model to explain the non-linear trend in the species growth profile. The mathematical expression of this growth equation is given by Here, N(t) represents the population density at any time t. The model parameters K, represent the carrying capacity and the curvature parameter of the theta-logistic model. Sibly et al. (2005) analyze the size-rgr relationship for 1780 time series from the Global Population Dynamics Database (GPDD) with this theta-logistic growth trait. The time series consists of four taxonomic groups, i.e., Aves, Insect, Fish, and Mammals. The author found that the density regulatory behavior in most species growth profiles is maintained by < 1 . But, for the large mammals, this density regulation is followed by the relation > 1 . Hence, it is clear that the theta-logistic growth trait is one of the stable dynamics in the species growth regulation. Apart from this growth regulation, some species follow the Allee growth law in their life cycle. Note that, Sibly et al. (2005) estimated this minor percentage as 0.2, and Bhowmick et al. (2015) evaluated this figure to be 3.92% for cooperative species, where cooperation is a particular case of an Allee effect. (1)

Discrete analog of theta-logistic model
We have already mentioned that the theta-logistic growth dynamics well explain the size-rgr relationship in most of the species. Sibly et al. (2005) use the deterministic continuous setup to elucidate this scenario. It is worthy of mentioning that this continuous deterministic form of the theta-logistic law does not relate to species extinction. But, the one-dimensional discrete Ricker map of the thetalogistic growth equation possesses several characteristics, viz., 2-period cycle, 3-period cycle, period doubling, limit cycle, chaos, etc., for a specific range of parameter values. Although the association of the chaotic dynamics with species extinction is still debatable, we believe that these unpredictable, chaotic dynamics certainly increase the probability of extinction very much. So, there is an urge to project a spotlight on the discrete dynamics of the theta-logistic growth law. The discretized version of theta-logistic law is given by where the notations have their usual meaning.

Discrete system and chaotic regime
Chaotic oscillations are generally observed after attaining the asymptotic size. Hence, the maximum population size or the carrying capacity plays the role of attractor in any ecosystem. Thus, from the mathematical perspective, one can say that carrying capacity creates a basin of attraction after which the orbital period of the growth trajectory becomes infinite. This would lead to chaos in the dynamic system. In this connection, we provide the following theorem to elucidate the analytical property of the bifurcation diagram. Proof Let us consider f (N(t), ) in Eq. 2 by In order to find the maximum point, we need to solve the Note that, Thunberg (2001) proposed that any onedimensional map from an interval I to itself will be unimodal (one-humped) if ∃ any point c ∈ I = [a, b] such that f (N(t), Θ) has a unique maximum value at c, i.e., the function is monotonically increasing on [a, c) and then decreasing on the interval (c, b]. Here the one-dimensional map of the theta-logistic model 2 is confined in the closed interval [0, K]. Moreover, the specified result mentioned above depicts that N(t) ∈ [0, K] . This ensures that the bifurcation diagram should possess the one-humped characteristics. ▪ The numerical output of the bifurcation diagram (see the right panel of the Fig. 1) also depicts that it possesses the one-humped characteristics. The stability analysis in difference equation (Kot 2001) shows that the discrete thetalogistic system 2 becomes stable if it maintains the condition 0 < r K < 1 . In addition, the growth Eq. 2 undergoes damping oscillation for 1 < r K < 2 , where the orbital period of the one-humped mapping 2 revolves within the power of 2. The beginning of the limit cycle is observed in the system 2 when r K > 2 , i.e., the one-dimensional mapping is now beginning to follow an odd-order period. The phase of the limit cycle then gradually converts to a chaotic stage when the period of the orbit becomes three (Li and Yorke 2004).
Since the bifurcation diagram is produced by varying the parameter values, it would become necessary to identify the threshold value of the parameter responsible for attaining the chaos. The Lyapunov exponent (LE) should be the formidable measure in this case. It is defined as where N 0 be the fixed point of the mapping N(t + 1) = f (N(t), Θ) (Kot 2001). The positive magnitude of LE concerning the bifurcation parameter predicts the onset of chaos in any deterministic discrete system. The numerical analysis of the expression 4 gives that critical value of the parameter for the discrete theta-logistic map. The horizontal red line in the left panel of the Fig. 1 is used for the separation of the positive and negative zone of the LE, which also helps to recognize the threshold value of the IGR parameter to accomplish the chaos. This bifurcation analysis is being carried out on the abundance data of Oncorhyncus gorbuscha (GPDD ID-1844) available on the GPDD website. The diagram reflects that the pink salmon can achieve their stable state with a very low reproduction rate.
Moreover, the chaotic regime also occurs with a minimal level of IGR. But, the reproduction rate of the pink salmon is not too low at the Dean Channel of British Colombia, Canada, from which the data is collected (Armstrong et al. 2018). We believe the deterministic process is unable to highlight the true story beyond the species growth regulation. In this connection, we develop the following section to nurture the species growth regulation in any noisy interactive system.

Stochastic setup
The growth dynamics of each species are affected by the internal and external disturbances (Chakraborty et al. 2017).
The internal fluctuations among the species are very nuanced dynamics which are well maintained through the densitydependent mechanism. But, the involvement of external affairs in any species growth can not be well judged only through its density dependence; it requires a further assessment. However, it is observed that in a large-sized population, fluctuations in the observed number of individuals are small. But, it seems that due to the small population size, internal fluctuation plays the role of a positive catalyst in the extinction of the species (see Fig. 2). In growth curve modeling, it is very much essential to incorporate the unexpected changes into the dynamical system (Chakraborty et al. 2017). Stochastic differential equations (SDE) is the natural choice to model any dynamical system when the species confront any randomness (Bishwal 2008;Anderson 2013). In addition, environmental fluctuation is a linchpin in any ecosystem. Therefore, most of the natural phenomena do not follow strictly deterministic laws, rather they oscillate randomly about some average value so that the deterministic equilibrium is no longer an absolutely fixed state (Arditi and Ginzburg 1989). This analyzes the deterministic growth curve models under stochastic perturbations as a promising and important field in the present decade. Any SDE can be formulated by introducing randomness into a deterministic dynamical system described by an ordinary differential equation (ODE). The simplest form of an SDE is given by where G(y(t)) is a function representing the deterministic response of the system, W(t) is a Wiener process representing the stochastic driving, and is a scaling constant representing the intensity of the noise. The term dW t is the mathematical demonstration of the physical noise. Any external perturbations directly affect the species reproductive rate in such a way the parameter IGR (r) will no longer be a constant one. As a consequence, the parameter r behaves like a random variable with certain mean r and standard deviation dW (t) dt . This leads the population size of the species (N(t)) to a random variable. It is worthy of mentioning that the incorporation of such types of characteristics makes the deterministic theta-logistic model (1) the stochastic one, which the following SDE can explore.
Similarly in the discrete setup the above Eq. 6 can be written as

Remark 1
In a deterministic setup, there exist two equilibrium points for the theta-logistic model; one is the extinction state, i.e., 0, and the other is near carrying capacity. Among these, 0 is the unstable equilibrium, and another is stable. But when we give a sufficient perturbation into the system, then the stability pattern changes. The zero stable equilibrium also acts as a stable one, and there may be a false interpretation. It can also be explained from Fig. 3. The diagram exhibits that the simulated sample paths of Microtus arvalis (GPDD ID 9295) follow a to and fro movement in between the two fixed points, i.e., zero (0) and their asymptotic size. But the mean profile, i.e., the blue dashed line, preserves the deterministic sigmoidal property of the underlying growth dynamics.
(7) Fig. 1 Bifurcation diagram of the population size (N(t)) of Oncorhynchus gorbuscha (GPDD ID-1844) expressed by the deterministic growth Eq. 2 is present on the right panel of the figure. The left panel figures indicate the Lyapunov exponents. Here we use the estimated values of K, , i.e., 12068.92, 0.1, respectively with the initial population size 6900 to draw this diagram (see the sub-figures a, and b). Note that in the sub-figures c, d; e, f; g, h, we gradually increase the magnitude of by 50%, 75%, and 90% respectively to observe the dynamical change in the equilibrium points ◂ Fig. 2 Bifurcation diagram of the population size (N(t)) of Oncorhynchus gorbuscha (GPDD ID-1844) expressed by the stochastic growth Eq. 7 with respect to the IGR (r). Here we use the estimated values of K, , , i.e., 13,176.7, 0.65, and 3.29 respectively with the initial population size 100 to draw this diagram (see the sub-figure a). Note that in the sub-figures b, c, and d we gradually increase the magnitude of by 50%, 75%, and 90% respectively to observe the dynamical change in the equilibrium points It is worth mentioning that we use the bifurcation analysis for predicting the species growth status for the discrete setup. Moreover, for the continuous process, i.e., Eq. 6, we follow the standard protocol of the stochastic stability analysis. Now, for the sake of simplicity let us consider the transformation y(t) = N(t) K − 1 in the Eq. 1. Now, the modified version of 1 can be expressed as Note that the substitution of = 1 in Eq. 8 provides the same setup for the logistic model (Golec and Sathananthan 2003). So, in the spirit of Golec and Sathananthan (2003), we will perform our further analysis based on Eq. (8). Note that the expression 8 is the linearized version of the original thetalogistic growth law. It is a well-known fact that the linearization technique helps to remove the non-linearity from any dynamical system (Saha et al. 2013;Bhowmick et al. 2016). We also follow this similar approach on the parameter , which is responsible for the non-linearity in the theta-logistic model and thus obtain Eq. (8).
Let us now suppose that growth rate IGR (r) is a random variable with r(t) =r+ ∈ t , where r > 0 , ∈ t is a Gaussian white noise with a time-varying intensity 2 such that ∈ t = dW(t) . So, the stochastic version of (8) can now be expressed by the following Ito-type initial value problem We consider the initial conditions y 0 > −1 , where needed, where y(t) = y(t, y 0 ) . The principle of convergence for Thus, the stability of zero solution of (8) will guarantee the stability of the equilibrium N(t) = K . The Theorem (5.1) deals with the stochastic stability of the zero solution of Eq. (9).
Theorem 5.1 Let us consider the stochastic logistic model 9. If holds for t ≥ 0 , the equilibrium position y = 0 is stochastically stable i.e. for every ≥ 0 we have when y 0 → 0.

If in addition,
for t ≥ 0 , the equilibrium position is stochastically asymptotically stable in large on the interval (− , ∞) in addition for every y 0 ≥ −1 , we have The simulated time series of the stochastic thetalogistic model (6)  Proof Let us consider the Lyapunov function V(t, y(t)) such that it has continuous partial derivatives up to 2 nd order in y(t) and to the first order in t. Moreover, the process y(t) satisfies an Ito differential (8). Thus, So the function (13) has an Ito differential where the first part denotes the Markov generator (LV(t, y(t))) with a and b as the drift vector and diffusion matrix respectively.
Here, and Thus, Here, LV(t, y(t)) = −Ky 2 (t) r − 1 2 2 K . This analytical form is indicated as the Markov generator. Now, for the sake of stochastic stability, the Markov generator should be negative (Khasminskii 2011). So, the required condition is given by, Here the solution y(t) = 0 , i.e., a line contained in a domain t × U . Moreover the Lyapunov function V(t, y(t)) is positive definite and LV(t, y(t)) ≤ 0 for y ≠ 0 , then y(t) = 0 be stable in probability, which is nothing but that the solution be stochastically stable, which proves the first part of the theorem. Moreover, inf t>0 V(t, y(t)) → ∞ as |y(t)| → ∞ . Hence, a sufficient condition for the solution y(t) = 0 be stable in large is that ∃ a positive definite function V(t, y(t)) such that LV(t, y(t)) is negative definite and inf t>0 V(t, y(t)) → ∞ as |y(t)| → ∞ . This shows that the solution be asymptotically stochastically stable in large over the interval (− , ∞) . ▪

Remark 2
The model 9 reflects the impact of the density dependence with the stochastic analog in such a way that the ecologists can easily infer about the "variation thresholds" beyond which the species have a chance of extinction. Based on the density-dependent selection process, MacArthur (1962) conclude that the carrying capacity of any species is always maximized due to the principle of evolution. Although the conclusion made by MacArthur (1962) is acknowledged by Roughgarden (1979), Charlesworth et al. (1994), etc., recently Lande et al. (2009) shows that this proposition is not valid for the stochastic environment. The external fluctuations produce a trade-off between the selection of species with an increased intrinsic growth rate at the small abundance and the selection for larger carrying capacity at the high population density, summarizing the notion of the r − K strategy (Pianka 1970;Desharnais and Costantino 1983;Boyce 1984). The ubiquitous influence of the natural selection theory ("r" and "K" selection strategy) on the species growth trait conveys Fagan et al. (2001) to categorize the population dynamics into several groups. The association between the environmental noise intensity ( ) and the intrinsic growth rate (r) derived by Fagan et al. (1999) will deliver a distinction between these classifications. So it is very much necessary to put an environmental restriction delineating a correlation between the species IGR and the external fluctuation.
The application of the stochastic calculus upon the model 9 provides an environmental bound to the density regulation. Since the contribution of all the parameters is evident in the growth of any species, so the bound should be a function of all the growth parameters. We found from the Theorem 5.1 the analytical expression of such kind of environmental restriction. The relation (10) depicts that the environmental noise intensity has a strong positive correlation with the species IGR (r). The negative association between the density regulation and the environmental parameter ( ) elucidates that the low density regulated species have a chance to dispute with the more level of external fluctuations. Moreover, the carrying capacity is also inversely related to environmental disturbances; i.e., it is quite challenging for the species with a small amount of carrying capacity to sustain in a system with a high level of environmental fluctuation. It is relatively infrequent that the large carrying capacity species can better withstand the external noise (Lande et al. 2009).
Although the categorization of the species dynamics by Fagan et al. (2001) is regulated by these characteristics, the relation between the external noise and the species IGR described by the author lacks the necessary information about the carrying capacity and the density regulation.

Data
In this article, we use the population abundance data of different species from the GPDD (NERC Centre for Population Biology 2010). Fegraus et al. (2005)  Consequently, like most of the authors (Fagan et al. 2001;Sibly et al. 2005;Bhowmick et al. 2015), we incorporate the time series data of these first four taxonomic classes in our analysis.
Remark 3 Very few of the time series in the GPDD represent exact census values (those that do come largely from captive or island introduced populations that may not be representative of species trends as a whole). It is worthy of mentioning that the status of a species in a captive environment (extinction/threatened, etc.) may not always be the same as the local states of the species in a wild environment. But, the behavior of a species in a wild environment can be predicted well using the data of captive individuals based on the machine learning mechanism such as random forest, support vector machine, and artificial neural network, where Rast et al. (2020) used acceleration data as a tool for prediction. It is also true that the captive population data may not always be a good representative of the species' overall trend in the wild environment for a specific locality. But, if a species is showing threatened/extinction status under the captive environment, this information can be used as an early warning tool for inferring the localized status of the species in a wild environment. Less genetic variation is one of the prime reasons for a captive species moving towards extinction states. For example, Fraser (2008) studied that the stability of the salmonid species decreases gradually under the captive environment for the genetic diversity. But, at the same time, some of the species have their recovery mechanism. They can return from their extinction state even if there is a reduction in genetic variation. So, studying growth traits in a captive environment may be a meaningful tool for identifying the whole species trend in the wild.

Remark 4
The sampling protocol in most of the time series of the GPDD occurred by count. It is worthy of mentioning that the count data and the relative density are indeed two different measures of the population dynamics. But, they have a relationship. Relative abundance is defined as the count per unit area (Krebs 1978). So, mathematically one is a scalar multiple of the other. The estimates of model parameters may be different for the two cases, but the estimated population trajectories must be pretty similar. So, it hardly matters whether we are analyzing the proposed methods either on count or abundance data. Moreover, there are many instances in which the researchers used the raw GPDD data directly to focus on the local extinction (Saha et al. 2013;Sau et al. 2020) and behavioral issues, viz., cooperation (Bhowmick et al. 2015) of the species.

Data filtration
We have already mentioned that Sibly et al. (2005) used 1780 population time-series data from the GPDD database to study the density regulatory behavior on the species growth dynamics. The analysis shows that the fitness profile of most of the species (1585 time series, i.e., 89.04%) follows the concave upward trend with < 1 . However, 934 and 124 time-series data remain with the negative and zero estimates of the density regulatory parameter, respectively. The studies of Sibly et al. (2005) and Bhowmick et al. (2015) reveal that the negative estimate of any growth parameter leads to the immediate extinction of the population. Moreover, the zero estimates of the curvature parameter turn the growth profile of any species in an unbounded fashion, which does not meet the criteria of density regulation. So, it is better to neglect these datasets for further calculation. Note that in the rest of the 527 time series, we observe that there exist 7 cases with GPDD ID 1389, 6162, 6320, 6378, 6511, 6649, 6875, whose growth profile is much more synergistic with the cooperative attitude than the theta-logistic trait (Bhowmick et al. 2015). Consequently, in the remaining 520 time series, we again find 19 datasets representing the taxonomic classes other than the Mammals, Aves, Insects, and Bony fishes. As we mentioned above, these four taxonomic groups are most abundant, so it will be better to pursue one part of the analysis with these 501 time series from the database. Additionally, we will use the remaining 3367 GPDD IDs for identifying any other time series, whose RGR profile will follow the concave upward trend. The summary of this whole filtration process is present in Table 1 of the supplementary file.

Estimation process
The estimation process associated with the discrete-time analog becomes difficult due to the non-explicit nature of the conditional dynamics in the sampling process of SDEs. So, applying standard techniques, like grid-search, least square method, etc., is quite tough for pursuing the regression work with the SDE. In this connection, we need to construct a different framework that will delineate the estimation process when the underlying system becomes stochastic. However, any SDE can be generally expressed by resents the corresponding slowly varying continuous component, i.e., the drift coefficient and the rapidly fluctuating continuous random component, i.e., the diffusion coefficient of the SDE (17) respectively. Note that the unknown component Θ(⊆ ℝ) is the parameter space, which needs to be estimated. However, we find a list of authors who put their great effort into developing different estimation processes for the stochastic setup (Dacunha-Castelle and Florens-Zmirou 1986; Dohnal 1987; Gallant and Long 1997;Nicolau 2002Nicolau , 2004Alcock and Burrage 2004;Uchida and Yoshida 2005;Brouste and Iacus 2013).
Here we use R software (R Core Team 2019) with the version 3.6.1 for pursuing the whole estimation process. A couple of packages already exists in this software to perform the estimation process. We enlist the name of these packages with their functional properties in Table 1. In our case, we choose the package Sim.DiffProc from this list. Note that it does not approximate the transition density function of (17) (Guidoum and Boukhetala 2020) rather the sample path N(t) is being approximated in such a way that the discrete analog will produce a precise likelihood function.
Here the fitsde function is acting like a precursor in implementing the pseudo-maximum likelihood estimator (henceforth, PMLE) through the Sim.DiffProc package. The principal arguments require to drive the fitsde function is the univariate time series object (ts) and the initial guess (start) for commencing the process. The dataset collected from the GPDD database represents a discrete time-series object sampled over different time windows. Now, the estimation process is being conducted by the routine pmle, where one needs to specify the underlying method. The standard methods that exist to pursue this process are the Euler (Yoshida 1992; Florens-Zmirou 1989), Ozaki (Ozaki 1992), Shoji and Ozaki (Shoji and Ozaki 1998), Kessler (Kessler 1997), etc. We implement four of these methods in the estimation procedure. As a consequence, the output returns the magnitude of the estimated model parameters for all of these mentioned processes. Now we require standard model selection criteria to select the best estimate. Here we use Akaike Information Criterion (AIC) as a statistical measure for the model selection purpose, which is being reconstructed by Uchida and Yoshida (2005) for the diffusion processes. Note that the outcomes of this research preserve the property of AIC, i.e., the lower the AIC value, the better will be the model. In our case, we initially started this estimation process with 537 time-series data. However, the convergence of the PMLE takes place for only 483 timeseries data, i.e., the execution of the estimation process runs smoothly for 89.94 % of the whole datasets. We mention the estimated values of the model parameters for 483 time series in the Supplementary Material. There exists an alternative measure of parameter estimation, which is available for the stochastic Gompertz model in the research article of Knape and de Valpine (2012). (2012) introduce the measurement uncertainty in the discrete Gompertz model. The research article used the additive process error, which is manifested on 627 population time-series data of GPDD. The authors termed the process error as the uncertainty in the system. As discussed by the authors, Gompertz models can be transformed into a simple log-linear population structure. This might provide a decent first-order approximation for many populations except those undergoing strongly non-linear dynamics. In this case, the parameter estimates are biased, and magnitudes are unknown and difficult to estimate. We deal with the stochastic theta-logistic model, which possesses a highly non-linear structure. Note that the estimation process with lag-1 autocorrelation structure of errors using the Kalman filter is intractable for the stochastic theta-logistic model.

Sim.DiffProc Estimation process is conducted by constructing pseudo-maximum likelihood method by the function fitsde
Guidoum and Boukhetala (2020)

Remark 6
We have obtained unreliable estimates of the model parameters while fitting the stochastic theta-logistic model to some of the time series from the GPDD data. The parameter estimates are unreliable if it produces an unrealistic value of the ecologically meaningful model parameters.
In our study, we obtain the negative values of the parameter in many cases, which is not feasible under any ecological system. There are the four main reasons for these unreliable estimates.
• We observed both maximum and minimum RGR values for the same population size for various time series of species abundance. The estimation of model parameters based on the above unrealistic RGR data is impossible by the pseudo-likelihood procedure. We face such difficulty in the case of several population time series with GPDD ID 6471, 6497, and 6503. • Another example of some unrealistic time series has been observed for the GPDD ID 8573 and 8653, where more than 50% of the data are clustered around the neighborhood of the initial population size. We do not process in estimating the model parameters for such time series. • We face similar problems for two erratic time series with GPDD ID 9215 and 6508. • We have fitted the size-RGR regression based on the logistic model for the three GPDD ID 6678, 9488, and 9856. In this case, the RGR values are realistic, so they are reliable. But, we are unsuccessful in obtaining the reliable estimates of the model parameters by the pseudolikelihood method due to the failure of convergence criteria.

Mechanisms of species extinction
Theoretically, null population size is defined as the extinction state of the species. However, the definition of extinction has different variants in the domain of ecology. In the present discussion, we will explore this classification with clarity: (i) The offspring production of a species will be impossible to achieve a population size 1. The population size with unity ultimately compels the species to an extinction state as mentioned by several authors (Saha et al. 2013;Sau et al. 2020;Saether et al. 1998). (ii) Allee effect in population can also be treated as one of the extinction vortices of the species. The fitness of the species must be negative when the population abundance falls below a critical threshold (Allee threshold) value (Courchamp et al. 2008;Saha et al. 2013). The obvious consequence of this event is population extinction. This phenomenon is known as strong Allee effect in ecology. (iii) The seminal contribution by Rosenzweig (1971) and many subsequent studies (Huffaker et al. 1963;Luckinbill 1973;Veilleux 1979;Fussmann et al. 2000) suggest that an enrichment that is perceived as beneficial for the growth of populations might have the potential to cause destabilization of the ecosystem and eventually extinction of the predator-prey populations. (iv) Last but not least, environmental fluctuation is the most vital component for population extinction, where the population size can be reduced drastically due to sudden unexpected natural calamities (Lande 1987(Lande , 1993. We have described the above four vortices by a flowchart as depicted in Fig. 4. We define the extinction concept from the facet of the simulation process. We generate the population time series for a future time window based on a simulation experiment's history data. We proclaim the species is on the verge of extinction; if at least 50% cases from the entire time series, the population abundance falls below the initial population size of the cohort.

Proposed simulation algorithm and probability of extinction
We consider very few data among 447 time series to pursue the simulation work. The selection process is manifested by the magnitude of the percentage of variation (henceforth, PV). This measure is defined as the percentage increment of the theoretical magnitude of the noise intensity parameter ( th ) from its estimated value ( ̂ ). Here, we use the Theorem 5.1 to enumerate the magnitude of th . Note that the species will be in relatively safer zone, if ̂ does not exceed the theoretical bound, th (i.e., ̂≤ th ). Therefore, those species now became the cause of concern for which this relationship is disrupted. Hence, we continue our further analysis with all these species datasets.
Here, we observe that for each of the four taxonomic classes, the PV range is too wide between the theoretical and estimated magnitude of . So, we partitioned the range of PV into three categories, i.e., (i) low range of variation; (ii) moderate range of variation; (iii) the high range of variation. We mention this classification in Fig. 5. However, it is quite clear from the classification that the species belonging to the last two clusters will be more prone to extinction than the first category. The deviation between ̂ and th is playing the major role. In this connection, we consider an interval I = [a PV , b PV ] to select those species datasets, where a PV denotes 0.5% of the maximum PV, and b PV indicates the maximum percentage of increment. Thus, the time series should be selected so that the PV magnitudes should belong to the last two classifications and that particular interval I. Hence, the number of datasets is now being reduced to 71, i.e., 15.92%, which are on the verge of extinction.
The GPDD species abundance data were collected on different time windows. Overlapping windows with a wider span are hardly available, so comparing extinction status based on non-overlapping time windows is not meaningful. To check the overlapping and non-overlapping time windows, we have drawn a diagram depicted in Fig. 6. From this figure, we can conclude that the assessment of extinction status based on the Insect data is meaningful compared to other taxonomic groups. A significant overlapping period between 1970 and 1990 is available in the GPDD database for the Insect community. To overcome this problem for the different taxonomic groups, we will propose a simulated prediction model for determining the extinction status of the species for the current time period (2021). Note that the simulation process will now be manifested on the above 71 time-series data. The protocol of the simulation experiment is given in the following: 1. Generate a time series data for q time points following the growth Eq. 6, where q = |2021 − t f | and t f be the year up to which the data is collected. 2. Repeat the process for m times to generate a m × q (m > q) order matrix, where each row represents a single sample path. 3. Check the "extinction status," as described in Sect.
"Mechanisms of species extinction" for each of the m realizations. 4. Identify the realizations for which our rule of "extinction" is applicable and mark them by number 1 and others by 0. 5. Let p denote the number of cases marked by 1. So, the probability of extinction (henceforth, PE) for these m realizations is p m .
6. Repeat the processes 1-5 for s (> m) times to generate a distribution of PE. 7. Finally, we have identified the mean of this distribution as an estimate of PE for any particular time series data.
Remark 7 Knape and de Valpine (2012) rightly pointed out that it would be difficult to confidently identify the model parameters and the evidence of density dependence while applying measurement uncertainty in time series. But instead, this is more a problem of statistics than the structure of the GPDD data while dealing with a short and sparse time series. We used the pseudo-likelihood method of estimation instead of the general likelihood for the reliable estimates of the model parameters for the stochastic theta-logistic model. For the deterministic case, we have evaluated the bootstrap confidence intervals of the model parameters by permuting the errors of regression in several bootstrap samples. Finally, we increase the population time series by the simulation protocol to predict the extinction probability, which lengthens the existing abundance data. Hence, we believe these three tools that (a) (b) (c) (d) Fig. 6 The specified time domain of the species of four major taxonomic classes, i.e., Bird, Mammal, Fish, and Insect from the GPDD, which is being used in our analysis we have used will indeed nullify the effect of small and sparse data to predict the species future growth status.

Discussion
The dynamical relationship between the population density and its fitness can be well understood by the density regulation in any species growth (Sibly et al. 2005(Sibly et al. , 2007Ross 2009;Bhowmick et al. 2015). The classical growth model describing the density dependence is the logistic growth equation. However, the linear decrement of relative growth rate (henceforth, RGR) with the increasing population size puts a restriction on the explanation of density dependence through the logistic growth law (Ross 2009). In this connection, Sibly et al. (2005) put a significant effort into justifying that the species-inherent growth dynamics can not be well described by this linear association. The authors explain that the fitness profile of most of the species is depicted through the non-linear structure (i.e., either convex or concave) depending upon its density regulation. The empirical study of Sibly et al. (2005) delineates that the density regulation in most of the species is maintained through the theta-logistic model. The authors described that most species have a fundamental characteristic to show the non-linear fitness profile with a concave upward trend. Consequently, it seems that at the low population density, the fitness level of the species is going to 0. The species then have to produce more offspring to reach their asymptotic size. Note that the research work of Sibly et al. (2005) is demonstrated based on the continuous setup. However, we have already mentioned in "Aim of the present study" that in any ecosystem, two types of growth generations are observed, i.e., overlapping and non-overlapping. The first one can be expressed through the continuous growth equation and the latter with the discrete growth law. Since we focus on the thetalogistic growth regulation besides the continuous system, we must spotlight the discrete growth dynamics. So, we construct Sect. "Discrete analog of theta-logistic model", where we take the consortium of bifurcation analysis to describe the inherent growth mechanism of the discrete growth law. The characteristics of any discrete system can be well understood through the bifurcation analysis. In this connection, we propose the Theorem 4.1 to elucidate the intricate property of the discrete theta-logistic growth dynamics. The numerical simulation on this discrete map produces the bifurcation diagram Fig. 1, which shows the attainment of all the phases, i.e., 2-period cycle, limit cycle, chaos, etc., of the discrete theta-logistic dynamics. In this diagram, we consider Oncorhyncus gorbuscha as a testbed species to describe the bifurcation analysis. Note that these cases' consideration is driven by choice of IGR and density regulated parameter values. The study of Clark et al. (2010) reveals that the two parameters r, can not act as the independent parameters in any growth law. Instead, the compilation of these two parameters will always provide the necessary information about the species' asymptotic size as well as the stability of the concerned system (Saha et al. 2013;Kundu et al. 2018). Note that the increase in the density regulation leads any species towards the Malthusian growth. Consequently, if the species do not optimize their offspring production rate, a chaotic situation will ensue due to the limited resources. Hence, the species have to minimize the IGR so that stability is maintained throughout the ecosystem. Figure 1 elucidates this scenario. It is shown that if the IGR values are increased with the enhancement of the density regulation, the chaotic regime will be found in the discrete growth dynamics, leading any species to extinction.
Moreover, a hump has been formed for each period-doubling limit cycle case in the case of the bifurcation diagram of the discrete one-dimensional map. This period-doubling limit cycle ultimately converts to a chaotic regime in the long run while the bifurcation diagram progresses. Such period-doubling limit cycles are visible in Fig. 1. The natural question arises whether the bifurcation diagram remains unaltered if we introduce the noise in the deterministic system. For illustration, we have chosen a fish Oncorhyncus gorbuscha, which is known to be a species with non-overlapping generation due to its anadromous character (Kot 2001). Hence the discrete theta-logistic model is appropriate to describe the growth profile of the species. So the bifurcation diagram concerning the growth parameter IGR is depicted in Fig. 2 under a stochastic setup. The diagram exhibits that under the stochastic environment, the stable character of the species growth dynamics is maintained up to a comparatively high magnitude of IGR (r = 3.1) than the deterministic process. Moreover, just one hump is visible in the stochastic case compared to multiple humps of deterministic setup, which implies that species are more prone to reach the chaotic zone under the noisy interaction quickly. Moreover, when we increase the noise intensity by 50%, 75%, and 90%, the species rapidly moves towards the chaotic state (see panels b, c, and d of Fig. 2). So there is a strong association between the chaotic regime and the environmental variation. Sibly et al. (2005) identified 527 population time series out of 1780 GPDD ID for which the estimated value of the density regulation parameter theta lies in the interval (0, 1) . The author used the grid-search technique for their analysis. We initially obtained 501 time series, with ̂ ∈ (0, 1) , by analyzing the same data through the non-linear least square. Note that the non-linear least square estimates with optimum statistical properties are always preferable to the estimates obtained through the grid-search method.
But, we have identified 36 more time series with estimated theta that lie in the open interval (0, 1) when analyzing the additional 3376 GPDD ID. So combining the results of two studies, we can conclude that the percentage of concave upward size-RGR relationships persists in a 10% time series 1 3 (537 out of 5156) where ̂ lies in the interval (0, 1) . We incorporate stochasticity in the deterministic model through the parameter IGR. We estimate the parameter by the pseudolikelihood method, where the cases in the population time series with less than 1 are further reduced. In this case, we obtain 203 time series out of 5156 (4%) for which the estimated theta is contained in the interval (0, 1) . We also identify the number of species in Table 2, which shows the distinction between the concave upward and convex upward growth traits due to the mapping from the deterministic to stochastic system. Moreover, we estimated the confidence interval of the model parameters by the general frequentist approach, which are depicted in the Supplementary Materials.

Different avenues of population viability analysis
Population viability analysis (henceforth, PVA) is one of the popular methods in ecology to explore the species extinction status under some specified conditions. Several ecologists acknowledge the escalating role of PVA in species conservation. The concept of PVA initially comes from the perception of minimum viable population size (henceforth, MVP). It is assumed that the manifestation between MVP and PVA is first noticed in the work of Shaffer (1981). The author used the stochastic population dynamics simulation model with certain demographic structures to enumerate the risk of extinction and the MVP for the Grizzly Bears. However, the pros and cons of the PVA are first noticed in the article of Boyce (1992), where the authors discuss necessary components to drive the PVA. Later many authors put a significant effort in providing a proper shape to the concept of PVA.
The scrutiny of this literature depicts that the role of PVA in predicting the species extinction status is still a debatable issue due to its level of accuracy. The studies of Taylor (1995), Brook et al. (1997), Reed et al. (1998), andLudwig (1999) show that the data deficiency and the presence of ambiguity over the population structure may sometime lead to the unacceptable results in the population persistence. On the contrary, Brook et al. (2000) performed a retrospective analysis using 21 wildlife populations and concluded that the outcomes of PVA are surprisingly accurate. Note that this study is being provoked to criticism due to the usage of short time windows. Apart from this, Schiegg et al. (2005) also showed the predictive accuracy of PVA by the population abundance dataset of red-cockaded woodpecker (Picoides borealis). The authors found that the predictability of PVA turned out to be accurate in a short period of 5 years. Consequently, Crone et al. (2013) discovers the poor accuracy level of PVAs due to the presence of variation in the environmental conditions in between the data collection and forecasting periods. So, it is recommended to include the environmental changes in the forecasting model structure for better accuracy. In this connection, several authors, viz., McCarthy et al. (2001), Münzbergová andEhrlén (2005), Pe'er et al. (2013), etc., suggest a list of methods to enhance the accuracy level of PVA. Note that, whatever be the degree of precision, the analysis of PVA is conducted based on the computer simulation.
IUCN has summarized five criteria (A to E) for evaluating the threatened status of many species. These criteria are based on the population reduction, the extent of occurrence, areas of occupancy, the decline of small populations, and number of mature individuals in declining tiny populations. Finally, a quantitative computer-based simulation is an essential tool for predicting the species threatened status.
The application of PVA is noticed in species conservation by evaluating the extinction probability of any species within a specified time window. This enumeration method in PVA depends on the two major components, i.e., (i) data availability and (ii) ecology of the species. Münzbergová and Ehrlén (2005) shows that the output of PVA can create a wrong impression on the species conservation due to the inappropriate collection of the datasets. In this regard, the authors propose a method describing that the sampling of the dataset at any stage of the species life cycle should be identical for the better accuracy of PVA outputs. In addition, the ecology of the species also plays a key role in providing a precise estimate of PVA, which the model selection criteria can well demonstrate. IUCN prescribed five different types of mathematical models to nurture the species extinction norms, viz., (i) occupancy model (Sjögren-Gulve and Hanski 2000), (ii) scalar dynamic (unstructured) model (Dennis et al. 1991;Burgman et al. 1993), (iii) stagestructured model (Akçakaya 2000), (iv) individual-based model (Lacy 2000), and (v) metapopulation model. Our entire work is designed on the scalar dynamic model structure. The dataset estimated by Sibly et al. (2005), where ̂ ∈ (0, 1) 188 259 447 The dataset, which are beyond the list of Sibly et al. (2005), where ̂ ∈ (0, 1) 15 21 36 Dennis et al. (1991) proposed a method of estimation of growth and extinction parameters for endangered species based on the stochastic exponential growth model. This seminal contribution is the first systematic and analytical study of population viability analysis using sophisticated statistical concepts and tools. However, most of the species growth rate is density dependent, so the choice of exponential model is hardly a likely event in reality. The same author has used an alternative statistical computing method, called data cloning, to calculate the maximum likelihood estimates and its standard errors for complex ecological models (Lele et al. 2007). Although the vital component of the method involves the Bayesian setup for MCMC computation, this framework is only used as a tool for the likelihood calculation. In contrast to the classical Bayesian analysis, the inferences are entirely invariant to the choice of the prior distributions.
Several software packages exist, viz., GAPPS, INMAT, RAMAS/AGE, VORTEX, etc., to explore the population viability analysis in recent times. However, there has not been any systematic attempt to compare predictions by several softwares based on a single dataset. There has been a significant change in the output of PVA analysis if we use the density dependence model for species growth with demographic and environmental stochasticity. Nevertheless, it is essential to note that the modeler has no idea what is going on in software packages for proposing the result of any study. This is typically known as the black box syndrome.

Extinction probabilities of the selected species
It is already mentioned in Sect. "Proposed simulation algorithm and probability of extinction" that based on the standard data filtration and simulation process, we choose 71 timeseries profiles for describing their conservation status. There are 5 Aves, 19 Insect, 15 Mammals, and 31 Fish data present among these profiles. We include our prediction results about this set of species in the Tables 2-7 of the supplementary file. These tables are formulated to serve two significant purposes. The first one is the enumeration of the probability of extinction of these species concerning 2021, which are enlisted in the fifth column of each table. The second one is the categorization of the species according to the IUCN guidelines. It is worthwhile to mention that at present there are five criteria in IUCN, i.e., criteria A, B, C, D, and E, respectively, to predict any species conservation status. All these criteria except the last one are generally enumerated based on the species' demographic characteristics, area of occupancy, habitat place, etc., which are available in many kinds of literature. However, criterion E is overlooked in most of the studies. The plausible reason for this ignorance is the quantitative analysis, which is the fundamental characteristic of criterion E. This conservation measure is also subdivided into three labels, i.e., Endangered (EN), Critically Endangered (CR), and Vulnerable (VU). The guidelines of these three labels are present in the IUCN report. We also marked some species with these labels and put them into the 8th column of Tables 2-7 present in the supplementary file.
Note that IUCN does not include any invertebrates like insects in the conservatory analysis. Recently, Fox et al. (2019) prepared a list of insects by labeling the conservatory position according to IUCN protocol. However, our enlisted insects in Table 3 (see the supplementary file) do not overlap with the set of species mentioned in the article of Fox et al. (2019). So, we prepare such kind of list in Table 3   Table 3 Conservation guidelines for the species with high extinction threat * We do not find the species abundance report of Metapolophium festucae in United Kingdom; ** The species is very sensitive to the environmental fluctuations so the species need special care but this task is not easy at all. A further thorough and deep study would be beneficial to identify the management actions for the conservation of the species. The conservation biologists can explore this in future a The fish Pink Salmon has the anadromous character, which is synergistic with the non-overlapping generation of the supplementary file for the Insect data following the IUCN guidelines. In our entire analysis, we observe 40% cases where our prediction regarding the conservatory label meets the IUCN result, and we observe 60% cases where it does not. We also make a column at the end of the table  to elucidate this scenario except for the aforementioned  Table 3. Here, we use two abbreviations, i.e., "M" for the species where our result matches IUCN and "MM" for the mismatch cases. It is worthy of mentioning that the primary objective of our research is to focus on the local extinction status of the species rather than the evolutionary extinction. Most of the previous studies by several authors (Saha et al. 2013;Bhowmick et al. 2015;Sau et al. 2020) and much other research focused on the local extinction issues of the species available in the GPDD database. Note that these previous elaborate studies for identifying the evidence of the Allee effect in the Herring population and its extinction status emphasized on the data of four localized zones such as the Baltic Sea, Atlantic Sea, and Canadian and Iceland regions. Here, the authors predict the time to extinction and the probability of extinction of this species in these three regions but not the evolutionary extinction of Herring fish.

Remark 8
The species abundance data in GPDD are pretty old to 2021. So, the prediction about the current conservatory status of any species would be more accurate if the time gap between the final calendar year of any time-series profile and 2021 is reduced. In this connection, we perform the cluster analysis of the 71 time series to explore several compartments of the difference in the time windows (DTW). The output of this analysis is present in the Supplementary File, where the GPDD IDs are present at the single branch of this dataset. The figure shows the different groups of DTW (b) (c) Fig. 7 The frequency distribution of the estimated probability of extinction to 2021, for the three major taxonomic groups, i.e., a Insect, b Fish, and c Mammal among which we select the cluster with a reasonable magnitude of DTW, i.e., 50. We also make a column in each of the Tables 2-7 of the supplementary file by providing the header "Reliability percentage" under which we keep the DTW value along with the deviation percentage from level 50. Note that we do not provide any reliability percentage in the Table 5 of the supplementary file due to the sizeable DTW magnitudes as the prediction process becomes cruder with the increasing DTW level.

Remark 9
We got the estimate of the noise intensity parameter ( ̂ ) for the aforementioned 483 time-series data. We draw the frequency histogram Fig. 8 of these estimated magnitudes according to the taxonomic class, i.e., Aves, Mammal, Fish, and Insect, respectively. The diagram reflects that the tolerability level of the noise intensity is not too high for these four groups of species. However, there are some variations in these tolerance levels. The figure shows that the taxonomic class Aves possesses the lowest tolerance level, and the Insect is the highest among these four. So, suppose we prepare a rank list about these tolerance levels. In that case, it should follow the trend of Insect, Fish, Mammal, and Bird respectively.
We enumerate the extinction probability of the selected 71 time-series profile by our proposed simulation method as mentioned in Sect. "Proposed simulation algorithm and probability of extinction". These extinction probabilities are attached at the fifth column of Tables 2-7 of the supplementary file. Figure 7 depicts the histogram of these probability magnitudes for three major taxonomic groups, i.e., Insect, Fish, and Mammals. Note that in the case of the selected birds, the probability values turn out to be 0; that is why we exclude the Aves group from the histogram Fig. 7 (see Table 2 in the supplementary file). The diagram exhibits that the extinction probabilities can be distinguished into two distinct zones like low and high, in the case of Fish and Mammal. However, for the Insect, we observe low-, moderate-, and high-risk zones about the probability magnitudes. It is worth mentioning that the high extinction probability certainly significantly raises any species' risk status (see Fig. 8). In this connection, we choose at most three species from the Insect, Fish, and Mammals with the high extinction probability to discuss their conservation management. But, in the case of the taxonomic class Mammal, we find only one such species. The taxonomic name of each selected species is present in the second column of Table 3. In order to assess the effectiveness of our prediction model in the conservation scheme of any species, we provide a detailed discussion on the population abundance of these selected species in several provinces of the world. This is given in the Supplementary File. Based on this analysis, we also prepare a column in Table 3, mentioning the conservation guidelines of each species. It is needless to say that the dataset of Metapolophium festucae in the GPDD is collected from the Hereford and Worcester provinces of the UK. However, we do not find any report about the abundance of this insect in the UK to support our findings. We also observe different PE values of any individual for the dataset of distinct zones. As a consequence, a variation is generated in the conservation status of any particular species. The taxonomic names of such species are Oncorhyncus nerka (sockeye salmon) and Solea vulgaris (Dover sole) with the GPDD ID 2019 and 2041 respectively. So, we exclude these two IDs from the third column of Table 3. Although sockeye salmon stands out to be the highly threatened species among several places of British Colombia and Alaska, the extinction probability has been 0 for the dataset of Stellako River. The survey of Welch et al. (2007) on the Stellako River support this fact. Similarly, in the case of Dover sole, we observe a high extinction level at the ICES VIID, North Sea region. Nevertheless, the species is abundant in the Bay of Biscay. Koutsikopoulos and Lacroix (1992) have a closed parity with our findings for the region Bay of Biscay. So, the predicted probability of extinction for the species Solea vulgaris is naturally low (0.003) in this region. The plausible reason for this kind of disparity is the estimation of the local extinction of many species.

Conclusion
The study of Sibly et al. (2005) reveals that most of the species follow the theta-logistic growth trait with the convex downward trend. This entire research work is formulated based on the deterministic setup. In contrast, the involvement of external noise in any growth system is inevitable. Here, we frame the theta-logistic model with the stochastic analog in two directions, i.e., the discrete and continuous setups. The discrete model analysis provides that the extinction state's attainment happens more quickly in the case of a stochastic setup than the deterministic version. Although the role of chaos in species extinction is debatable, literature surveys suggest that chaos with stochasticity accelerates the extinction of species. Similarly, in the case of the continuous version, we performed a theoretical study on the stochastic theta-logistic model to provide a critical value of the noise intensity parameter.
This environmental bound will act as the sustainability criteria of any species environmental tolerability. In this connection, we use the data of four major taxonomic groups, i.e., Bird, Insect, Mammal, and Fish, from the GPDD database and classify the species based on the environmental sensitivity, i.e., low, moderate, and high. The highly sensitive species have a low tolerance level, associated with the small magnitude of environmental noise intensity parameter . Moreover, the simulation prediction model on these four taxonomic classes also shows that the overall extinction probability of the Bird is almost negligible for the current time window.