The Prediction of Population Dynamics Based on the Spatial Distribution Pattern of Brown Planthopper ( Nilaparvata lugen Stal . ) Using Exponential Smoothing – Local Spatial Statistics Sri

This study aims to predict the population dynamics of Brown Planthopper (BPH) in highly endemic areas of Central Java province, Indonesia. The research was conducted by modifying the method proposed by Legendre and Fortin (1989), through three stages. Those were predicting BPH attacks using Exponential Smoothing Holt Winter, analyzing spatial structure using I, C and Z test on Local Statistic, and making the connectivity inter the periodic predictions of planting season. The results showed that, the studied areas will experience the hotspots phenomenon based on the analysis by the method of Moran's I, Geary's C and Getis Ord Statistic. The analysis of Local Moran's and Getis Ord showed that, four counties namely Boyolali, Klaten, Karanganyar and Sragen experienced a local migration current from region to region around them, whereas other counties are independent. The migration current was influenced by topography, biotic interactions, and anthropogenic factor. Viewed from the spatial scalability in the studied areas, there are four categories of BPH population distribution; point, site, local, and landscape. BPH local migration interregion happened in the County of Klaten, Boyolali, Karanganyar and Sragen. It was caused by some factors: (1) the local climate, (2) the repetition of the use of rice plant variety in a long time, (3) the use of insecticide intensively (3-4 times in one planting period/season), and (4) the irrigation, allowing the spread of BPH larvae and eggs into its surroundings.


Background Research
The population dynamics of Brown Planthopper (Nilaparvata lugens Stal.), hereinafter referred to the BPH, can be modeled and predicted using statistics approach (Damos & Soultani, 2012;Hunter & Price, 1998).BPH population dynamics are influenced by three components, of which the source of food, the natural enemies, and the biotic factors such as weather and topography (Hunter & Price, 1998).Those three components can be a limiting factor of the BPH population development in a particular time and zone.The significant indicator that can be used as an indication of the limiting factors in the BPH population dynamics is the concentration of the population following the spatial and temporal patterns (Perry, 1997).Basically, BPH population dynamics models can be used as a signal occurrence identification of factors that affect population fluctuations and cycles, the reconstruction of population distribution, and the empirical test of the natural population dynamics phenomenon (Ellner, Seifu, & Smith, 2002).Population dynamics can be classified into two categories : descriptive quantitative, is the identification of changes in the total population that can be used to measure the trend and predictions for the future and identification of the various physical and biological factors as the population limiting factors in order to develop a strategy of population control (Juliano, 2007).
Several countries have taken the advantage of population dynamics data for predicting the future population tendency.They have also carried out the strategy of BPH population control through the surveillance program of plant pest, which was called as insect of crop.One of the prediction methods used in pest surveillance program is Holt-Winters Exponential Smoothing, which is later on called as HW (Pal & Gupta, 1994).Based on empirical studies, HW method was chosen for surveillance with some considerations, which are: simpler, adaptive and easy to be developed in the application, more accurate and reliable compared to other methods, and effective for seasonal spatiotemporal data classification on various scalability (Shmueli & Fienberg, 2005;Unkel & Farrington, 2010;Lu, Zeng, & Chen, 2008;Burkom, Muphy, & Shmuelli, 2006).HW method is appropriate for seasonal and trend pattern data, such as BPH surveillance data (Ai, 1999;Suhartono, 2008).
BPH is one kind of a rice plant (Oryza sativa L.) pest, which is the most destructive and has a major impact on food tenacity around the world.In Asia, BPH attacks occurred in Bangladesh, Brunei, Burma (Myanmar), China, Hong Kong, India, Japan, Cambodia, Korea, Laos, Malaysia, Nepal, Pakistan, Philippines, Singapore, Sri Lanka, Taiwan, Thailand, and Vietnam (Catindig et al., 2009;Kuno & Dick, 1984).In 2011, BPH attacks in Indonesia covered 228,133 ha, the highest attack intensity occurred in East Java with an area of 147,857 ha and Central Java with an area of 50,390 ha.In Central Java, the minor attack intensity covered 32 474 ha, 3867 ha for middle attack, 804 ha for severe attack, 13,245 ha for puso (Humas, 2012;Deptan, 2012).
Widely BPH invasion was caused by the migration capability in a remote geographical distance following seasonal wind flow (Monsoon).Physiologically, BPH is capable of flying at the speed of 5-11 m/sec and in an altitude of 1000-3000 AD (Wada et al., 2008;Seino et al., 1987).BPH migration and distribution will lead to changes of the complex spatial and temporal structure comprising topography, climate, anthropogenic and biotic interactions factors (Wang, 2009).Distribution of population caused by topography factor is irrigation, climatic factors (duration of solar radiation, temperature, rainfall and humidity) that support the proliferation of BPH (Win et al., 2011).The population distributions caused by anthropogenic factors are the action of fertilization and using insecticides (Win et al., 2011).Furthermore, distribution of population caused by biotic interactions factor, for example regeneration ability in temperate region can only occur for three generations, whereas when migrating to the tropics, it is able to reach 12 generations per year (Bottrell & Schoenly, 2012).
All this time, the distribution of population modeling has been done using descriptive statistical approach which comprises three methods: K-Ripley, Moran's I, and the autocorrelation function (Aukema et al., 2008).Moran's I method can be used to explain the distribution of a group of individuals in one population of hosts or predators and give an overview of the frequency distribution of the population (Ellner et al., 2001;Ellsbury et al., 1997).This study aims to predict the population dynamics based on the spatial distribution patterns of BPH population in the province of Central Java Indonesia which becomes high endemic of BPH.The research was carried out by modifying the method proposed by Legendre and Fortin (1989) as a contribution in this study.Population dynamics prediction method that uses spatial pattern approach was done in three stages: the prediction of BPH attacks using HW, the analysis of the spatial structure using test I, C and Z on Local Statistics and the connectivity inter the periodicity planting season prediction (Legendre & Fortin, 1989).This paper is organized as follows: Section 1 describes research background.Section 2 describes previous research as a reference of this study.Section 3 describes the theoretical background of the spatial distribution of spatial insect, spatial distribution insect pattern, and local spatial statistics.Section 4 explains the methods and stages of the research.Section 5 describes experiments containing research data, BPH attack prediction, the prediction of spatial distribution, and visualization of the spatial distribution predictions.Section 6 provides the conclusions and the future research.

The Related Works
Insect population dynamics phenomenon was originated from the existence of geographical factor as the limiting evolution process, like the ocean, mountains, and glaciers that were formed 400 million years ago.The acceleration of insect population dynamics had occurred when human motilities have increased drastically in the last 1000 years.It occurs through the formation of lineages in a new location (the invasion) either naturally or anthropogenically.The invasion of insect is divided into three phases, namely: the arrival phase is the process by which individuals move into new areas outside their native range, the growth phase is the process by which a population grows to a sufficient level so that extinction is highly unlikely, and spreading phase is attacking species range expansion into new areas (Liebhold & Tobin, 2008).Insect population dynamics in various scenarios such as changes in weather conditions, changes in the landscape and the abundance of food resources in the environment can be modeled using spatial and temporal approach (Gruebler, Morand, & Daenzer, 2007).
The insect distribution modeling through several approaches comprising are: non-linear method, the estimation of disperse index, the frequency distribution and Geostatistics.In the following decade, Geostatistical methods got more attention from ecologists as the solution in population dynamics (Tobin, Fleischer, & Pitts, 1999).
Geostatistics method use for the development of Integrated Pest Management model using GIS technology (Geographical Information System).This model was applied for several purposes including: modeling the pest and predator insect spatial distribution, detection and monitoring of the pest insect spatial structure in a broader scalability, the determination of high-risk of pest areas based on Agroclimatological information, and compilation of time and intensity of pest emergence prediction.The classification of attack potency was assessed based on the pest spatial and temporal tendencies comprising: the pattern of eating habits, the potency of economic impact, the pattern of population dynamics, the method of spreading and recommendations to handle the attack (Dminić et al., 2010).Minh et al. (2002) developed an early warning model of BPH potential risks in Trungan, Thotnot District, and Cantho City Vietnam by combining several methods including multivariate regression, interpolation and Geostatistics.The study was conducted using data collected from 120 locations and 10 observation periods.The basic principle of this study was to determine the pattern of the relationship of environmental rainfall factors, maximum and minimum air temperature, air humidity and the prevalence of BPH attacks.This model can provide information about predictions of the spatial distribution occurrence prior to the BPH attack.Another observed factor is the high of water surface, the density of natural enemies and periodicity fertilization.Song et al. (1994) conducted the use of geographic information systems study to analyze the spread movement of rice plants (O.sativa) insect in South Korea.The study was conducted using surveillance data from 152 observation stations of pest plant with data from the period of 1981 to 1991.The analysis of BPH spatial distribution patterns was conducted based on the results of the observation value interpolation process throughout the observation stations.This study shows that the BPH population dynamics are influenced by the temperature and migration variable.Prasetyo et al. (2012) had conducted the research of BPH endemics determination using GISA, LISA, and Getis Ord methods in Central Java, Indonesia based on the historical data in 2001-2010.The studied areas consist of 7 counties, constituting the high BPH endemics areas.According to that research, it can be seen that the pattern of spatial object connectivity, such as BPH population centralization/clustering, the rainfall, and the geographical position of the area influences the distribution of BPH population in the whole studied areas.

Insect Spatial Distribution
Spatial distribution is the most dominating characteristic in the life of the insect population.Its attribute is dynamic.As a result, spatial and temporal variability will form different population structure.The knowledge of population structure both spatially and temporally will be able to provide clues to the information, which are: the spread of insect, identification of population dynamics and population density.This information will be useful in pest control framework and understanding of ecological processes that occur on a local scale (Debouzie & Thioulouse, 1986;Pata et al., 2010).The study of the insect population spatial distribution and other biophysical factors may help us to reveal several things related to the life cycle of insect, namely: the characteristics of the environment and the spread of the individual, the development of habitat manipulation strategy especially the beneficial species, and the design of the sampling area determination, assuming that the appropriate design can help to overcome the deviations caused by spatial heterogeneity (Holland et al., 1997).The focus on the determination of the spatial pattern was to see if a population of organisms has a random, homogeneous, or the combination of random and homogeneous distribution.However, this concept is extended to see how great the size of the organisms' population convergence is (Perfecto & Vandermeer, 2008).The spatial pattern is one of the important indicators in identifying organisms that form the dispersion of disease vectors.One type of vector-borne disease will not spread even be endemic if that disease focuses on one location only.Thus, the changes of spatial pattern will affect the changes in spreading and increasingly expanded the organism attacks (Klas, 1965).Generally, the study of spatial patterns for the ecological analysis has four objectives: testing spatial autocorrelation, spatial structure test, the causal variables test, interpolation mapping and the structure of spatial autocorrelation function.Every research goal has different methods and results, as shown in Table 1.The approach of the spatial patterns allows the predictive modeling and detailed mapping to be compiled in order to get a better understanding of the formation of an endemic pattern of a disease.The method of determining spatial patterns of endemicity can be done by the measurement of the studied area (Chadsuthi et al., 2010).The methods of characterization of the spatial patterns based on the research objective divided two categories: the measurement of population distribution pattern by the NNA and QA method, detecting the spatial pattern of organisms attack by the method of SAA.This can be seen in Table 2 (Chaikaew et al., 2009).Neighbor Analysis (NNA) The measurement of spatial pattern for determining population distribution Spatial Autocorrelation Analysis (SAA) Detecting spatial patterns of organism attack based on the attributes of the two scenes.

Local Spatial Statistics
Local spatial statistics was developed to meet the needs of the measurement and analysis of observation data result connectivity in a small area.That information will be useful for: (1) the identification of the population concentration or hotspot, (2) the assessment of the data stationarity structures, and (3) the identification of object distances that is out of reach of the population but has connectivity to population centers (outliers  4. In this equation, is identified area with the georeference = 1, 2,… .Each value is associated by variable value of the research in studied area represented by notation.notation is spatial weight vector, with the value defined as the distance among area.Local Morans'I function according to Du and Chen (2003) is on equation 3.
̅ and ̅ is deviations from the mean value of the observed study variables.notation is a weight matrix element which will be used for determining spatial proximity among areas.
notation is the weighted mean of the deviation around.Local Geary function according to Du and Chen (2003) is on equation 4.
∑ (4) Local Geary was used to calculate the difference between the squared deviations in location i and location j.

Research Procedure
The analysis of time series on surveillance data requires a minimum of 30 periods of observation (Ylioja et al., 1999).The criterion of number of observations area of the research that aims for the exploration of the spatial structure is at least 30 regions (Fortin & Dale, 2005).In accordance with this approach, this study used data of 240 observations and 124 observation area.The research was conducted in four steps, as shown in Figure 1

The Pr
Generally, Gadu plan of July to year of 200 3).Trend i future atta

The Prediction of Spatial Distribution
The prediction of BPH population distribution per county based on the spatial pattern according to the method of Zhang et al. (2009) is in Table 4.  4, the potential occurrence of local migration and population distribution from one area into surrounding areas through a variety of media that have connectivity can be identified.However, the high endemicity has not occured yet because of the humidity factor, as an indicator of low population increase.The low humidity was caused by the low rainfall in the studied area.The analysis of Getis-Ord Statistic is used to detect local spatial association indicated by the formation of hotspots.Hotspots represent the concentration of the number of occurrences of high BPH attacks.On the other hand, concentration of low value, called as coldspots, represents the concentration of a low number of BPH attacks occurrence (Table 5).

The Visualization of Spatial Distribution Prediction
In the analysis of Local Moran's and Getis Ord map, four counties were selected because they showed high spatial connectivity in terms of local migration current from one area to its surrounding based on BPH attack widespread (in hectares).Those counties were Boyolali, Klaten, Karanganyar and Sragen.The other counties showed that BPH attacks did not have spatial connectivity, or in other words, they were independent.The area on the Local Moran's map, which is worth High-High (HH), is called as positive autocorrelation and indicated by the red colour.
The HH area has the high BPH value, surrounded by its nearby areas which have high value as well.It is indicated that there is the local movement and migration of BPH population to nearby areas in HH area.On the other hand, the area on Local Moran's map, indicated by the white colour, has High-Low (HL) value.It means that, the area with high value is surrounded by its neighbors which have low value.This low value area is potential for the occurrence of BPH population distribution from the high value area, called as negative autocorrelation.While the area in yellow on the Local Moran's map is worth insignificant, it means that, it is not significant toward the distribution of BPH population in the nearby areas.spatial connectivity in terms of local migration current from one area to surrounding based on the BPH attack widespread data (in hectares).The counties were Boyolali, Klaten, Karanganyar and Sragen.The other counties did not have spatial connectivity, or they were independent each other.Some factors that support local migration process based on spatial component connectivity are topography, biotic and anthropogenic interactions.From the spatial scalability in the studied area point of view, there were four categories of BPH population distribution: point, site, local, and lanscape.Looking at it from the spatial connectivity concept, BPH local migration interregion happened in the County of Klaten, Boyolali, Karanganyar and Sragen.It was caused by some factors: (1) the local climate (the rainfall, the temperature and the air humidity), (2) the repetition of the use of rice plant variety in a long time, (3) the use of insecticide intensively (3-4 times in one planting period/season), and (4) the irrigation networking, allowing the spread of BPH larvae and eggs into its surroundings.
In our future works, we will conducts the spatial modeling and BPH migration prediction based on the spatial connectivity on regional scalability using spatial statistic approach.The prediction of BPH population migration will be done by using the seasonal BPH data and climate data classification according to local season.
, step I is the classification of BPH outbreak data and rainfall data into three planting season, which are : main planting season, gadu planting season and dry planting season.Step II is the prediction of BPH potential attacks using the exponential smoothing holt-winters method.Step III is the modeling and mapping of the spatial distribution using local spatial statistical approaches, including: Local Moran's (I), Local Geary (C) and the Getis-Ord Statistic (G and G *). Step IV is the visualization of the distribution in the form of: choropleth map, G and G * map, Local Moran's and Geary Local Map and Boxplot.Final step is the analysis and interpretation of results based on the type of spatial patterns which consists of the three forms, namely: the pattern of concentration (Cluster Pattern), the pattern in random (Random Pattern) and the pattern of spread (Dispersed Pattern) as in Table 3(Zhang et al., 2009).

Figure 1 .
Figure 1.The research stages

Table 1 .
Legendre and Fortin (1989)atial pattern for the ecological analysis based on the research objectives according toLegendre and Fortin (1989)

Table 2 .
The characteristic of spatial pattern according toChikaew et al. (2009) (Du & Chen, 200396)atistics function comprising Getis-Ord Statistics called as G and G*, Local Moran's called as I and Local Geary called as C(Getis & Ord, 1996).Local spatial statistics function, that has oftenly been used in practice, comprises G (d) statistics, Local Moran's and Local Geary(Du & Chen, 2003).Getis-Ord Statistics function can be seen in equation 1 and 2, and Local Geary function in Equation

Table 4 .
Zhang et al. (2009)he prediction of six period-Moran's I and Geary's C spatial pattern per county according toZhang et al. (2009)

Table 5 .
Hotspots detection using the analysis of Getis-Ord Statistic