An Autoregulatory Model of Forest Insect Population Dynamics and Forest Stand Damage Dynamics in Different Habitats: An Example of Lymantria dispar L.

: This paper addresses the problem of constructing a mathematical model of population density dynamics and the dynamics of forest areas damaged by spongy moth ( Lymantria dispar L.) outbreaks in the United States, Europe, Russia, and Japan. The key variable of the model is either the pest population density or the area of forests damaged by spongy moths during a season. This variable can be considered proportional to the total current pest abundance in the study area. For the purposes of modeling, data from a number of different authors was used (see bibliography), as well as data from surveys conducted at the egg or caterpillar stage. The complexity of modeling the dynamics of L. dispar abundance is largely due to the fact that, when studying the dynamics of spongy moth population density, the values of external factors such as parasites, predators, and the amount of available food are often unknown. A simple model was proposed using only two types of data: population density and monthly weather characteristics. Our analysis demonstrated that, even in the absence of knowledge regarding the characteristics of ecosystem components interacting with the spongy moth population (parasites, predators, and the state of forage trees), it is possible to introduce models that characterize the regulatory processes in the population in terms of (i) the presence of negative and positive feedbacks in the system and (ii) the influence of external weather factors. The system under investigation was described as an autoregressive


Introduction
One of the main methods for understanding the patterns and mechanisms of forest insect population dynamics is mathematical modeling.From the point of view of a systematic approach, in order to design a population dynamics model of a particular insect species, it is necessary to take into account all the factors of this dynamics: the current population density, the influence of parasites and predators, the state and amount of available food, weather conditions, and so on.It is theoretically possible to write an equation or system of equations that takes into account all these factors.Systems of ordinary differential equations have been used to design such models since the classic works of A. Lotka and V. Volterra [1,2] to their modern interpretations of the late 20th century [3][4][5][6][7].Each of those equations described the dynamics of a separate ecosystem component, i.e., a regulatory factor that depends on the density of the modeled insect population.Those regulatory factors (e.g., parasites, predators, food) were taken into account while constructing the insect population dynamics model.However, such construction of dynamics models raises problems in describing modifying factors, e.g., weather, because it is incorrect to introduce a separate equation into the system for weather (weather does not depend on other regulating factors), and thus we have to introduce the equations of weather influence on model parameters.But this leads to the introduction of additional equations and additional free model parameters, which greatly complicates the verification of the model.
The problem of designing a model is also complicated by the fact that most of the characteristics of insect dynamics factors are not actually measured during field population counts.Most often, only the current population density and meteorological data are determined.Therefore, it is impossible to verify the proposed systems of equations and their related characteristics (population densities of parasites and predators, feed indicators).
Lymantria dispar (Lepidoptera: Erebidae) was used as a model object in this work.This species was chosen for several reasons.Firstly, due to the wide area covering several naturalclimatic zones, it is possible to test the universality of the developed model on populations of one species inhabiting various climatic conditions.Second, the spongy moth is a classical eruptive species with a wide amplitude of population fluctuations [8][9][10][11], which increases the resolution of time series analysis in estimating population density.Third, the spongy moth is a species with a very well-studied ecology, which can make it easier to interpret the resulting information on population dynamics and patterns.In particular, the influence of modifying factors such as weather, regulating factors (producers and second-order consumers), and the interaction of these factors are well described for this species [12].Fourthly, the studied species is an economically important pest for many countries in the northern hemisphere and is included in the list of 100 most dangerous invasive species, so the obtained results on population density modeling will be in demand for applied services (pest control, quarantine service).Thus, using this insect species, we want to validate an ADL model to describe the population dynamics of a species on an area-wide scale.Although the construction of population dynamics models for this species has received some attention in the literature [13], the fundamental novelty of our study is that this model requires a minimal input set (only the number of specimens counted or a quantitative estimate of damage) and that the performance of the model will be tested at the area-wide scale.

Objects of the Study
Lymantria dispar (spongy moth) was chosen as the object of modeling in the present work.The choice is based on this species' wide distribution area.The spongy moth Lymantria dispar L. is one of the most dangerous pests with high invasive potential for deciduous forests.The geographic range of this species includes North America, European countries, China, Japan, and countries of the former Soviet Union and of North Africa [10,[14][15][16].The breadth of the species' area (this species can exist in a wide range of ecological conditions) makes it possible to assess factors of the species' dynamics.The population size of the spongy moth is prone to fluctuations of several orders of magnitude, leading to population outbreaks and the destruction of vast areas of boreal forests.Due to this fact, the species is quite actively studied, and extensive literature on the dynamics of this species is available.A huge number of publications are devoted to the patterns of (and factors affecting) population size dynamics of this species, which are reported to be influenced by a large number of factors, both modifying (primarily weather factors) [17][18][19][20][21][22][23][24] and regulatory (parasites, predators, quality and volume of food, and bacterial and viral diseases) [13,[25][26][27].Nonetheless, in different habitats of the spongy moth, different factors affect population size dynamics, and lengthy field studies are necessary to analyze the population density dynamics of the spongy moth over a long period of time in a given habitat and to identify the influence of these factors.
In the present study, population counts from various sources were used to design and verify the models.The data on continuous counts of spongy moth population density used for modeling is given in Supplement S1.The length of these series was not less than 15 years, which allowed us to obtain model characteristics that are significant at the level of p = 0.05.All presented data is characterized by the presence of both periods when there were population outbreaks, and periods during which the density of populations was low.It should be noted that a set of data on the population dynamics of this species, covering almost all habitat regions around the world, is presented for the first time.Part of the data is the author's material; part was collected from literary sources.
Two types of data were used to analyze the spongy moth's population density dynamics and damage to plants within its geographic range: 1.
Several time series of spongy moth population density in a given habitat.These time series were analyzed for forests in Russia: in the Saratov, Chelyabinsk (two habitats), Moscow, and Khabarovsk regions in different years [9,[28][29][30].

2.
Several time series of area S of seasonal pest damage in a small area.These time series were analyzed for such territories as the states of Massachusetts, New Hampshire, and New Jersey in the USA [13], data on area S of damaged forests in the Tyumen region (Russia) and in Japan (prefectures Hokkaido and Niigata), and data on the average level P of damage to tree crowns on test plots in the northwestern part of France (46.00-49.00• N, 0-1.00 • E).Spongy moth habitat and observation points are shown in Figure 1.Below, we demonstrate the construction of similar models for two different types of initial time series.

Methods
Since it appears that a model of population dynamics for an individual population must be constructed using only field-measured parameters, when there is little or no information on numerous ecosystem components, an analog to the discrete Cobb-Douglas  Below, we demonstrate the construction of similar models for two different types of initial time series.

Methods
Since it appears that a model of population dynamics for an individual population must be constructed using only field-measured parameters, when there is little or no information on numerous ecosystem components, an analog to the discrete Cobb-Douglas equation [31][32][33], widely used in economic modeling problems, can be proposed in order to describe population dynamics: where i is some point in time; X(i − j) is the population density at time (i − j); W(t)dt some accumulated over time period τ value of modifying (e.g., weather) characteristics (or the average value 1 dt of these characteristics over some time period τ), k is an order parameter characterizing the influence of the system's past on its present state.
If there is a time series {X(i)} of accounting data, then, as follows from (1), in order to develop a model, we need to know the order k of Equation ( 1), the values α 1 , . .., α k and the value W(τ).
To account for the influence of the past on the current value of population density, we need to know the "memory" of the system-the length of time k at which, according to Equation (1), the past affects the current state.
From multiplicative Equation (1) we can go to additive form by taking the logarithm of both sides of Equation (1).Then, we obtain: Expression ( 2) is characterized as an ADL (autoregressive distributed lag) equation [34,35] and includes an autoregressive component and a component describing the effect of modifying factors.
In Equation (2), the terms on the right-hand side can be presented as a set of feedbacks that are affected by external and stochastic factors.The ratio of positive and negative feedback determines the temporal dynamics of the modeled system.The conceptual scheme of the model ( 2) is shown in Figure 2. Feedback loops reflect the influence of an object's past states on its current state.For the studied systems, it was assumed that the feedback loops in question are linear and that a feedback signal can be expressed as , with delay τ and feedback amplitude a > 0 for positive feedback and a < 0 for negative feedback.The number of feedback loops of different signs (positive or negative) was unknown, but in the general  Feedback loops reflect the influence of an object's past states on its current state.For the studied systems, it was assumed that the feedback loops in question are linear and that a feedback signal can be expressed as z(i) = a•y(i − τ), with delay τ and feedback amplitude a > 0 for positive feedback and a < 0 for negative feedback.The number of feedback loops of different signs (positive or negative) was unknown, but in the general case, taking into account the feedback loops, a steady-state system can be described as autoregressive (AR).It can be written via the following equation: Quantity a j = ∂ ln X(i) ∂ ln X(i−j) characterizes the susceptibility of the population in year i to a change of population density in "year i − j".In this case, the influence of feedback loops at time t was proportional to variable y(i − j), and the type of feedback (positive or negative) was determined by the sign of a j .If a j is positive, then this parameter denotes linear positive feedback between generations i and i − j.If aj is negative, then aj describes linear negative feedback between these generations.
If the population densities are known, the order k of the autoregressive (AR) component is estimated from the value of the partial autoregressive function [36].
For a known value of k the AR-component can be viewed as a linear regression equation with respect to the unknown parameters a 0 , ..., a k .These parameters can be found using the traditional method of solving linear regression equations.
An autoregressive model to describe population dynamics was proposed earlier [37], but in this model, the lag k = 2 was fixed and the influence of modifying factors was not taken into account.If the modifying factors do not have a strong influence on the population dynamics and the lag k in the system does not differ significantly from 2, such a model seems to be able to describe the observed population dynamics.However, if k ̸ = 2 and the influence of modifying factors are significant, the model may not be correct.
In the approach proposed in this paper, the lag k is not fixed in advance and is estimated for a specific population.The effect of modifying factors is also accounted for.Unfortunately, there is no good theory for choosing a modifying factor and the period during which this factor affects population dynamics.It is possible to propose a random search method, in which some weather indicators and the duration τ of the influence of these indicators on population dynamics are set based on general considerations.Then, a model is calculated for each of the proposed indicators, and an indicator that maximizes the coefficient of determination for Equation ( 2) is picked.
Thus, based on the classical multiplicative Cobb-Douglas model, it is possible to construct a model describing population dynamics under conditions of limited information about the state of the ecosystem.
An important characteristic of a model is its stability.From the point of view of automatic control theory, the stability of any system (including the insect system in the forest) is related to the feedback influence [38].In the presence of strong and fast negative feedback, the system quickly returns to normal after an external influence.In contrast, strong and fast positive feedback leads to strong deviations from the norm.In the presence of both negative and positive feedback, the situation becomes uncertain, and complex indicators that take into account the influence of both negative and positive feedback are needed to assess the stability of the system.
In automatic control theory, such a complex indicator is the stability margin η [39][40][41].The greater the stability margin of a system, the lower the probability of "outliers" in its characteristics (in this case, the rises and falls of insect populations during outbreaks).The stability margin is calculated from the values of the AR-model coefficients characterizing both positive and negative feedback.
During the analysis of these time series, to reduce variance in the time series, data were log-transformed, whereas to reduce noise in a time series, its data were processed using the Hunn filter, which cuts off a high-frequency component (over 0.25 Hz) of a series [42] Forests 2024, 15, 1098 Thus, to build model 2 (see Equation ( 2)), it is necessary to estimate the order of autoregression k and the values and signs of coefficients in Equation ( 2).The order of autoregression in the system and the delay in the set of feedback loops were estimated by means of the PACF.The order of autoregression in the system was estimated using the order of a maximally significant PACF [36].
If the time series {y(i)} and order k of the autoregression of Equation ( 2) are known, then this equation can be considered a linear regression equation with unknown parameters a 0 , . . ., a j , . . ., a k .In this case, feedback coefficients in the AR equation can be found by solving Equation (2) as a regression equation.
Equation ( 2) characterizes feedback loops in the system.Nonetheless, population behavior is affected by modifying factors that do not depend on current and past population density.The impact of the modifying factors can be described by switching from the AR equation to an ADL (autoregressive with a distributed lag) equation [34,43] where W(i) is a characteristic of a modifying factor at time i, and ε is measurement error.For a chosen known modifying factor, susceptibility b j = ∂y(i) ∂ ln W(i−j) of the population to this modifying factor can be determined by treating Equation ( 4) as a regression equation, similarly to the above approach to Equation (2).
Unfortunately, no theory exists that can establish a relationship between variations in insect population densities and weather conditions (temperature and precipitation on specific days and months).To determine which modifying factor at which moment of seasonal dynamics can affect the population in the absence of a developed theory describing the effect of such factors, the following procedure was employed: the weather characteristics from several months were selected as indicators, with the use of these indicators in model (4) resulting in achieving the maximum coefficient of determination.Nonetheless, despite this selection, some weather indicators remain insignificant, even with a p ≤ 0.10.The hydrothermal coefficient (HTC) was found to be the most significant indicator for model (4) of L. dispar.
A hydrothermal coefficient (HTC)-the ratio of precipitation to the average temperature of a month-was selected for calculations.HTC is low in dry and warm weather but increases in cold weather, involving a large amount of precipitation.It is clear that the spongy moth does not follow the calendar, and the period of the impact of weather on populations obviously does not correspond to calendar dates; however, weather variability can be roughly estimated using monthly HTC data.
Next, for this factor, its values were found in certain periods during a season.Then, regression equations of Equation (4) type were analyzed for various periods during a season.HTC for the current year was chosen as the modifying factor with the highest Radj 2 .
For ADL models of population dynamics of forest insects, model stability margin η is an important characteristic that makes it possible to assess a change in the state of a population during various alterations of its environment and of population characteristics themselves [39,40].In this regard, under the conditions of uncertainty about the state of the environment and about the state of the modeled population of forest insects, it is necessary to evaluate the model stability margin.The methodology for calculating the stability margin is given in Supplement S2.

Results
Modeling of insect population dynamics based on time series of population counts.Time series of population density dynamics for spongy moths in the Saratov region of Russia (eggs tree −1 ) in the period 1977-2002 [28] are used for calculating autoregulatory model parameters.We will start the analysis with calculations of the PACF of the smoothed logarithmic time series (Figure 3).
Modeling of insect population dynamics based on time series of population counts.Time series of population density dynamics for spongy moths in the Saratov region of Russia (eggs tree −1 ) in the period 1977-2002 [28] are used for calculating autoregulatory model parameters.We will start the analysis with calculations of the PACF of the smoothed logarithmic time series (Figure 3).As presented in Figure 3, the PACF is a second-order function, i.e., Equation ( 2) should be regarded as an AR(2) equation.
Because the order of autoregression according to Figure 3 can be assumed to be 2.0, the equation describing population dynamics can be written as where HTC(i, j) is a ratio of average monthly precipitation to average temperature in the j-th month of the i-th year.Table 1 presents the calculations of the coefficients for Equation (4) and their statistical significance.As presented in Figure 3, the PACF is a second-order function, i.e., Equation ( 2) should be regarded as an AR(2) equation.
Because the order of autoregression according to Figure 3 can be assumed to be 2.0, the equation describing population dynamics can be written as where HTC(i, j) is a ratio of average monthly precipitation to average temperature in the j-th month of the i-th year.Table 1 presents the calculations of the coefficients for Equation (4) and their statistical significance.As shown in Table 1, the regulation of spongy moth population density in the Saratov region features the presence of two feedback loops: positive feedback between the current population density and population density in the preceding season, and negative feedback between the current density and the density two years before.Positive feedback apparently reflects reproduction in the population: a greater density of parental individuals leads to a higher density of offspring individuals.As for the negative feedback, information about the state of the population is limited, and we know nothing about the effect of parasites or food on the population; therefore, it can be hypothesized that negative feedback that emerges after two generations characterizes (with a delay of 2) the effect of parasites or changes in the state of food resources.Indeed, for the population in the Saratov region in 2006, at the maximum of population density, the reproduction coefficient reached 4.09 while the proportion of parasitized pupae was 5% to 24%; however, this coefficient decreased to 0.304 in 2008, 2 years after the peak of the outbreak, when the proportion of parasitized pupae reached 50%-60% [28].
The susceptibility, b, of the population to changes in HTC should be negative, which means that the lower the HTC (that is, warmer and drier weather), the more rapidly population density increases.If the sign of coefficient b attached to a modifying factor is positive, then this should mean that damper and colder weather is better for the pest.Nevertheless, available literature data contradict this statement, and therefore the model is incorrect, and b should be <0.
The correctness of the model was assessed first by looking at the value of the adjusted coefficient of determination R adj 2 in Equation ( 4): the closer R adj 2 is to 1.0, the greater proportion of variance of variable y(i) is explained by the model.Fisher's F test and t test were performed to evaluate the significance of the model's coefficients [44].
Graphs of time series of log-transformed source data from population monitoring (eggs tree −1 ) and of model data are given in Figure 4.
viduals leads to a higher density of offspring individuals.As for the negative feedback, information about the state of the population is limited, and we know nothing about the effect of parasites or food on the population; therefore, it can be hypothesized that negative feedback that emerges after two generations characterizes (with a delay of 2) the effect of parasites or changes in the state of food resources.Indeed, for the population in the Saratov region in 2006, at the maximum of population density, the reproduction coefficient reached 4.09 while the proportion of parasitized pupae was 5% to 24%; however, this coefficient decreased to 0.304 in 2008, 2 years after the peak of the outbreak, when the proportion of parasitized pupae reached 50-60% [28].
The susceptibility, b, of the population to changes in HTC should be negative, which means that the lower the HTC (that is, warmer and drier weather), the more rapidly population density increases.If the sign of coefficient b attached to a modifying factor is positive, then this should mean that damper and colder weather is better for the pest.Nevertheless, available literature data contradict this statement, and therefore the model is incorrect, and b should be <0.
The correctness of the model was assessed first by looking at the value of the adjusted coefficient of determination Radj 2 in Equation ( 4): the closer Radj 2 is to 1.0, the greater proportion of variance of variable y(i) is explained by the model.Fisher's F test and t test were performed to evaluate the significance of the model's coefficients [44].
Graphs of time series of log-transformed source data from population monitoring (eggs tree −1 ) and of model data are given in Figure 4.

According to Table 1, R adj
2 of the model is 0.90, that is, the proposed model explains 90% of the data variance, and this performance of the model is, of course, quite good.
It should be noted that the proposed model not only describes the amplitudes of a time series of data is but also consistent with data about the phase, as evidenced by the values of the cross-correlation function (CCF) of the data series and model series (Figure 5).According to Table 1, Radj 2 of the model is 0.90, that is, the proposed model explains 90% of the data variance, and this performance of the model is, of course, quite good.
It should be noted that the proposed model not only describes the amplitudes of a time series of data is but also consistent with data about the phase, as evidenced by the values of the cross-correlation function (CCF) of the data series and model series (Figure 5).We can conclude from Figure 5 that the data series and the model series are synchronous, judging by the CCF being close to 1.0 when shift k = 0.
A similar procedure was performed on data about the population density of the  We can conclude from Figure 5 that the data series and the model series are synchronous, judging by the CCF being close to 1.0 when shift k = 0.
A similar procedure was performed on data about the population density of the spongy moth in the Chelyabinsk, Moscow, and Khabarovsk regions.Summary Table 2 characterizes models from Equation (4) for these habitats and indicates that in all these cases, the regulation of population density is mediated by one positive and one negative feedback loop.In this context, the values of the susceptibility of a population's current state to regulation are similar among different habitats.As an illustration of the proposed procedure, Figure 6 shows a time series of logtransformed data and a time series of model data about spongy moth population density dynamics in the Khabarovsk region.As an illustration of the proposed procedure, Figure 6 shows a time series of log-transformed data and a time series of model data about spongy moth population density dynamics in the Khabarovsk region.Modeling based on records of damaged areas.A similar procedure was applied to a time series of the areas of forests damaged by the spongy moth.Obviously, the choice of the area for this analysis is quite arbitrary, and one could use, for example, not an entire US state but rather a part of it (then the area of damage would be smaller).If we suppose that the regulatory processes are identical between the entire study area and its parts, then the model's coefficients (except for the free term of model 4) will not depend on the chosen size of the area for calculations in the model.
In this model of damage area dynamics, positive feedback suggests that the damage to forests in year i may depend on the density of the insect individuals in that year, and this density can be considered proportional to the area S(i − 1) of damage in "year i − 1".Negative feedback in the model of damage area dynamics may reflect restoration of tree foliage in year i after defoliation in year i − 2.
A time series of the area damaged by spongy moths in Massachusetts is presented in Figure 7.In this state, starting in 1980, there has been a linear negative trend in the area of damage, and model 4 describes a steady-state process; therefore, our analysis of the area of dam-aged forests began with the selection of a trend and switching to a detrended series (curve 3 in Figure 7).Modeling based on records of damaged areas.
A similar procedure was applied to a time series of the areas of forests damaged by the spongy moth.Obviously, the choice of the area for this analysis is quite arbitrary, and one could use, for example, not an entire US state but rather a part of it (then the area of damage would be smaller).If we suppose that the regulatory processes are identical between the entire study area and its parts, then the model's coefficients (except for the free term of model 4) will not depend on the chosen size of the area for calculations in the model.
In this model of damage area dynamics, positive feedback suggests that the damage to forests in year i may depend on the density of the insect individuals in that year, and this density can be considered proportional to the area S(i − 1) of damage in "year i − 1." Negative feedback in the model of damage area dynamics may reflect restoration of tree foliage in year i after defoliation in year i − 2.
A time series of the area damaged by spongy moths in Massachusetts is presented in Figure 7.In this state, starting in 1980, there has been a linear negative trend in the area of damage, and model 4 describes a steady-state process; therefore, our analysis of the area of damaged forests began with the selection of a trend and switching to a detrended series (curve 3 in Figure 7).Assuming that the time series DSS(T) can be regarded as autoregressive, let us calculate the PACT of this series (Figure 8).Assuming that the time series DSS(T) can be regarded as autoregressive, let us calculate the PACT of this series (Figure 8).Assuming that the time series DSS(T) can be regarded as autoregressive, let us calculate the PACT of this series (Figure 8).As depicted in Figure 7, order k of the autoregression is 2. Accordingly, we will examine the ADL model where HTC for April of a current year is regarded as a modifying factor (Table 3).When the weather indicator HTC(i, 4) is used in the model, R adj 2 is the highest among all the average monthly weather indicators.Note that the sign of the coefficient that is attached to HTC is correct (−), meaning that the lower the HTC, the greater the area of damage.For the HTC of May or June, the sign of the coefficient attached to this variable is incorrect (+).A time series of data on damage area and the model from Table 4 are illustrated in Figure 9.When the weather indicator HTC(i, 4) is used in the model, Radj is the highest among all the average monthly weather indicators.Note that the sign of the coefficient that is attached to HTC is correct (−), meaning that the lower the HTC, the greater the area of damage.For the HTC of May or June, the sign of the coefficient attached to this variable is incorrect (+).A time series of data on damage area and the model from Table 4 are illustrated in Figure 9.The data series and model series are synchronous, as evidenced by the CCF value (Figure 10).The data series and model series are synchronous, as evidenced by the CCF value (Figure 10).Similarly to the aforementioned analysis of the data about the area of damage in Massachusetts, the models' coefficients were calculated for damage in states New Hampshire and New Jersey (USA), in prefectures Hokkaido and Niigata (Japan), in different regions of France, and in the Tyumen region (Russia); Table 4.  Similarly to the aforementioned analysis of the data about the area of damage in Massachusetts, the models' coefficients were calculated for damage in states New Hampshire and New Jersey (USA), in prefectures Hokkaido and Niigata (Japan), in different regions of France, and in the Tyumen region (Russia); Table 4.
The calculations in Table 4 allow us to compare the properties of spongy moth populations in different habitats.
Order k of autoregression is 2 for all analyzed habitats (both for population densities and for data on the area of damage to plants), with one positive and one negative feedback loop (except for forests in the state of New Jersey, where k = 3).In this context, the development of an outbreak significantly depends on the weather in April or May (Table 5).
The average temperature for all habitats is 10.6 ± 0.53 • C. For some habitats, this temperature is reached in April, and for others, in May.
Are there differences in feedback characteristics between the two time series (data on population density dynamics and data on the area of damaged forests)?Figure 11 shows data from model 4 regarding the spongy moth's population density dynamics and the dynamics of the area of forests damaged by this pest.

Discussion
Calculations on the basis of Tables 1-4 confirm that these conditions are satisfied for all the models in question, and we can say that both population density dynamics models and defoliation area dynamics models can be considered stable.
Methods for analyzing population density dynamics of forest insects and for the construction of ADL models imply that the model type is known and the model's parameters do not change over time.On the other hand, a real population and its environment cannot be modeled with absolute accuracy; they can change in unpredictable ways and may be affected by various perturbing factors.
In such systems, there will always be deviations from an ideal model for various reasons, including:

−
Alterations of the parameters of the population and of its natural environment owing to various circumstances; − The population's and its natural environment's properties that are not taken into account in the model; − The lag-during an interaction of the population with its environment-is not taken into consideration in the model; − A change in the stable state of the population; Figure 11 indicates that in the {a 1 , a 2 } plane, the positions of clusters differ between the population dynamics model and the damage area dynamics model.Nevertheless, according to the t test, there are no significant differences (p ≥ 0.05) between mean values <a 2 > = −0.842± 0.04 from population dynamics models and mean values <a 2 > = −0.838± 0.12 from damage area dynamics models.By contrast, according to the t test, there are significant differences (p < 0.05) between means <a 1 > = 1.507 ± 0.06 from population density dynamics models and means <a 1 > = 1.247 ± 0.09 from damage area dynamics models.

Discussion
Calculations on the basis of Tables 1-4 confirm that these conditions are satisfied for all the models in question, and we can say that both population density dynamics models and defoliation area dynamics models can be considered stable.
Methods for analyzing population density dynamics of forest insects and for the construction of ADL models imply that the model type is known and the model's parameters do not change over time.On the other hand, a real population and its environment cannot be modeled with absolute accuracy; they can change in unpredictable ways and may be affected by various perturbing factors.
In such systems, there will always be deviations from an ideal model for various reasons, including: -Alterations of the parameters of the population and of its natural environment owing to various circumstances;   Figures 12 and 13 show that the stability margin is directly proportional to a1 and inversely proportional to a2.The values of the coefficients of determination of the linear model linking a1 and a2 to η indicate the significance of these relationships.The mean stability margin <η> = 0.08 ± 0.02 for population dynamics models is significantly (p < 0.05) less than <η> = 0.21 ± 0.03 for damage area dynamics models.Such differences may be due to the greater variability of real-world monitoring data used to determine population density dynamics in comparison with the real-world data used to determine the dynamics of the area of damage.Thus, the damage risks calculated from population density data will be somewhat higher than the risks computed from data on the area of damaged forests.
It follows from Figures 12 and 13 that the smaller the absolute values of coefficients a1 and a2, the greater the value of the stability margin η.Since the index η characterizes the stability of the regulation process, in the limit, the population will be stable at large values of η, i.e., at very small values of both coefficient a1, which characterizes the reproduction of the population, and coefficient a2, which describes the lag in the dynamics of the population and is associated with the impact of parasites and predators on it.On the contrary, the population will be unstable at large values of a1 and a2.
However, we modeled the population dynamics of the unstable species L. dispar and found that when a1 > 1.1 and

5
. 0 2 ≥ а , the value of η ≤ 0.10 in half of the studied pop-   Figures 12 and 13 show that the stability margin is directly proportional to a1 and inversely proportional to a2.The values of the coefficients of determination of the linear model linking a1 and a2 to η indicate the significance of these relationships.The mean stability margin <η> = 0.08 ± 0.02 for population dynamics models is significantly (p < 0.05) less than <η> = 0.21 ± 0.03 for damage area dynamics models.Such differences may be due to the greater variability of real-world monitoring data used to determine population density dynamics in comparison with the real-world data used to determine the dynamics of the area of damage.Thus, the damage risks calculated from population density data will be somewhat higher than the risks computed from data on the area of damaged forests.
It follows from Figures 12 and 13 that the smaller the absolute values of coefficients a1 and a2, the greater the value of the stability margin η.Since the index η characterizes the stability of the regulation process, in the limit, the population will be stable at large values of η, i.e., at very small values of both coefficient a1, which characterizes the reproduction of the population, and coefficient a2, which describes the lag in the dynamics of the population and is associated with the impact of parasites and predators on it.On the contrary, the population will be unstable at large values of a1 and a2.
However, we modeled the population dynamics of the unstable species L. dispar and found that when a1 > 1.1 and

5
. 0 2 ≥ а , the value of η ≤ 0.10 in half of the studied population dynamics series.As an outbreak species, population stability should be low.It Figures 12 and 13 show that the stability margin is directly proportional to a 1 and inversely proportional to a 2 .The values of the coefficients of determination of the linear model linking a 1 and a 2 to η indicate the significance of these relationships.The mean stability margin <η> = 0.08 ± 0.02 for population dynamics models is significantly (p < 0.05) less than <η> = 0.21 ± 0.03 for damage area dynamics models.Such differences may be due to the greater variability of real-world monitoring data used to determine population density dynamics in comparison with the real-world data used to determine the dynamics of the area of damage.Thus, the damage risks calculated from population density data will be somewhat higher than the risks computed from data on the area of damaged forests.
It follows from Figures 12 and 13 that the smaller the absolute values of coefficients a 1 and a 2 , the greater the value of the stability margin η.Since the index η characterizes the stability of the regulation process, in the limit, the population will be stable at large values of η, i.e., at very small values of both coefficient a 1 , which characterizes the reproduction of the population, and coefficient a 2 , which describes the lag in the dynamics of the population and is associated with the impact of parasites and predators on it.On the contrary, the population will be unstable at large values of a 1 and a 2 .

Figure 1 .
Figure 1.Spongy moth habitat and observation points (red) on the world map.

Figure 3 .
Figure 3.The PACF of the smoothed logarithmic time series of population density dynamics for the spongy moth in the Saratov region of Russia.1: the PACF; 2: the standard error of the PACF.

Figure 3 .
Figure 3.The PACF of the smoothed logarithmic time series of population density dynamics for the spongy moth in the Saratov region of Russia.1: the PACF; 2: the standard error of the PACF.

Figure 4 .
Figure 4. Time series of log-transformed population monitoring data (curve 1) and of model data (curve 2).

Figure 4 .
Figure 4. Time series of log-transformed population monitoring data (curve 1) and of model data (curve 2).

Figure 5 .
Figure 5.The cross-correlation function (CCF) of the data series and model series.1: the CCF; 2: the standard error of the CCF.

Figure 5 .
Figure 5.The cross-correlation function (CCF) of the data series and model series.1: the CCF; 2: the standard error of the CCF.

Figure 6 .
Figure 6.A series of log-transformed data (1) and a model series (2) of spongy moth population density dynamics in Khabarovsk region.

Figure 6 .
Figure 6.A series of log-transformed data (1) and a model series (2) of spongy moth population density dynamics in Khabarovsk region.

Figure 9 .
Figure 9. Series DSS(T) (curve 1) and the model series (curve 2) for the area affected by spongy moth damage in Massachusetts.

Figure 9 .
Figure 9. Series DSS(T) (curve 1) and the model series (curve 2) for the area affected by spongy moth damage in Massachusetts.

Figure 10 .
Figure 10.The CCF of series DSS(T) (curve 1) and of the model series (curve 2).

Figure 11 .
Figure 11.Coefficients a1 and a2 for models of population density dynamics (1) and models of damage area dynamics (2).

Figure 11
Figure 11 indicates that in the {a1, a2} plane, the positions of clusters differ between the population dynamics model and the damage area dynamics model.Nevertheless, according to the t test, there are no significant differences (p ≥ 0.05) between mean values <a2> = 04 .0 842 .0 ± − from population dynamics models and mean values <a2> =

Figure 11 .
Figure 11.Coefficients a 1 and a 2 for models of population density dynamics (1) and models of damage area dynamics (2).

Figures 12 and 13
Figures 12 and 13 depict the relation of coefficients a 1 and a 2 from model 3 with the stability margin for population density dynamics models and for damage area dynamics models.

ForestsFigure 12 .
Figure 12.The relation between coefficient a1 and stability margin η for population density dynamics models (1) and damage area dynamics models (2).

Figure 13 .
Figure 13.The relation between coefficient a2 and stability margin η for population density dynamics models (1) and damage area dynamics models (2).

Figure 12 .
Figure 12.The relation between coefficient a 1 and stability margin η for population density dynamics models (1) and damage area dynamics models (2).

ForestsFigure 12 .
Figure 12.The relation between coefficient a1 and stability margin η for population density dynamics models (1) and damage area dynamics models (2).

Figure 13 .
Figure 13.The relation between coefficient a2 and stability margin η for population density dynamics models (1) and damage area dynamics models (2).

Figure 13 .
Figure 13.The relation between coefficient a 2 and stability margin η for population density dynamics models (1) and damage area dynamics models (2).

Table 1 .
Computed coefficients of Equation (4) and their statistical significance for the spongy moth population in the Saratov region (Russia).The Hydrothermal coefficient HTC(i, 4) is the hydrothermal coefficient for April of the current year.

Table 2 .
Characteristics of models from Equation (4) for different habitats of the spongy moth in Russia.

Table 4 .
ADL(2,1): the model of dynamics of the area affected by spongy moth outbreaks in different habitats.

Table 4 .
ADL(2,1): the model of dynamics of the area affected by spongy moth outbreaks in different habitats.

Table 5 .
Average air temperatures in a month included in the ADL model of dynamics.