Next Article in Journal
Feedback Control Techniques for a Discrete Dynamic Macroeconomic Model with Extra Taxation: An Algebraic Algorithmic Approach
Next Article in Special Issue
Modeling Study of Factors Determining Efficacy of Biological Control of Adventive Weeds
Previous Article in Journal
Cooperative Guidance Strategy for Active Spacecraft Protection from a Homing Interceptor via Deep Reinforcement Learning
Previous Article in Special Issue
Survival Analysis of a Predator–Prey Model with Seasonal Migration of Prey Populations between Breeding and Non-Breeding Regions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Autoregression, First Order Phase Transition, and Stochastic Resonance: A Comparison of Three Models for Forest Insect Outbreaks

1
V.N. Sukachev Institute of Forest SB RAS, Krasnoyarsk 660036, Russia
2
Krasnoyarsk Scientific Center SB RAS, Krasnoyarsk 660036, Russia
3
Institute of Biophysics SB RAS, Krasnoyarsk 660036, Russia
4
Department of Ecology and Nature Management, Siberian Federal University, Krasnoyarsk 660041, Russia
*
Author to whom correspondence should be addressed.
Mathematics 2023, 11(19), 4212; https://doi.org/10.3390/math11194212
Submission received: 12 September 2023 / Revised: 5 October 2023 / Accepted: 6 October 2023 / Published: 9 October 2023

Abstract

:
Three models of abundance dynamics for forest insects that depict the development of outbreak populations were analyzed. We studied populations of the Siberian silkmoth Dendrolimus sibiricus Tschetv. in Siberia and the Far East of Russia, as well as a population of the pine looper Bupalus piniarius L. in Thuringia, Germany. The first model (autoregression) characterizes the mechanism where current population density is dependent on population densities in previous k years. The second model considers an outbreak as analogous to a first-order phase transition in physical systems and characterizes the outbreak as a transition through a potential barrier from a low-density state to a high-density state. The third model treats an outbreak as an effect of stochastic resonance influenced by a cyclical factor such as solar activity and the “noise” of weather parameters. The discussion focuses on the prediction effectiveness of abundance dynamics and outbreak development for each model.

1. Introduction

Analyzing processes related to forest insect population dynamics is crucial for two reasons. Firstly, the simplicity of population estimation opens up opportunities to study patterns of population dynamics and the factors influencing these dynamics. Secondly, forest disturbances are increasingly common worldwide, and they are growing in frequency and intensity [1]. These disturbances have significant impacts on both local and global carbon stocks in terrestrial ecosystems [2]. Insect damage to forests can result in forest mortality and trigger large-scale forest carbon emissions [3]. For example, the bark beetle Dendractonus frontalis damaged more than 15 million ha of forest lands in the western United States [4]. Attacks by the bark beetle have resulted in the destruction of 10 to 20 Tg C year−1 of living biomass [5,6] and have caused an associated decline in NEP, estimated at about 3 Tg C per year from 2000 to 2009 [5].
In Canada, major bark beetle infestations affecting 12 million ha were recorded in the 2000s, killing about 20 Tg C year−1 of living biomass, resulting in an estimated average annual reduction in NEP of 13.5 Tg C year−1 [7].
Large-scale carbon balance impacts from defoliation by the gypsy silkmoth (Lymantria dispar L.) have affected more than 5 million ha of forests in the northeastern United States since the 1980s [4]. Defoliation by the gypsy silkmoth resulted in a 41% reduction in NEP in mixed hardwood stands in New Jersey between 2005 and 2007 [8]. Regionally, NEP reductions can also last for decades when significant tree mortality occurs [7,9,10,11].
In the territory of Russia, in 2021, signs of forest damage were detected in an area of 76,806 thousand hectares [12].
Therefore, it is of high importance to develop forecasting methods for modeling the population dynamics of pests and their outbreaks, which would allow to take appropriate counteractions to reduce the pest population and its impact on the forest.
The analysis of processes in complex systems, specifically those involving insect populations, should follow a systems approach. It is necessary to analyze all components of the system with respect to the insect species under study [13,14,15,16,17,18]. Usually, to describe such a system, since the time of A. Lotka and V. Volterra, differential equation systems are used, where each equation describes the dynamics of a separate component [19,20,21]. However, the system analysis postulates fall apart upon the initial encounter of a researcher with the object of study, specifically a population of a forest insect species. While theoretical models of population dynamics can incorporate any desired variables, measuring relevant community characteristics poses a challenge in real forest insect population modeling [22,23]. Long-term forest measurements typically only allow for the registration of population densities [24,25]. Accounting for system components such as parasite and predator populations, as well as the volume and quality of available forage resources, proves to be exceedingly difficult. Thus, research opportunities on the dynamics of the studied forest insect species are limited in the field, and the system as a whole “escapes” from system model analysis. Furthermore, incorporating modifying factors, such as weather, as separate equations in the system is impossible due to the self-generation of individuals in the population that would occur as a result. Describing the impact of weather on the coefficients of a system’s equations requires introducing further equations that account for weather factors. As these additional equations are nonlinear, they contribute to increased model system dimensionality and difficulties in obtaining analytical solutions [26].
There are two types of population dynamics observed with forest insects [27,28]. The prevailing type features stable dynamics with low population densities in most species of insects. However, a few insect species, such as the pine looper Bupalus piniarius L. and Siberian silkmoth Dendrolimus superans Tschetv. (which is a dangerous pest of Siberian taiga forests), can be characterized as a bistable system, which can be in one of two phases: the stable dispersal phase and the outbreak phase. In the stable dispersal phase, the density x1 of the population is extremely low and its individuals utilize only a small part of the available food objects. The outbreak phase occurs only when the pest population density surpasses a critical value (xr), eventually reaching the value x2, where the entire population consumes all available food objects in the given territory [28,29,30,31,32], and finally returns from state x2 to state x1.
In this paper, we aim to develop minimalist models of forest insect population dynamics based only on density data of the studied species, without using the properties of regulating and modifying factors. To model a crucial event such as a population outbreak of forest insects, we consider the outbreak on the one hand as an autoregressive process and on the other hand as a bistable system, where the outbreak is viewed as an analog of phase transitions of the first kind (TPI) in physical systems, and finally as a realization of the effect of stochastic resonance (SR) in the ecosystem.

2. Materials and Methods

Siberian silkmoth Dendrolimus sibiricus Tschetv. is a species of the family Lasiocampidae. The caterpillars feed on needles of almost all coniferous species found within its range. It prefers larch Larix sibirica LDb, often also damaging fir Abies sibirica Ldb. The species is widespread in Russia within the Urals, Western and Eastern Siberia, and the Far East. Outside Russia, the Siberian silkmoth is distributed in Mongolia, Kazakhstan, Korea, and northeastern China. The southern border of its range is at 40° N.
The Siberian silkmoth is the most dangerous pest of coniferous forests in Siberia and in the Russian Far East [33,34,35,36,37,38,39,40,41,42]. For instance, during the last outbreak (years 2014–2017), forest stands of Siberian fir Abies sibirica Ledeb. were damaged in the Krasnoyarsk region in a territory of 800 thousand hectares.
The high efficiency of Siberian silkmoth outbreaks in dark coniferous forests is determined primarily by the unusually high density of its populations. The maximum infestation of fir stands reaches more than 20 thousand caterpillars per tree. However, in overconsolidated centers, the pest starts to lack food. The reproduction rate decreases due to the aggravation of competitive relations. Fecundity (1.5 times) and survival rate (2.5 times) become significantly lower than in the previous phase.
In the depression phase, when the population density sharply decreases, the extremely harsh conditions of the Siberian silkmoth’s existence during the depression period lead to a minimal population survival rate (less than 1%) [26].
The total duration of the gradational cycle of the Siberian silkmoth in the dark coniferous southern taiga forests of the Krasnoyarsk region averages 14 years. The phase of increasing numbers in the population has the longest time interval (4–7 years), which is due to the manifestation of long and complex processes of “breakaway” from natural enemies and reorganization of intrapopulation mechanisms. The maximum phase lasts no more than two years, during which intensive stand damage is observed. The final stages of the outbreak (thinning and depression) are transient (two to three years), and the duration of recovery and stabilization is from two to six or more years [26].
The dynamics of Siberian silkmoth populations is largely regulated by a complex of entomophages. In Siberian silkmoth reserves, they affect 60–90% of eggs, 10–40% of caterpillars, and 50–70% of pupae [39]. During the phase of increasing insect numbers, the proportion of affected eggs varies from 16.6 to 54%, caterpillars from 3.6 to 10.9%, and pupae from 4.8 to 22.3% [36]. During the climax of the outbreak, parasites destroy 70–90% of eggs, 65–95% of caterpillars, and 25–60% of pupae [34]. However, practically all parasites of the Siberian silkmoth are affected by superparasites, and specialized entomophages are most intensively affected, e.g., eggeater Telenomus gracilis Mayr, caterpillar parasite Rhogas dendrolimi Mats., pupal parasite Masicera zimini Kol. [26,36]. This significantly increases the inertia of the complex of natural enemies and reduces their regulatory effect.
The complexity of monitoring and predicting the abundance dynamics of this species is that outbreaks can occur in some zones in a remote sparsely populated territory of 20 million square kilometers in total; as a consequence, for over more than a century of observations, researchers have managed to obtain only two fairly long time series of the dynamics of this species [37,41].
Pine looper B. piniarius is a dangerous pest of pine Pinus sylvestris L. The caterpillars of the pine looper damage the needles of the Scots pine. Adult insects appear in the pine stands from the last week of May (southern), but mostly begin in June (middle) and in late June to early July in the northern areas. Larvae undergo four ages and complete their development in late summer or early fall. Pupation usually takes place in October. Outbreaks of this species occur in the pine forests of Europe (in Scotland, England, the Netherlands, Germany, Poland, Sweden, and Latvia) and in forest stands in Russia (the Urals, the Altai, and South Siberia) [42,43,44,45,46,47,48,49,50,51,52,53,54,55,56]. The abundance dynamics series of this species in Thuringia, as published by F. Schwerdtfeger (1939, 1952, 1968), is one of the longest series (80 years) known in forest entomology.
For these species, three types of models were examined to describe their abundance dynamics.

2.1. The Autoregression (AR) Model

Traditionally, during a description of population dynamics in ecosystems, oscillatory processes have been analyzed in terms of the Lotka–Volterra model, and the presence of fluctuations has been regarded as a consequence of interactions between populations of a prey species and of species of predators and parasites [57]. In our opinion, however, it is impossible to use such models to describe real outbreaks of insects because data on the density of populations of parasites or predators are usually unavailable to the researcher. Nonetheless, because a population of any animal species can be considered an autoregulating system with feedback loops, this kind of general model can be described in terms of a so-called transfer function with feedback loops of different strengths and directions and with different delays [58]. In this context, during an analysis of general properties of temporal dynamics of animal populations, the factors that affect a given population can be ignored, and it is sufficient to assess (i) the susceptibility of the population to the influence of the feedback loops and (ii) the delay of the response of the system to its state at the start [59,60,61]. In this case, if the aim is to build a model of a prey species, then the verification of such a model does not require information about the abundance of predators. To construct an AR model of abundance dynamics of specific animal species, it is necessary to evaluate the sign and sensitivity of the feedback loops.
For the systems under study, it is assumed here that the analyzed feedback loops are linear, depend on the population density y(t), and have the following form: w ( t ) = a x ( t τ ) , with delay τ and feedback amplitude a > 0 for a positive feedback loop and a < 0 for a negative feedback loop. The number of feedback loops having different signs and different delays is unknown for specific species, but, in the general case, taking into account feedback loops for a temporally stationary AR system, the following equation can be written:
y ( t ) = j = 1 k a j y ( t j ) + μ
where y(t) is the density of insect population; k is the order of the AR; μ is the measurement error; and a0, …, aj, …, ak are coefficients characterizing the intensity of feedback loops having delay j.
Quantity a j = y ( t ) y ( t j ) in Equation (1) characterizes the susceptibility of the population in year t to changes in its abundance in year (tj). In this context, the influence of feedback loops at time point t is proportional to variable y(tj), and the type of feedback (positive or negative) is determined by the sign of aj. If aj is positive, then this parameter characterizes a linear positive feedback loop between years t and (tj). If aj is negative, then aj characterizes a linear negative feedback loop between abundance values in years t and (tj).
Therefore, to build a model, it is necessary to estimate order k of AR. k characterizes the number of feedback loops in the system and signs and values of coefficients in Equation (1). A partial autocorrelation function (PACF) is usually employed to estimate the order of AR and a delay in a feedback system. AR order k in the system is estimated from the maximum value of k for which the PACF is significant [62]. If time series {y(t)} and order k of the AR in Equation (1) are known, then this equation can be treated as a linear regression equation with respect to unknown parameters a0, …, aj, … ak. In this case, values of the feedback coefficients in the AR equation can be found by analyzing Equation (1) as a linear regression equation with unknown variables, and then its solution can be found using traditional least squares methods [63].

2.2. Models of an Outbreak as a First-Order Phase Transition

Taking into consideration the challenges surrounding surveys on population density, it appears that models are necessary to describe the processes occurring in populations that use only minimal information about the state of the system and a minimal number of variables and free parameters in the model.
In this paper, to describe the processes associated with the development of population outbreaks of phyllophagous insects, we propose to use ecological analogs of the model of first-order phase transitions in physical systems. A classic example of a first-order phase transition in physical systems is, e.g., boiling (the transition of a substance from a liquid state to a gaseous state that occurs when the critical temperature Tc of the environment is reached [64]).
According to the phenomenological model, phase transitions in physical systems are described via thermodynamic potential G, and condition G min is an optimization principle that regulates the phase transition process [65,66]. Stable states of a population are characterized by a minimum (local or global) of potential function G(x). Potential function G(x) for a bistable system with two stable (or metastable) states is characterized by the presence of two local minima of G(x) values (potential pits).
For an ecological system, such a potential function can be estimated using an approach based on investigation into the statistics of the times the system spends in various states [67]. With this approach, to describe the model of an outbreak, let us introduce the concept of potential function G(x) of the state of the system, which we define as the reciprocal of the density function of distribution f(x) of density values x of the population for the long period T, sufficient for all possible states of the system to materialize. Then, possible states of the system can be described by means of potential function G(x), which is the inverse of probability p(x) of the state of the system with density x of the population:
G ( x ) = 1 p ( x )
If the probability of occurrence of state xm during a sufficiently large observation time T is small, then the value of potential function G(xm) will be large. If, during period T, state xm occurs often enough, then G(xm) is small.
The value of G(x) potentially depends on a large number of various factors, and it is difficult to express the dependences of G(x) on these factors. Nevertheless, G(x) for a population can be classified depending on the number, values, and positions of local minima and maxima of this function.
Function G(x) can have one global or several local minima. Potential function G(x) with one global minimum at x = x1 will characterize a system with one stable state. The existence of the system near stable state x1 is associated with the implementation of negative feedback loops in the system when its state deviates from x1. A bistable potential function of the state with two local minima G(x1) and G(x2) (where x1 << x2) and one local maximum G(x12) between the local minima will characterize a system in which switching from state to state is possible. If G(x12) is large, this means that the state in question is observed very rarely and that the system goes from state x1 to x2 and back, quickly skipping state G(x12), which will be observed quite rarely.
By analogy with physical models of first-order phase transitions, we will examine ecological phase transitions characterized by jumps of the “forest–insects” system from a state with low population density and low damage to trees in a forest stand to a state with high density and severe damage to the forest stand, and back. The more x(t) deviates in absolute value from the x1 value at some time point t, the less is the probability of reaching this value and the higher is the probability that the “forest–insects” system in this state will not exist for a long time.
From the standpoint of the theory of phase transitions, external factors affecting the insect population can be considered analogous to external fields.
A general shape of the potential function and its change under the influence of an external field are depicted in Figure 1.
At point x = xc, function G(xc) is G0, which is greater than the values of G(x1) = G(x2) at points x1 and x2. The G0 value characterizes the height of the potential barrier between the two stable states. Values of function G(x) near x1 (or near x2) are usually described as a potential well, and the depth of the potential well is defined as the difference between G0 and G(x1) or G(x2). Point x = xc is the point of the phase transition from the state with minimum population density x1 to the state with maximum population density (the state of an outbreak), i.e., density x2. Thus, we can say that the shape of the empirical potential function describing abundance dynamics of forest insects can be explained using the notion that population processes are first-order phase transitions, manifesting themselves as large jumps in population density and subsequent rapid fluctuations around stable states.
Function G(x) in Figure 1 can be described with the following parameters:
  • Local minima of potential function G(x1) and G(x2) and corresponding values x1 and x2 of population density (in this context: x1 << x2);
  • Local maximum of function G(xc): the height of the potential barrier when potential function G(x) reaches the local maximum;
  • With the help of these “base values” additional parameters can be defined;
  • Difference Δ x = x 2 x 1 , which characterizes the range of population densities.
  • Difference Δ G = G ( x 2 ) G ( x 1 ) between values of potential functions G(x2) − G(x1);
  • Depths of potential wells: Δ G 1 = G ( x r ) G ( x 1 ) and Δ G 2 = G ( x r ) G ( x 2 ) ;
  • Half-width of the potential barrier x 1 c at which G ( x 1 c ) = G ( x c ) 2 is the half-height of the potential barrier;
  • Absolute values of derivatives d G d x of the potential function to the left and to the right of point x = xc. These derivatives can be viewed as the susceptibility of state function G(x) to a change in density x of the population. At sufficiently small values of d G d x , we will say that potential G(x) is “soft”; at large values of d G d x , the potential function will be characterized as “hard”. Approximately, the values of the derivatives can be replaced by quantities G ( x c ) G ( x 1 ) x c x 1 and G ( x c ) G ( x 2 ) x 2 x r .
  • Current susceptibility χ = x ( i + 1 ) x ( i ) (or on a logarithmic scale: χ ln = ln x ( i + 1 ) ln x ( i ) ) characterizes the current reproduction rate of the population.
The influence of external factors within the model of first-order phase transitions for the “forest–insects” system is taken into account as an analog of some external field h (for example, weather parameters). If external field h > 0 increases, then the depth of the left local maximum in Figure 1 will decrease in proportion to the magnitude of the external field, and, at some time point tr, the amplitude of fluctuations of normalized population density at low population density will become comparable or even greater than difference G ( x c ) G ( x ( t r ) . Under this scenario, the system will switch to the state with density x2 and an outbreak will occur (curve 3 in Figure 1). A change in the sign of the external field and/or an increase in the impact of regulatory factors (for example, an increase in the pressure of parasites and/or a decrease in the volume and quality of available food) will cause a reverse switch in the population, i.e., cessation of the outbreak and a return of the system to the stably rarefied state.
A distinctive feature of phase transition models is that they do not depend on time [68]. Such an assumption substantially simplifies the model. Now, to describe the system, it is necessary to find only stable minima of function G. Of course, with the representation in the form of a potential function, information about temporal dynamics of population density is lost. Nevertheless, the totality of the main and additional parameters that thoroughly characterize the shape of potential function G(x) gives us clues to the dynamics of transitions between the phases of the system.
Thus, to build a model of insect abundance dynamics, it is necessary to evaluate the shape of the potential function of an insect species by means of long-term data on population density both in a stably rarefied state and in an outbreak phase, to determine the influence of weather factors on pest abundance dynamics, and to evaluate the nature of population density fluctuations near stable states.

2.3. Stochastic Resonance (SR) Model

This model manifests itself in nonlinear systems having different characteristic frequencies when the system is subjected to a weak periodic perturbation and external noise is added [69,70,71,72,73]. The amplitude of the periodic external perturbation is assumed to be so small that it rules out passage through the barrier in the absence of noise. Nonetheless, a small periodic external perturbation leads to small modulation of the height of the potential barrier. At low noise intensity, the average transition times are quite long and greatly exceed the periodicity of the modulation signal. The term “resonance” is used in the SR model because the switch from one state to another depends on the amplitude of the noise factor. The highest probability of switching from state x1 to state x2 is observed at some amplitude Aopt, and at amplitudes much larger or smaller than Aopt, the probability of the jump is much less [67]. Accordingly, at a certain noise level in a stochastic bistable system, a modulating signal will be amplified maximally.
Population density dynamics in the two-well potential U ( G ) = a x 2 2 + b x 4 4 under the action of white noise ξ(t) having intensity D and periodic force f ( t ) = A cos ( ω t + ϕ ) are described by the following equation [67]:
d x d t = x x 3 + A cos ( ω t + ϕ ) + 2 D ξ ( t )
SR can produce intriguing and counterintuitive dynamic effects in living systems. Many authors have stated that noise acts not only as a source of disorder but also as an integral feature of the dynamics of natural systems, including a system of insects [74,75,76,77,78].
What kind of specific mechanism of SR can exist for an insect population? It can be hypothesized that the role of an external modulating signal can be played by geophysical disturbances described by solar activity cycles and characterized by means of the Wolf number (http://www.wdcb.ru/stp/solar/sunspots.ru.html, accessed on 12 September 2023). Solar activity has been repeatedly proposed as a generator of outbreaks, and weather fluctuations during a season may act as external noise.
The hypothesis of a direct or indirect link between solar activity and animal population dynamics, proposed over a century ago [79,80], has been tested repeatedly [81,82,83,84,85,86,87,88,89,90,91,92]. Hypothesis testing frequently was reduced to estimating cross-correlation functions between time series of Wolf numbers and time series of population dynamics of selected animal species. However, this approach raised the problem of explaining the lack of synchronization of local population dynamics at the planetary level in a situation where they should have been regulated by global solar activity.
Nevertheless, because solar activity is planetary in nature, if insect abundance dynamics are modulated by solar activity, then cycles of abundance dynamics of all species of forest insects on the planet should be synchronized with cycles of Wolf numbers, and outbreaks of all species should occur simultaneously. This idea, however, is not supported by observations. By contrast, under the SR regimen, population outbreaks will occur with the periodicity of solar activity only where the noise level is such that it ensures the implementation of SR. As soon as noise intensity on this territory starts to deviate from the level optimal for the manifestation of SR, the cyclicity of outbreaks will be disrupted. Thus, SR can lead to a “smearing” of outbreaks in time and space, which is observed in the wild.

3. Results

3.1. AR Models of Dynamics of Forest Insect Abundance

At the first stage of the construction of the AR model for abundance dynamics of the Siberian silkmoth D. sibiricus Tschetv., the data series was transformed.
Steps of the transformation of the dynamics series of Siberian silkmoth abundance x (individuals per tree) included switching x ( i ) ln x ( i ) from the population density scale to the logarithmic scale. At the next step of the transformation, high-frequency filtering of series {ln x(i)} was performed with the help of the Hann filter [93]:
L ( i ) = 0.24 ln x ( i 1 ) + 0.52 x ( i ) + 0.24 x ( i + 1 )
Next, order k of AR of series {L(i)} was assessed using a PACF (Figure 2).
As shown in Figure 2, the order of the AR model for transformed series {L(i)} can be assumed to be 2. Then, the AR equation of transformed series {L(i)} of the Siberian silkmoth abundance dynamics can be written as
L ( i ) = a 0 + j = 1 2 a j L ( i j )
where a0–a2 are some coefficients.
Because variables L(i) in Equation (5) are known, this equation can be considered a linear regression equation with respect to the unknown parameters a0–a2. The solution of such an equation by the least squares method is well known, and Table 1 presents values of parameters a0–a2 for this population.
As one can see in Table 1, interactions between different generations of this species from the standpoint of regulation theory can be characterized as a positive feedback loop between adjacent generations and a negative feedback loop between ith and (i − 2)th generations.
In Figure 3, values of the transformed series and model series are given for abundance dynamics of the Siberian silkmoth in the Lower Yenisei region.
To assess the quality of the model, four criteria are utilized: coefficient of determination R2 has to be close to 1.0, values of the coefficients in regression Equation (5) have to be significant according to the t test, the transformed survey data series and the model series have to be synchronous (the synchronicity is evaluated using characteristics of a cross-correlation function (CCF) between these series), and a series of residuals between values of the survey data series and the model series have to be uncorrelated and to have properties of white noise. The latter condition is assessed, for example, by the Dickey–Fuller test [94]. If the model meets all these criteria, we can say that it adequately describes the dynamics of the population under study.
As presented in Table 1 and Figure 3, the accuracy of model 5 is very high: coefficient of determination Radj2 of the model is 0.96, and the model’s coefficients are significant at p < 0.02.
For the Siberian silkmoth population in the Russian Far East, order k of AR, just as in the Siberian population, is 2 (see Figure 2). Table 1 lists the coefficients and their significance for the Siberian silkmoth model in the Far East, and Figure A1 presents curves of dynamics of the transformed field data and model calculations.
As one can see, after a comparison of data on abundance dynamics of the Siberian silkmoth in Table 1, coefficients of the models for these two populations are quite similar both in the sign and in absolute value. The explanation for the positive feedback loop between the densities in years i and (i − 1) is extremely simple: the bigger the parent generation, the larger the number of offspring.
The negative feedback loop between abundance levels in years i and (i − 2) can be attributed to a delay in the influence of parasites and predators on the prey population [61].
For AR series with order k = 2, the cyclic behavior is more characteristic, and the spectrum of such series contains maxima [63]. Spectra of the series for the Far East and the Lower Yenisei {L(i)} are displayed in Figure A2. As readers can see in the figure, the frequency of population outbreaks of this pest in both taiga forests near the Lower Yenisei and in forests of the Far East is 11–14 years. However, the spectrum of the abundance series is quite wide, indicating high variance in the pest abundance dynamics and making simple prediction of the development of outbreaks problematic.
Thus, to build models of population outbreaks of the Siberian silkmoth by means of only data from population surveys, data from the counts for 6–8 years are needed. It is practically impossible to obtain such data in Siberian taiga as a whole, and the presented model provides an understanding of the processes underlying the dynamics but does not allow for a practical forecast of outbreak risks.
To construct a model of abundance dynamics of the pine looper, the same sequential procedures are performed: series filtering, calculation of the PACF (for this species, the AR order is 3; Figure 2), computation of model coefficients, construction of abundance dynamics series for the transformed and model data, and calculation of the CCF and of the spectrum. To compute the model parameters, we use survey data between years 1881 and 1930. Then, from the coefficients of the model (Table 1), we calculate model data on density for the 1931–1939 period, and these data are compared with transformed field data (Figure A3). The proposed model makes it possible to accurately predict population density for 9 years into the future (curve 3 in Figure A3).
The highest spectral density for the analyzed series from the pine looper is seen at frequency f = 0.086 (with periodicity T = 11.6 years) (Figure A2).
As readers can see, very good agreement between the field and model data is observed for the series under study. Thus, the development of an outbreak according to the AR model is determined by the influence of regulatory factors, and the systems investigated here have long-term memory, which is what causes the cyclic fluctuations of population densities.

3.2. The Model of First-Order Phase Transitions for the Siberian Silkmoth Population in the Lower Yenisei River Region

Using the data from the transformed series and from the high-frequency component-free time series of abundance dynamics of the Siberian silkmoth in the Lower Yenisei River region, potential function G(ln x) is constructed (Figure 4).
It is evident that the empirical state function is in good agreement with the theoretical model, is bistable, and has two local minima ln x1 ≈ −1 and ln x2 ≈ 5 and a potential barrier near density ln xc ≈ 3. A similar shape is shown by potential functions for the Siberian silkmoth population in the Far East and for the pine looper (Figure 4).
A question arises: is it possible to estimate the risk of a phase transition from a stably rarefied state to an outbreak state using the characteristics of the state functions? The risk of an outbreak may be related to extreme values of function G(x), e.g., local minima and maxima. In Table 2, characteristics of the state function are given for the analyzed phyllophagous populations.
From the data in Table 2, it follows that the type of potential function differs among different populations of phyllophages. The hardness of the potential function is high for the Siberian silkmoth both in the Lower Yenisei region and in the Far East, i.e., with sufficiently small increases in population density, the probability of reaching the critical population density for this species will be maximal, and, even with a sufficiently small rise of abundance, the risk of an outbreak will be high. For the pine looper, the potential function can be regarded as soft.
To assess the relation between population density and the risk of a pine looper outbreak in the 1880–1940 period, let us examine Figure 5, where the horizontal axis is the logarithm of current density ln x(i), and the vertical axis is susceptibility χ ln .
As illustrated in Figure 6, the cluster with the upper left and lower right boundaries 3.5 3 0 0 , which characterizes the state of the population immediately before an outbreak, is quite compact, and when the population reaches this state, an outbreak will inevitably occur.
Another derived parameter that can be employed for outbreak risk assessment is semicritical density ln x1c, i.e., population density when state function G(ln x1c) is equal to half a maximum of G(ln xc) (Table 2). At the semicritical density value for the Siberian silkmoth (2.3–2.5), the probability of reaching a critical value at which an outbreak can occur is twice as high as the risk of reaching state ln xc.
Fortunately for taiga forests, the Siberian silkmoth cannot be considered the most dangerous species in terms of the proposed risk index. Rather, this insect can be classified as a species with a rather low probability of outbreaks, as confirmed by rather long average times between adjacent outbreaks of 11–14 years [36,39]. It is possible that the absence of the most dangerous species is due to evolutionary selection. An extremely dangerous species is able to destroy its food base and then will itself be eliminated from the ecosystem.
Thus, the risk of an outbreak can be linked to a population’s susceptibility and semicritical density near the potential barrier. If the susceptibility of a population varies little over time, then, other things being equal, the current risk of an outbreak is related to the ratio between current and semicritical densities. The closer the current population density is to the semicritical density for a given species (and the greater its susceptibility), the higher the risk of an outbreak at the current population density, especially with warm weather acting as an external field.

3.3. An Outbreak as a Consequence of SR

For bistable systems, SR results from the combined action of a cyclic factor and random noise. For the analyzed populations of forest insects, the factors leading to SR are solar activity (characterized by a series of Wolf numbers W (or series dW of first-order differences of Wolf numbers) for a given year) and fluctuations of the environmental temperature (described by a stationary broadband series of first-order differences of adjacent values of air temperature throughout a year). Resonance will be observed at certain levels of the power of the noise signal. In this case, SR is likely present if CCF(k) between the abundance dynamics series and the Wolf number series is close to 1.0 at shift k = 0 between these series (Figure 6).
Can we regard the series of first-order differences of air temperature for the pine looper in Thuringia (according to the weather station closest to the source area in Halle, Germany (+51:30:53, +011:57:02)) as noise?
Figure 7 presents this time series.
The calculated spectrum of time series is broadband, and as a generalized characteristic of the time series of first-order differences, total spectral power S = 0 0.5 s ( f ) d f was next used.
Figure 8 shows combined data on Wolf numbers W(t), spectral power S(t), and transformed density L(t) of the moth population in Thuringia.
As shown in Figure 8, the spectral power of the external noise factor is stationary in time, and series W(t) and L(t) are synchronous (see the CCF in Figure 6).
For the Siberian silkmoth population, instead of the time series of Wolf numbers W, series of first-order differences dW of the Wolf numbers are employed. Although data on current population density of the Siberian silkmoth before 1954 and after 1968 are unavailable, the years when outbreaks occurred are known [36]. From these data, it is possible to build a dot series of outbreak years and a series of first-order differences dW of Wolf numbers (Figure 9).
If we assume that SR manifested itself over the entire period from 1954 to 2020, then population outbreaks of the Siberian silkmoth should have been observed in the years 1978, 1988, 1998, and 2011 (Figure 9). In actuality, outbreaks were observed near the Lower Yenisei in the years 1998 and 2014, and two outbreaks that were possible according to parameters of dW (years 1978 and 1988) did not take place. Why did this happen? If we compare the mean of noise S = 8.33 and error σx = 0.30 during 1952–1974 with the mean of S = 7.27 and error σx = 0.24 of noise during 1975–1993, then, according to the t test, we will obtain the t value of 2.74, significant at p < 0.01. This finding suggests that the power of noise in the years when an outbreak materialized is significantly greater than the power of the noise in the years when, despite the presence of maxima of dW, outbreaks did not occur. This finding is consistent with a parabolic dependence of SR on the power of noise [67].
For the Siberian silkmoth population in the Far East, the series of Wolf numbers and the series of abundance dynamics are in antiphase (Figure 6). Nevertheless, series S(t) is stationary, and for the entire observation period, the dates of outbreak maxima coincide with the dates of Wolf number minima.

4. Discussion

The ecological mechanisms causing population outbreaks are similar among all three models. Positive and negative feedback loops that are responsible for the emergence of an outbreak according to the autocorrelation model, are responsible in the first-order phase transition model for entering the zone near the potential barrier (positive feedback) and for returning to state x1 (negative feedback). In the SR model, the combined effects of the cyclic impact of solar activity and of stochastic weather factors lead to a decrease in the height of the potential barrier and to the switching of the system from a stably rarefied state to an outbreak state. However, a definite correlation between solar activity and insect dynamics does not lead to success, because the possible relationships between solar activity and population dynamics are more complex and, as shown in the text, it is necessary to include additional variables in the analysis, and it is also necessary to take into account the possibility of phase shifts between the studied series. Additional research is required to determine the causes of such phase shifts.
Nonetheless, to predict outbreaks by means of these models, it is necessary to perform the monitoring using diverse parameters. Surveys of the abundance of pine looper populations in pine stands that are close to the park type are quite simple, but an assessment of population density of the Siberian silkmoth in taiga forests is extremely difficult. Furthermore, the models do not allow to determine where specifically to expect an outbreak. The impact of solar activity is global, whereas outbreaks are local, and it is not clear where exactly an outbreak will occur in the next period favorable for insect development. Previously, the present authors have proposed approaches that can help to identify insect-susceptible zones by means of remote indicators of the state of forest stands [95]. It is possible that combined use of abundance dynamics models and data from remote probing will enable investigators to solve the problem of outbreak forecasting.
Although there is likely to be a considerable effect of solar activity on abundance dynamics of the species of forest insects under study, for each of the three studied populations, the type of the relation between outbreaks and solar activity is distinct. For the pine looper, dates of outbreak maxima are synchronized with dates of Wolf number maxima; for the Siberian silkmoth in the lower reaches of Yenisei, dates of outbreak maxima are synchronized with dates of maxima of first-order differences in the Wolf number series; and for the Siberian silkmoth in the Far East, dates of outbreak maxima are in antiphase with dates of maxima in the Wolf number series. The reason for these differences is still unclear. Apparently, a further analysis is needed involving data on abundance dynamics of other species. It is probable that some characteristics of host plants are responsible for the phase shift between the dynamics of solar activity and abundance dynamics of individual insect populations; these plant characteristics are not taken into account in the models presented here owing to the lack of such data. Nevertheless, it can be hypothesized that these local phase shifts lead to the observed desynchronization—at the planetary scale—between the indicators of solar activity and abundance dynamics of local animal populations.

5. Conclusions

The three population dynamics models presented are not exclusive to each other but instead reflect distinct facets of insect function, each characterized by unique analysis and forecasting capabilities for population dynamics and outbreak development.
The AR model is able to describe the population dynamics and predict the development of outbreaks in situations where population counts are carried out over a long period of time in a given habitat. However, due to the limited number of surveys in the taiga area of 2.7 million km2, the AR model is valuable for its description of the positive and negative feedbacks in the forest ecosystem, which is crucial for understanding the fundamentals of the population dynamics of forest insects.
The stationary potential function of the first-order phase transition model allows for assessing outbreak risk in a specific locality with population data from two consecutive years.
The stochastic resonance model, very similar to the AR model, does not enable the identification of polulation outbreak zones. However, by analyzing the global dynamics of solar activity and local weather conditions, it is possible to predict outbreak development in a region of uniform weather, even if that region is large.
In general, all three models can be applied to solve problems of population monitoring and prediction of forest insect outbreaks. Since the species that cause outbreaks are not limited to those discussed in this article, models of the dynamics of other species can be developed using the techniques suggested in this article.

Author Contributions

Conceptualization, V.S.; methodology, V.S., A.K., O.T. and Y.I.; software, A.K.; validation, A.K. and Y.I.; formal analysis, A.K., V.S. and O.T.; investigation, A.K. and O.T.; data curation, V.S. and A.K.; writing—original draft, V.S.; Writing—review and editing, A.K. and Y.I. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the Russian Science Foundation (grant #23-66-10015).

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Figure A1. Transformed natural data (1) and model calculations (2) of the Siberian silkmoth population dynamics in the Far East.
Figure A1. Transformed natural data (1) and model calculations (2) of the Siberian silkmoth population dynamics in the Far East.
Mathematics 11 04212 g0a1
Figure A2. Spectra of time series of the Siberian silkmoth for the Far East (1) and the Lower Yenisei (3) and for the pine looper in Thuringia (2).
Figure A2. Spectra of time series of the Siberian silkmoth for the Far East (1) and the Lower Yenisei (3) and for the pine looper in Thuringia (2).
Mathematics 11 04212 g0a2
Figure A3. Dynamics of transformed natural data and model calculations of the pine looper population dynamics. 1, transformed natural data; 2, model (Yrs. 1885–1930); 3, forecast (Yrs. 1931–1939).
Figure A3. Dynamics of transformed natural data and model calculations of the pine looper population dynamics. 1, transformed natural data; 2, model (Yrs. 1885–1930); 3, forecast (Yrs. 1931–1939).
Mathematics 11 04212 g0a3

References

  1. Van Mantgem, P.J.; Stephenson, N.L.; Byrne, J.C.; Daniels, L.; Franklin, J.F.; Fulé, P.Z.; Harmon, M.E.; Larson, A.J.; Smith, J.M.; Taylor, A.H. Widespread increase of tree mortality rates in the western United States. Science 2009, 323, 521–524. [Google Scholar] [CrossRef] [PubMed]
  2. Chapin, F.S., III; Matson, P.A.; Mooney, H.A. Principles of Terrestrial Ecosystem Ecology; Springer: New York, NY, USA, 2012; 436p. [Google Scholar]
  3. Peterson, D.L.; Vose, J.M.; Patel-Weynand, T. Climate Change and United States Forests; Advances in Global Change Research; Springer: Dordrecht, The Netherlands, 2014; 57p. [Google Scholar]
  4. Hicke, J.A.; Allen, C.D.; Desai, A.; Dietze, M.; Hall, R.J.; Hogg, E.H.; Kashian, D.; Moore, D.J.P.; Raffa, K.F.; Sturrock, R.N.; et al. Effects of biotic disturbances on forest carbon cycling in the United States and Canada. Glob. Chang. Biol. 2012, 18, 7–34. [Google Scholar] [CrossRef]
  5. Ghimire, B.; Williams, C.A.; James, G.; Collatz, G.J.; Vanderhoof, M.; Rogan, J.; Kulakowski, D.; Masek, J.G. Large carbon release legacy from bark beetle outbreaks across western United States. Glob. Chang. Biol. 2015, 21, 3087–3101. [Google Scholar] [CrossRef] [PubMed]
  6. Hicke, J.A.; Meddens, A.J.H.; Allen, C.D.; Kolden, C.A. Carbon stocks of trees killed by bark beetles and wildfire in the western United States. Environ. Res. Lett. 2013, 8, 035032. [Google Scholar] [CrossRef]
  7. Kurz, W.A.; Dymond, C.C.; Stinson, G.; Rampley, G.J.; Neilson, E.T.; Carroll, A.L.; Safranyik, L. Mountain pine beetle and forest carbon feedback to climate change. Nature 2008, 452, 987–990. [Google Scholar] [CrossRef]
  8. Clark, K.L.; Skowronski, N.; Hom, J. Invasive insects impact forest carbon dynamics. Glob. Chang. Biol. 2010, 16, 88–101. [Google Scholar] [CrossRef]
  9. Albani, M.; Moorcroft, P.R.; Ellison, A.M.; Orwig, D.A.; Foster, D.R. Predicting the impact of hemlock woolly adelgid on carbon dynamics of eastern United States forests. Can. J. For. Res. 2010, 40, 119–133. [Google Scholar] [CrossRef]
  10. Dymond, C.C.; Neilson, E.T.; Stinson, G.; Porter, K.; MacLean, D.A.; Gray, D.R.; Kurz, W.A. Future spruce budworm outbreak may create a carbon source in eastern Canadian forests. Ecosystems 2010, 13, 917–931. [Google Scholar] [CrossRef]
  11. Stinson, G.; Kurz, W.A.; Smyth, C.E.; Neilson, E.T.; Dymond, C.C.; Metsaranta, J.M.; Blain, D. An inventory-based analysis of Canada’s managed forest carbon dynamics, 1990 to 2008. Glob. Chang. Biol. 2011, 17, 2227–2244. [Google Scholar] [CrossRef]
  12. Review of Sanitary and Pathological Condition of Russian Federation Forests for 2021; Russian Forest Protection Center: Moscow, Russia, 2022; 117p. (In Russian)
  13. Ludwig, D.; Jones, D.; Holling, C.S. Qualitative analysis of insect outbreak systems: The spruce budworm and forest. J. Anim. Ecol. 1978, 47, 315–332. [Google Scholar] [CrossRef]
  14. Hastings, A.; Abbott, K.C.; Cuddington, K.; Francis, T.; Gellner, G.; Lai, Y.C.; Morozov, A.; Petrovskii, S.; Scranton, K.; Zeeman, M.L. Transient phenomena in ecology. Science 2018, 361, eaat6412. [Google Scholar] [CrossRef] [PubMed]
  15. Hastings, A.; Abbott, K.C.; Cuddington, K.; Francis, T.B.; Lai, Y.C.; Morozov, A.; Petrovskii, S.; Zeeman, M.L. Effects of stochasticity on the length and behavior of ecological transients. R. Soc. Interface 2021, 18, 20210257. [Google Scholar] [CrossRef] [PubMed]
  16. Abbott, K.C.; Dwyer, G. Using mechanistic models to understand synchrony in forest insect populations: The North American gypsy moth as a case study. Am. Nat. 2008, 172, 613–624. [Google Scholar] [CrossRef] [PubMed]
  17. Ives, A.R.; Abbott, K.C.; Ziebarth, N.L. Analysis of ecological time series with ARMA(p,q) models. Ecology 2010, 91, 858–871. [Google Scholar] [CrossRef]
  18. Sharma, Y.; Abbott, K.C.; Dutta, P.S.; Gupta, A.K. Stochasticity and bistability in insect outbreak dynamics. Theor. Ecol. 2014, 8, 163–174. [Google Scholar] [CrossRef]
  19. Lotka, A. Elements of Physical Biology; The Williams and Wilkim Co.: Baltimore, MD, USA, 1925; 495p. [Google Scholar]
  20. Volterra, V. Leçons sur la Théorie Mathématique de la Lutte pour la Vie; Gauthier-Villars: Paris, France, 1931; 214p. (In French) [Google Scholar]
  21. Murray, J.D. Mathematical Biology; Springer: New York, NY, USA, 2002; 575p. [Google Scholar]
  22. Liebhold, A.M.; Hajek, A.E.; Walter, J.A.; Haynes, K.J.; Elkinton, J.; Muzika, R.M. Historical change in the outbreak dynamics of an invading forest insect. Biol. Invasions 2022, 24, 879–889. [Google Scholar] [CrossRef]
  23. Allstadt, A.J.; Haynes, K.J.; Liebhold, A.M.; Johnson, D.M. Long-term shifts in the cyclicity of outbreaks of a forest-defoliating insect. Oecologia 2013, 172, 141–151. [Google Scholar] [CrossRef]
  24. Williams, D.W.; Fuester, R.W.; Metterhouse, W.W.; Balaam, R.J.; Bullock, R.H.; Chianesei, R.J. Oak defoliation and population density relationships for the gypsy moth (Lepidoptera: Lymantriidae). J. Econ. Entomol. 1991, 84, 1508–1514. [Google Scholar] [CrossRef]
  25. Liebhold, A.M.; Simons, E.E.; Sior, A.; Unger, J.D. Forecasting defoliation caused by the gypsy moth from field measurements. Environ. Entomol. 1993, 22, 26–32. [Google Scholar] [CrossRef]
  26. Isaev, A.S.; Khlebopros, R.G.; Kiselev, V.V.; Kondakov, Y.P.; Nedorezov, L.V.; Soukhovolsky, V.G. Forest Insects Population Dynamics; Publishing House of the Eurasian Entomological Journal: Novosibirsk, Russia, 2009; 115p. [Google Scholar]
  27. Barbosa, P.; Schultz, J. Insect Outbreaks; Academic Press: New York, NY, USA, 1987; 578p. [Google Scholar]
  28. Barbosa, P.; Letourneau, D.; Agrawal, A.A. Insect Outbreaks Revisited; Wiley-Blackwell: Oxford, UK, 2012; 480p. [Google Scholar]
  29. Berryman, A.A. Dynamics of Forest Insect Populations; Plenum Press: New York, NY, USA, 1988; 603p. [Google Scholar]
  30. Berryman, A.A. Population Cycles. Causes and Analysis. Population Cycles: The Case for Trophic Interactions; Oxford University Press: New York, NY, USA, 2002; pp. 3–28. [Google Scholar]
  31. Berryman, A.A. On principles, laws and theory in population ecology. Oikos 2003, 103, 695–701. [Google Scholar] [CrossRef]
  32. Isaev, A.S.; Soukhovolsky, V.G.; Tarasova, O.V.; Palnikova, E.N.; Kovalev, A.V. Forest Insect Population Dynamics, Outbreaks and Global Warming Effects; Wiley and Sons: Hoboken, NJ, USA, 2017; 286p. [Google Scholar]
  33. Prozorov, S.S. The silk moth in fir forests of Siberia. Proc. SibLTI. Krasn. 1952, 93–132. (In Russian) [Google Scholar]
  34. Boldaruyev, V.O. Population Dynamics of the Siberian Silk Moth and Its Parasites; Buryat Publishers: Ulan-Ude, Russia, 1969; 162p. (In Russian) [Google Scholar]
  35. Epova, V.I.; Pleshanov, A.S. Zones of Severity of Phyllophagous Insects in Asian Russia; Nauka: Novosibirsk, Russia, 1995; 147p. (In Russian) [Google Scholar]
  36. Kolomiyets, N.G. Parasites and Predators of the Siberian Silk Moth; Nauka: Novosibirsk, Russia, 1962; 172p. (In Russian) [Google Scholar]
  37. Kondakov, Y.P. Siberian silk moth outbreaks. In Population Ecology of Forest Animals in Siberia; Nauka: Novosibirsk, Russia, 1974; pp. 206–265. (In Russian) [Google Scholar]
  38. Rozhkov, A.S. Siberian Silk Moth; Nauka: Moscow, Russia, 1963; 175p. (In Russian) [Google Scholar]
  39. Rozhkov, A.S. Outbreak of the Siberian Silk Moth and Insect Control Measures; Nauka: Moscow, Russia, 1965; 178p. (In Russian) [Google Scholar]
  40. Soukhovolsky, V.G.; Tarasova, O.V.; Kovalev, A.V. A modeling of critical events in forest insects populations. Rus. J. Gen. Biol. 2020, 81, 374–386. (In Russian) [Google Scholar]
  41. Yurchenko, G.I.; Turova, G.I. Siberian and White-Striped Silkworms in the Far East; DalNIILCH Puplishers: Habarovsk, Russia, 2007; 98p. (In Russian) [Google Scholar]
  42. Barbour, D.A. The pine looper in Britain and Europe. In Dynamics of Forest Insect Populations: Patterns, Causes, Implications; Berryman, A.A., Ed.; Plenum Press: New York, NY, USA; London, UK, 1988; pp. 291–308. [Google Scholar]
  43. Bevan, D.; Brown, R.M. Pine Looper Moth; Forestry Commission Forest Record; HMSO: London, UK, 1978; 119p. [Google Scholar]
  44. Botterweg, P.F. Moth behaviour and dispersal of the pine looper, Bupalus piniarius (L.) (Lepidoptera, Geometridae). Neth. J. Zool. 1978, 28, 341–464. [Google Scholar] [CrossRef]
  45. Broekhuizen, N.; Evans, H.F.; Hassell, M.P. Site characteristics and the population dynamics of the pine looper moth. J. Anim. Ecol. 1993, 62, 511–518. [Google Scholar] [CrossRef]
  46. Broekhuizen, N.; Hassell, M.P.; Evans, H.F. Common mechanisms underlying contrasting dynamics in two populations of the pine looper moth. J. Anim. Ecol. 1994, 63, 245–255. [Google Scholar] [CrossRef]
  47. Demidko, D.A.; Trefilova, O.V.; Kulakov, S.S.; Mikhaylov, P.V. Pine Looper Bupalus piniaria (L.) Outbreaks Reconstruction: A Case Study for Southern Siberia. Insects 2021, 12, 90. [Google Scholar] [CrossRef]
  48. Kendall, B.E.; Ellner, S.P.; McCauley, E.; Wood, S.N.; Briggs, C.J.; Murdoch, W.W.; Turchin, P. Population cycles in the pine looper moth: Dynamical tests of mechanistic hypotheses. Ecol. Monogr. 2005, 75, 259–276. [Google Scholar] [CrossRef]
  49. Klomp, H. The dynamics of field population of the Bupalus piniarius L. Advans. Ecol. Res. 1966, 3, 207–306. [Google Scholar]
  50. Kolomiyets, N.G. Pine looper outbreaks in West Siberia and their practical significance. In Plant Protection; RIO LTA: St. Petersburg, Russia, 1977; Volume 2, pp. 61–63. (In Russian) [Google Scholar]
  51. Palnikova, E.N.; Sviderskaya, I.V.; Soukhovolsky, V.G. Pine Looper in Siberian Forests; Nauka: Novosibirsk, Russia, 2002; 252p. (In Russian) [Google Scholar]
  52. Smits, A.; Larsson, S. Effects of previous defoliation on pine looper larval performance. Agric. For. Entomol. 1999, 1, 19–26. [Google Scholar] [CrossRef]
  53. Smits, A.; Larsson, S.; Hopkins, R. Reduced realised fecundity in the pine looper (Bupalus piniarius) caused by host plant defoliation. Ecol. Entomol. 2001, 26, 417–424. [Google Scholar] [CrossRef]
  54. Smits, A. Performance of pine looper Bupalus piniarius larvae under population build-up conditions. Entomol. Exp. Appl. 2002, 104, 117–124. [Google Scholar] [CrossRef]
  55. Straw, N.A. The impact of pine looper moth, Bupalus piniaria L. (Lepidoptera; Geometridae) on the growth of Scots pine in the Tentsmuir Forest, Scotland. For. Ecol. Manag. 1996, 87, 209–232. [Google Scholar] [CrossRef]
  56. Subklew, W. Untersuchungen uber Bevolkerungsbevegung des Kiefernspanners Bupalus piniarius L. MFF 1939, 10, 10–51. [Google Scholar]
  57. Ginzburg, L.; Colyvan, M. Ecological Orbits: How Planets Moved and Populations Grow; Oxford Univ. Press: New York, NY, USA, 2004; 184p. [Google Scholar]
  58. Dorf, R.C.; Bishop, R.H. Modern Control Systems; Prentice Hall: Hoboken, NJ, USA, 2008; 1018p. [Google Scholar]
  59. Royama, T. Analytical Population Dynamics; Chapman and Hall: London, UK, 1992; 371p. [Google Scholar]
  60. Turchin, P. Complex Population Dynamics: A Theoretical/Empirical Synthesis; Princeton Univ. Press: Princeton, NJ, USA, 2003; 472p. [Google Scholar]
  61. Soukhovolsky, V.G.; Ivanova, Y.D.; Kovalev, A.V. Are Data on Predators Necessary When Modeling Prey Population Dynamics? Biol. Bull. Rev. 2023, 13, 216–227. [Google Scholar] [CrossRef]
  62. Box, G.E.P.; Jenkins, G.M. Time Series Analysis: Forecasting and Control; Holden-Day: San Francisco, CA, USA, 1970; 784p. [Google Scholar]
  63. Anderson, T.W. The Statistical Analysis of Time Series; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 1971; 704p. [Google Scholar]
  64. Bruce, A.D.; Cowley, R.A. Structural Phase Transitions; Taylor & Francis: London, UK, 1981; 326p. [Google Scholar]
  65. Lifshitz, E.M.; Pitaevsky, L.P. Statistical Physics. Part 2. The Theory of Condensed State; Butterworth-Heinemann: Boston, MA, USA, 1980; 387p. [Google Scholar]
  66. Landau, L.D.; Lifshitz, E.M. Statistical Physics; Elsevier: Oxford, UK, 2013; 544p. [Google Scholar]
  67. Anishchenko, V.S.; Astakhov, V.V.; Neiman, A.B.; Vadivasova, T.E.; Schimansky-Geier, L. Nonlinear Dynamics of Chaotic and Stochastic Systens; Springer: Berlin/Heidelberg, Germany, 2002; 446p. [Google Scholar]
  68. Benzi, R.; Sutera, A.; Vulpiani, A. The mechanism of stochastic resonance. J. Phys. A Math. Gen. 1981, 14, 453. [Google Scholar] [CrossRef]
  69. McNamara, B.; Wiesenfeld, K. Theory of stochastic resonance. Phys. Rev. A 1989, 39, 4854. [Google Scholar] [CrossRef]
  70. Lanzara, E.; Mantegna, R.N.; Spagnolo, B.; Zangara, R. Experimental study of a nonlinear system in the presence of noise: The stochastic resonance. Am. J. Phys. 1997, 65, 341–349. [Google Scholar] [CrossRef]
  71. Gammaitoni, L.; Hänggi, P.; Jung, P.; Marchesoni, F. Stochastic resonance. Rev. Mod. Phys. 1998, 70, 223–287. [Google Scholar] [CrossRef]
  72. Mantegna, R.N.; Spagnolo, B.; Trapanese, M. Linear and nonlinear experimental regimes of stochastic resonance. Phys. Rev. E 2000, 63, 011101. [Google Scholar] [CrossRef]
  73. Mei, D.C.; Du, L.C.; Wang, C.J. The effects of time delay on stochastic resonance in a bistable system with correlated noises. J. Stat. Phys. 2009, 137, 625–638. [Google Scholar] [CrossRef]
  74. Benzi, R.; Parisi, G.; Sutera, A.; Vulpiani, A. Stochastic resonance in climatic change. Tellus 1982, 34, 10–16. [Google Scholar] [CrossRef]
  75. Valenti, D.; Fiasconaro, A.; Spagnolo, B. Stochastic resonance and noise delayed extinction in a model of two competing species. Phys. A 2004, 331, 477–486. [Google Scholar] [CrossRef]
  76. Bjørnstad, O.N.; Grenfell, B.T. Noisy clockwork: Time series analysis of population fluctuations in animals. Science 2001, 293, 638–643. [Google Scholar] [CrossRef] [PubMed]
  77. Chichigina, O.A. Noise with memory as a model of lemming cycles. Eur. Phys. J. B 2008, 65, 347–352. [Google Scholar] [CrossRef]
  78. Vilar, J.M.G.; Solé, R.V. Effects of noise in symmetric two-species competition. Phys. Rev. Lett. 1998, 80, 4099–4102. [Google Scholar] [CrossRef]
  79. Chizhevsky, A.L. Physical Factors of the Historical Process. A Short Sketch; 1-st Gostypolitographia: Kaluga, Russia, 1924; 72p. (In Russian) [Google Scholar]
  80. Elton, C.S. Periodic fluctuations in the numbers of animals: Their causes and effects. Brit. J. Exp. Biol. 1924, 11, 119–163. [Google Scholar] [CrossRef]
  81. Huxley, J. Fluctuations of mammals and the sun spot cycle. Harpers Mag. 1927, 12, 42–50. [Google Scholar]
  82. Uvarov, B.P. Insects and Climate; The Entomological Society of London: London, UK, 1931; 247p. [Google Scholar]
  83. Wing, L.W. Migration and solar cycles. Auk 1934, 51, 302–305. [Google Scholar] [CrossRef]
  84. Anthony, H.E. Sun-spots and squirrel emigrations. Lit. Dig. 1934, 20, 30–31. [Google Scholar]
  85. Elton, C.S. Voles, Mice and Lemming: Problems in Population Dynamics; Oxford University Press: London, UK, 1942; 490p. [Google Scholar]
  86. Levengood, W.G. Factors influencing biomagnetic environment during the solar cycle. Nature 1965, 205, 465–467. [Google Scholar] [CrossRef]
  87. Markson, R. Geophysical influences on Biological Cycles. J. Interdiscip. Cycle Res. 1972, 3, 134. [Google Scholar]
  88. Chizhevsky, A.L. Earth Echo of Solar Storms; Mysl.: Moscow, Russia, 1973; 349p. (In Russian) [Google Scholar]
  89. Chernishev, V.B. Solar activity and insects. Probl. Space Biol. 1989, 65, 92–99. (In Russian) [Google Scholar]
  90. Sinclair, A.R.E.; Gosline, J.M.; Holdsworth, G.; Krebs, C.J.; Boutin, S.; Smith, J.N.M.; Boonstra, R.; Dale, M. Can the solar cycle and climate synchronize the snowshoe hare cycle in Canada? Evidence from tree rings and ice cores. Am Nat. 1993, 141, 173–198. [Google Scholar] [CrossRef] [PubMed]
  91. Myers, J.H. Synchrony in outbreaks of forest Lepidoptera: A possible example of the Moran effect. Ecology 1998, 79, 1111–1117. [Google Scholar] [CrossRef]
  92. Puskás, J.; Hill, L.; Kiss, M.; Nowinszky, L. Influence of Solar Activity and Geomagnetism on Moth Species Collected by Pheromone Traps in North Carolina. Biomed. J. Sci. Tech. Res. 2023, 51. [Google Scholar] [CrossRef]
  93. Hamming, R.W. Digital Filters; Duver Publications Inc.: Mineola, NY, USA, 1989; 284p. [Google Scholar]
  94. Wei, W.S. Time Series Analysis: Unvariative and Multivariative Methods; Addison Wesley: Boston, MA, USA, 2006; 614p. [Google Scholar]
  95. Kovalev, A.; Soukhovolsky, V. Analysis of Forest Stand Resistance to Insect Attack According to Remote Sensing Data. Forests 2021, 12, 1188. [Google Scholar] [CrossRef]
Figure 1. The shape of the potential function for the first-order phase transition model. 1, no external field; 2, a weak external field; 3, a strong external field.
Figure 1. The shape of the potential function for the first-order phase transition model. 1, no external field; 2, a weak external field; 3, a strong external field.
Mathematics 11 04212 g001
Figure 2. Partial autocorrelation function (PACF) of insect population dynamics transformed series {L(i)} of population dynamics. 1, silk moth (Lower Enisei); 2, silk moth (Far East); 3, pine looper; 4, Stand. Err.
Figure 2. Partial autocorrelation function (PACF) of insect population dynamics transformed series {L(i)} of population dynamics. 1, silk moth (Lower Enisei); 2, silk moth (Far East); 3, pine looper; 4, Stand. Err.
Mathematics 11 04212 g002
Figure 3. Transformed time series of field data on Siberian silkmoth density in the Lower Yenisei region (1); time series of calculated data in AR-model (2).
Figure 3. Transformed time series of field data on Siberian silkmoth density in the Lower Yenisei region (1); time series of calculated data in AR-model (2).
Mathematics 11 04212 g003
Figure 4. Empirical state function G(ln x) of the insect populations. 1, silk moth (Lowes Enissei); 2, pine looper; 3, silk moth (Far East).
Figure 4. Empirical state function G(ln x) of the insect populations. 1, silk moth (Lowes Enissei); 2, pine looper; 3, silk moth (Far East).
Mathematics 11 04212 g004
Figure 5. The relation between current density and susceptibility of a population of the pine looper. 1, the stably rarefied state of the population; 2, the state 1–2 years before an outbreak; 3, population in the outbreak decline phase; 4, phases of a density maximum and of a population crisis after the outbreak maximum; 5, the value (of the logarithm of density) corresponding to a potential barrier.
Figure 5. The relation between current density and susceptibility of a population of the pine looper. 1, the stably rarefied state of the population; 2, the state 1–2 years before an outbreak; 3, population in the outbreak decline phase; 4, phases of a density maximum and of a population crisis after the outbreak maximum; 5, the value (of the logarithm of density) corresponding to a potential barrier.
Mathematics 11 04212 g005
Figure 6. Cross-correlation functions between solar activity and insect population dynamics. 1, silk moth (Far East); 2, silk moth (Lower Enissei); 3, pine looper; 4, St. Err.
Figure 6. Cross-correlation functions between solar activity and insect population dynamics. 1, silk moth (Far East); 2, silk moth (Lower Enissei); 3, pine looper; 4, St. Err.
Mathematics 11 04212 g006
Figure 7. The time series of first-order differences in air temperature according to weather station Halle (Germany) at +51:30:53, +011:57:02.
Figure 7. The time series of first-order differences in air temperature according to weather station Halle (Germany) at +51:30:53, +011:57:02.
Mathematics 11 04212 g007
Figure 8. Dynamics of the series of Wolf numbers W(t) (1), transformed density L(t) of the pine looper (2), and spectral power S(t) of the series of first-order differences of daily temperature (in degrees K) (3).
Figure 8. Dynamics of the series of Wolf numbers W(t) (1), transformed density L(t) of the pine looper (2), and spectral power S(t) of the series of first-order differences of daily temperature (in degrees K) (3).
Mathematics 11 04212 g008
Figure 9. Dynamics of the series of first-order differences dW of Wolf numbers. 1, dW; 2, years of maximum outbreak snear the Lower Yenisei river.
Figure 9. Dynamics of the series of first-order differences dW of Wolf numbers. 1, dW; 2, years of maximum outbreak snear the Lower Yenisei river.
Mathematics 11 04212 g009
Table 1. Parameters of population dynamics equations of the Siberian silkmoth Dendrolimus sibiricus Tschetv. in the region of the Lower Yenisei, Far East and the pine looper Bupalus piniarius L. in Thuringia.
Table 1. Parameters of population dynamics equations of the Siberian silkmoth Dendrolimus sibiricus Tschetv. in the region of the Lower Yenisei, Far East and the pine looper Bupalus piniarius L. in Thuringia.
Variables in Model (5)Coefficients in Model (5)Std. Err.t-Testp-Value
Silk moth, Lower Yenisei
a00.7070.2582.7450.021
L(i − 2)−0.8690.114−7.6060.000
L(i − 1)1.5860.11214.1380.000
Radj20.96
F138.5
Silk moth, Far East
a0−0.7300.300−2.4320.026
L(i − 2)−0.7820.147−5.3330.000
L(i − 1)1.4260.1589.0220.000
adjR20.830
F46.60
Pine looper, Thuringia
a0−0.2980.170−1.7480.086
L(i − 3)0.5070.1254.0610.000
L(i − 2)−1.4710.196−7.4960.000
L(i − 1)1.8110.12414.5630.000
adjR20.885
F139.900
Table 2. Characteristics of the state function for the analyzed species of phyllophagous insects.
Table 2. Characteristics of the state function for the analyzed species of phyllophagous insects.
Species, Habitatln x1ln x2ln xcln x1c:
Semi-Critical
Density
Type of Potential (Soft/Hard) χ 1 = G ( x c ) G ( x 1 ) ln   x c ln   x 1
Siberian silkmoth, Lower Yenisei−17.632.37.3, hard
Siberian silkmoth, Far East−24.032.54.95, hard
Pine looper, Thuringia−1621.51.67, soft
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Soukhovolsky, V.; Kovalev, A.; Ivanova, Y.; Tarasova, O. Autoregression, First Order Phase Transition, and Stochastic Resonance: A Comparison of Three Models for Forest Insect Outbreaks. Mathematics 2023, 11, 4212. https://doi.org/10.3390/math11194212

AMA Style

Soukhovolsky V, Kovalev A, Ivanova Y, Tarasova O. Autoregression, First Order Phase Transition, and Stochastic Resonance: A Comparison of Three Models for Forest Insect Outbreaks. Mathematics. 2023; 11(19):4212. https://doi.org/10.3390/math11194212

Chicago/Turabian Style

Soukhovolsky, Vladislav, Anton Kovalev, Yulia Ivanova, and Olga Tarasova. 2023. "Autoregression, First Order Phase Transition, and Stochastic Resonance: A Comparison of Three Models for Forest Insect Outbreaks" Mathematics 11, no. 19: 4212. https://doi.org/10.3390/math11194212

APA Style

Soukhovolsky, V., Kovalev, A., Ivanova, Y., & Tarasova, O. (2023). Autoregression, First Order Phase Transition, and Stochastic Resonance: A Comparison of Three Models for Forest Insect Outbreaks. Mathematics, 11(19), 4212. https://doi.org/10.3390/math11194212

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop