Carabid beetles in arable land move between different habitats to exploit resources that vary in time and space. Understanding such movement is key to explaining how the pest control function of carabids in arable crop lands depends on the spatial configuration of crop fields and associated semi-natural habitats, but movement at and beyond field scale is not well understood. Here we use a model selection framework to identify and parameterize a parsimonious movement model, based on mark-release-recapture data in two adjacent arable crop fields, one planted with rye, and the other with oil radish. The simplest model assumes motility of beetles to be the same in the two crops, and it does not consider losses of beetles over time due to death or mark loss. These assumptions are relaxed either separately or together in competing models, resulting in a comparison between four models. All models consider the effect of spatially heterogeneous pitfall trapping on the size of the moving population. Models were fitted to data with Poisson likelihood, and Akaike’s information criterion was then used to rank the models. The model selection showed that including a parameter for loss of beetles due to mortality or mark loss resulted in the best approximation of the observed data. The data did not support the assumption of different motility between the two crops. We conclude that our extended model can be used to simulate beetle recapture in mark-release-recapture experiments, but further refinements to the model are needed. The inverse modeling framework for model identification and parameter estimation that was applied in this study proved effective to select the most promising model and parameter values

Allema, A.B.

Berghuijs, H.N.C.

Groot, J.C.J.

Hemerik, L.

Rossing, W.A.H.

Werf, W., van der

Name not available

Fokker-Planck model for movement of the carabid beetle Pterostichus melanarius in arable land: Model selection and parameterization

van der Werf, W.

NARCIS 

Wageningen Yield 

Herman N.C. Berghuijs1,2, Bas Allema2, Lia Hemerik1, Wopke van der
Werf3, Jeroen C.J. Groot2 & Walter A.H. Rossing2
1Biometris, Department of Mathematical and Statistical Methods, Wageningen
University, PO Box 100, 6700 AC Wageningen, The Netherlands, E-mail:
lia.hemerik@wur.nl; 2Farming Systems Ecology Group, Wageningen University, PO
Box 563, 6708 PB Wageningen, The Netherlands; 3Centre for Crop Systems Analysis,
Wageningen University, PO Box 430, 6700 AK Wageningen, The Netherlands
Carabid beetles in arable land move between different habitats to
exploit resources that vary in time and space. Understanding such
movement is key to explaining how the pest control function of cara-
bids in arable crop lands depends on the spatial configuration of crop
fields and associated semi-natural habitats, but movement at and
beyond field scale is not well understood. Here we use a model selec-
tion framework to identify and parameterize a parsimonious move-
ment model, based on mark-release-recapture data in two adjacent
arable crop fields, one planted with rye, and the other with oil radish.
The simplest model assumes motility of beetles to be the same in the
two crops, and it does not consider losses of beetles over time due to
death or mark loss. These assumptions are relaxed either separately or
together in competing models, resulting in a comparison between four
models. All models consider the effect of spatially heterogeneous pit-
fall trapping on the size of the moving population. Models were fitted
to data with Poisson likelihood, and Akaike’s information criterion
was then used to rank the models. The model selection showed that
including a parameter for loss of beetles due to mortality or mark loss
resulted in the best approximation of the observed data. The data did
not support the assumption of different motility between the two
crops. We conclude that our extended model can be used to simulate
beetle recapture in mark-release-recapture experiments, but further
refinements to the model are needed. The inverse modeling framework
for model identification and parameter estimation that was applied in
this study proved effective to select the most promising model and
parameter values.
Keywords: Pterostichus melanarius, dispersal, Poisson distribution, AIC,
model selection
Fokker-Planck model for movement of the cara -
bid beetle Pterostichus melanarius in arable
land: 
Model selection and parameterization
PROC. NETH. ENTOMOL. SOC. MEET. - VOLUME 23 - 2012 21
22
In mark-release-recapture (MRR) experiments, a large number of insects are
released into an area and recaptured over time. Data from this type of experi-
ments either consist of dispersal distances or the locations and time at which a
certain number of beetles is recaptured (Turchin 1998). Analysis of these data
can provide insight in insect population dispersal patterns. A common way to
analyse such MRR data is to consider animal population spread as a continuous
diffusion process (Inoue 1978, Kareiva 1982, Turchin 1998, Cronin et al. 2000,
Yamamura 2002, Reeve et al. 2008). Diffusion in one or more dimensions is des-
cribed by partial differential equations (PDEs). According to Yamamura et al.
(2002), one advantage of PDEs is their mathematical tractability, i.e. there are
analytical solutions available that describe animal density in time and space.
Unfortunately, analytical solutions of these models cannot take spatially hetero-
geneous removal by recapture into account. If a substantial part of the beetle
population is removed in the course of the experiment, the process of beetle
removal by catch must be accounted for. Fagan (1997) presented an experimen-
tal approach to do this. He conducted a MRR experiment in which he released
mantids and recaptured them at the boundaries of the experimental field. He
used a diffusion model to analyse his data assuming absorbing boundaries. This
approach, however, cannot be applied to analyse MRR data that consist of counts
of insects that are recaptured within the experimental field. Yamamura et al.
(2003) solved this problem by adding a loss term to the diffusion equation for
beetle removal due to catching. Their model, however, requires that traps are
placed uniformly in a lattice pattern. 
Usually, PDE-based movement models for insects are based on Fickian dif-
fusion. However, as pointed out by Turchin (1998), ecological diffusion should
be based on the Fokker-Planck equation. The distinction is important because
the Fokker-Planck equation results in a non-homogeneous steady state, where
local density is inversely proportional to motility, whereas ordinary Fickian dif-
fusion demonstrates continuous dissipating density gradients until the density
has become spatially homogeneous.
Here, we present a movement model that accounts for local removal due to
catching that can be applied for any spatial configuration of traps. Dispersal in
this model is based on the Fokker-Planck equation. We identified the model
using data from an MRR experiment with the common carabid beetle
Pterostichus melanarius Illiger in arable land. We wanted to know if, besides
removal due to trapping, also a removal term should be included in the model for
mortality and mark wear, and secondly, we wanted to know if a different motil-
ity should be implemented for two fields grown with different crop species. We
formulated four Fokker-Planck models to answer these questions: (1) without
mortality and without a crop specific motility; (2) with mortality and without a
crop specific motility; (3) without mortality and with a crop specific motility;
and (4) with mortality and with a crop specific motility. In all these models
MODEL FOR MOVEMENT OF A CARABID BEETLE
mortality includes mark wear. Model selection was used to rank models and
select those with the greatest support from the data (Hilborn & Mangel 1997,
Bolker 2008). Ranking was based on Akaike’s information criterion (Akaike
1974, van der Hoeven et al. 2005).
MATERIAL AND METHODS
Field experiment
A mark-release-recapture experiment with P. melanarius was performed on par-
cel 5 of the organic university farm Droevendaal (Wageningen, The
Netherlands) in 2009. The parcel measured 229 × 52 m2 and was bordered at all
sides by a 3-6 m wide grass margin and at the north side by a grass margin with
a hedgerow of 1.5-2.5 m tall shrubs. The MRR experiment was carried out in a 120
× 52 m study area within this parcel. One half of the area was grown with win-
ter rye (Secale cereale var. Admiraal) and the other half with oil radish (Raphanus
sativus var. Brutus). Both crops were sown in the first week of August 2009 to a
row spacing of 12.5 cm.
The beetles were released on 7 September 2009 at 6 p.m. (1015 in oil radish and
1015 in winter rye) along a 12-m-long line in both crops (Figure 1). Beetles
released in winter rye were marked with gold nail polish (OPI Nail lacquer, NL
H41) while beetles released in oil radish were marked with pink nail polish (OPI
Nail lacquer, NL B777). Trapping stations (Ø 8.5 cm) were placed in parallel
lines at 10, 20 and 30 m from the release lines. The released beetles had been col-
lected by pitfall trapping in a triticale field on the same farm during summer.
Beetles were sexed and only females were retained. They were stored in 300 ×
400 mm containers in a dark room at 4 °C until the time of release. Beetles were
HERMAN N.C. BERGHUIJS, BAS ALLEMA, LIA HEMERIK ET AL.
23
Figure 1. The set-up of the mark-release-recapture experiment with oil radish (left) and
winter rye (right). Dots indicate the trap positions. Dashed lines indicate the beetle
release lines.
fed with fly maggots (Lucillia caesar) during storage. The mark was applied a few
days before release. A plastic screen of 1.40 m was placed behind (as seen from
the release site) the traps at 20 and 30 m from the release line to increase trap-
ping probability. The trapping stations were emptied to count recaptured beetles
between September 8 and September 30, 2009, with an interval between catches
of 1, 2 or 3 days. Recaptured beetles were not returned to the field.
Model development
We compare four movement models. We simulated beetle dispersal according to
the Fokker-Planck equation and included a term for the fraction of beetles that
was trapped per unit of time because a substantial part of the marked beetles was
captured (in total 468 pink and 528 gold beetles). A loss rate was also included
because marked beetles could die, or more likely, loose their mark over time. The
PDE of the most complicated model (our model 4, see below) is given in equa-
tion 1. The other three PDEs are simplified versions of equation 1. 
(1)
for x ≤ 60 it holds that μ(x,y) = μ1 and for x > 60 it holds that μ(x,y) = μ2.
Here, N(x,y,t) is the number of beetles at location (x,y) at time t, μ1 is the
motility (m2 day−1) in oil radish, μ2 the motility in winter rye, ξ (day−1) is the loss
rate due to mark wear and mortality and α(x,y) (day−1) is the loss rate due to
recapture at location (x, y). Our simplest model, model 1 describes beetle disper-
sal with local recapture: μ1 = μ2 and ξ = 0. In model 2 the basic model is extend-
ed with spatial independent loss rate: μ1 = μ2 and ξ ≠ 0. In model 3 the basic model
is extended with a habitat specific motility: μ1 ≠ μ2 and ξ = 0. Finally, in model 4
the basic model is expanded with both the loss rate and a habitat specific motil-
ity: μ1 ≠ μ2 and ξ ≠ 0. For the catch we estimated the global catch efficiency ω
with α(x,y) ≡ φω μ(x,y)/ΔxΔy, with φ = 1 at the location of the trap and φ = 0 at
all other locations.
In order to simulate beetle dispersal in time and space, boundary conditions
and initial conditions have to be formulated. In the field experiment beetles were
released on a 12 m transect in winter rye and oil radish. In the model it is
assumed that the beetles were released from four equally spaced points on these
release lines.
Several boundary conditions were compared. Reflecting boundaries are not
appropriate because they do not allow beetles to leave the field, which they can
do in reality. On the other hand, zero boundary conditions are neither appropri-
ate, because they would allow beetles to leave the field, never to return, result-
ing in an overestimation in beetle losses from the field. We explored model
MODEL FOR MOVEMENT OF A CARABID BEETLE
24
][ )t,y,x(N)y,x(
y
)t,y,x(N)y,x(
x
)t,y,x(N)y,x(
t
)t,y,x(N ???? ??
?
??
?
??
?
?
2
2
2
2
behaviour for different boundary conditions and found that a ‘slow release
boundary’ around the field gave the most satisfactory and credible simulation
results. The slow-release boundary consists of a 1-m-wide strip around the field
with a low motility (μ3). This low motility results in an aggregation of beetles in
this margin, and a slow release back to the field from this boundary. The value
of μ3 was determined once by optimisation for the least restricted model (Eq. 1)
and then used as a constant for the optimization of the other models.
Optimisation was based on the Poisson negative log-likelihood (see below) and
resulted in a value of μ3 of 3 m2 day−1. 
Numerical approximation
None of the four considered models can be solved analytically, because φ is spa-
tially heterogeneous. Consequently, numerical approximations are required in
order to run simulations and analyse the field data. The PDEs were solved
numerically by the forward central difference method.
Estimation of the parameters
The most suitable model and the value of model parameters were identified
using an inverse modelling approach employing the evolutionary algorithm of
Differential Evolution (Storn & Price 1997). The Poisson negative log likelihood
(NLL) was used as a measure for goodness of fit. This measure was calculated
as follows:
(2)
Here, Ndata(x,y,t) represents the number of trapped beetles in the field exper-
iment at a given time and location, Npredicted(x,y,t) represents the number of
trapped beetles at a certain time predicted by the model corresponding to the
chosen model (with or without habitat specific movement).
Model selection
Based on the number of estimated parameters p and the NLL, models can be
compared by using Akaike’s information criterion (Akaike 1974, van der Hoeven
et al. 2005, Hilborn & Mangel 1997):
AICi = 2NLLi + 2pi (3)
Here, i is the index of the model (i = 1…4). The AIC considers all models with
ΔAIC < 2 equivalent. For model validity also AIC-weights can be calculated
within the set of four models considered. Here the AIC weight for model i (with
i = 1…4) is:
(4)
HERMAN N.C. BERGHUIJS, BAS ALLEMA, LIA HEMERIK ET AL.
25
? ???? ???
x y t
t,y,xNt,y,xNt,y,xNt,y,xN ])!(ln[)](ln[)()(NLL datapredicteddatapredicted
?
?
??
??
? 4
1i
i
i
i
)AICexp(
)AICexp(
weightAIC
RESULTS
Comparison of models with Akaike’s information criterion demonstrated that
models with beetle losses (models 2 and 4) were far superior over those without
losses, the former having AIC values that were approximately 238 smaller than
the latter (Table 1). The better fit of models with beetle losses is graphically
illustrated in Figure 2. According to the criterion that models differ if their AIC-
values differ more than 2 we can conclude that no difference exists between
models 2 and 4. This is confirmed by the AIC weights of 0.66 for model 2 and
0.34 for model 4, indicating support for both models. Nevertheless, the AIC
weights of models 2 and 4 show that model 2 which has a single parameter for
the motility is more likely to explain the observations in the field experiment
than model 4 which has crop-specific parameters for motility.
DISCUSSION AND CONCLUSIONS 
The results convincingly show that loss of beetles should be accounted for in the
model describing dispersal. We have also shown that the P. melanarius popula-
tion dispersal in the field experiment is more likely to be similar across crop
habitats than habitat specific. 
The model could be made more specific by accounting for the plastic screens
that were placed behind the traps at 20 and 30 m (but not at 10 m) to increase
catching rate. This can be achieved by assigning different parameter values for
trap locations with or without a screen. The model could also be expanded with
the addition of the grass margin at the northern side and the grass margin with-
out hedgerow at the southern side of the simulated field. The disadvantage of
this addition is that there are no recapture data from these habitats and that the
motility in that border can only be estimated indirectly from data collected in oil
radish and winter rye.
Although the speed of beetles’ dispersal seems not to be affected by different
habitat types, it does not necessarily mean that beetles do not have a preference
for one of the crop types.  In further research, it could be investigated whether
beetle dispersal is influenced by the interface between two habitat types, rather
than by the habitat type itself. This approach may provide more insight in
whether beetles are attracted or repelled by a habitat type if they approach the
interface between different habitat types.
One disadvantage of MRR experiments is that recapturing beetles dilutes the
population density in the field by preventing further dispersal of beetles once
they are recaptured. In the ‘boundary flux approach’ designed by Fagan (1997),
reaction-diffusion models are used to analyse beetle dispersal by considering
recapturing at the edge of the experimental field. This approach makes it no
longer necessary to account for beetle population density dilution by recapturing
beetles in the field. However, this type of experiment does not allow using data
from traps within the experimental field and indirectly estimates beetle densi-
MODEL FOR MOVEMENT OF A CARABID BEETLE
26
HERMAN N.C. BERGHUIJS, BAS ALLEMA, LIA HEMERIK ET AL.
27
Table 1. Results of the model selection for the four Fokker-Planck models describing
beetle dispersal in two crops, arranged by increasing AIC; Model: model number as
explained in text; p: number of estimated parameters in the model. Model parameters:
ξ: loss rate; μ1 and μ2: motility in oil radish and winter rye crops (if μ1 = μ2 then the
value is put in the middle); ω: dimensionless global catch efficiency of the traps; AIC:
Akaike’s information criterion; ΔAIC: difference in AIC with the best model; AIC
weight: a measure of belief in the current model (out of the considered four models).
Model p ξ (day-1) μ1 (m2 day-1) μ2 (m2 day-1) ω (-) AIC ΔAIC AIC weight
2 3 0.067 153 0.15 2009.6 0.0 0.66
4 4 0.067 138 174 0.15 2010.2 0.7 0.34
1 2 224 0.09 2248.1 238.5 0.00
3 3 198 255 0.09 2248.9 239.3 0.00
Figure 2. The cumulative number of beetles recaptured throughout the experiment
(open circles) and the simulated number of recaptured beetles (solid lines) as estimated
by 4 considered models with different diffusion constants (μ1 and μ2 for the two crops)
and loss rate (ξ). (a) Model 1: μ1 = μ2 and ξ = 0; (b) model 2: μ1 = μ2 and ξ ≠ 0; (c)
model 3: μ1 ≠ μ2 and ξ = 0; (d) model 4: μ1 ≠ μ2 and ξ ≠ 0.
28
ties over time within the experimental field. Yamamura et al. (2003) solved this
problem by adding a loss term in the diffusion equation for removal due to trap-
ping. Their solution, however, requires that traps are placed uniformly in a lat-
tice pattern. As far as we know, our study is the first one to identify a dynamic
movement model based on the Fokker-Planck equation from MRR data collect-
ed with an irregular grid of traps. 
REFERENCES 
Akaike, H. 1974. A new look at the statistical model identification. IEEE Transactions on
Automatic Control 19: 716-723.
Bolker, B.M. 2008. Ecological Models and Data in R. Princeton University Press, 396 pp. 
Cronin, J.T., Reeve, J.D., Wilkens, R. & Turchin, P. 2000. The pattern and range of
movement of a checkered beetle predator relative to its bark beetle prey. Oikos 90:
127-138.
Fagan, W. 1997. Introducing a ‘boundary flux’ approach to quantifying insect diffusion
rates. Ecology 78: 579-587.
Hilborn, R. & Mangel, M. 1997. The Ecological Detective: Confronting Models with
Data. Princeton University Press, 330 pp.
van der Hoeven, N., Hemerik, L. & Jansen, P.A. 2005. Balancing statistics and ecology:
on the lumping of experimental data for model selection. In: T.A.C. Reydon & L.
Hemerik (eds.) Current Themes in Theoretical Biology: a Dutch Perspective.
Springer, Dordrecht, The Netherlands, pp. 233-263.
Inoue, T. 1978. A new regression method for analyzing animal movement patterns.
Resource Population Ecology 20: 141-163.
Kareiva, P. 1982. Experimental and mathematical analysis of herbivore movement: quan-
tifying the influence of plant spacing and quality on foraging discrimination.
Ecological Monographs 52: 261-282.
Reeve, J.D., Cronin, J.T. & Haynes, K.J. 2008. Diffusion models for animals in complex
landscapes: incorporating heterogeneity among substrates, individuals, and edge
behaviours. Journal of Animal Ecology 77: 898-904.
Storn, R. & Price, K., 1997. Differential evolution - A simple and efficient heuristic for
global optimization over continuous spaces. Journal of Global Optimization 11: 341-
359.
Turchin, P. 1998. Quantitative Analysis of Movement - Measuring and Modeling
Population Redistribution in Animals and Plants. Sinauer Associates, Sunderland,
MA, USA.
Yamamura, K. 2002. Dispersal distance of heterogeneous populations. Population Ecology
44: 93-101.
Yamamura, K., Kishita, M., Arakaki, N., Kawamura, F. & Sadoyama, Y. 2003. Esti mation
of dispersal distance by mark-recapture experiments using traps: correction of bias
caused by the artificial removal by traps. Population Ecology 45: 149-155.
MODEL FOR MOVEMENT OF A CARABID BEETLE


Fokker-Planck model for movement of the carabid beetle Pterostichus melanarius in arable land: Model selection and parameterization

Abstract

Similar works

Full text

Available Versions

Name not available

NARCIS

Wageningen Yield