OR I G I N A L AR T I C L E
Local population dynamics of gray wolves Canis lupus and
Eurasian lynx Lynx lynx exhibit consistency with
intraspecific contest competition models
Jurģis Sˇuba1 | Yukichika Kawata2 | Andreas Lindén3
1Latvian State Forest Research Institute
“Silava”, Salaspils, Latvia
2Faculty of Economics, Kindai University,
Higashiosaka, Osaka, Japan
3Natural Resources Institute Finland,
Helsinki, Finland
Correspondence
Jurģis Sˇuba, Latvian State Forest Research
Institute “Silava”, Rigas Str. 111, Salaspils,
LV 2169, Latvia.
Email: jurgis.suba@silava.lv
Funding information
European Regional Development Fund,
Grant/Award Number: 1.1.1.2/
VIAA/3/19/511
Abstract
In Europe, the gray wolf and Eurasian lynx populations are recovering after vari-
ous levels of persecution. The two species differ in their social structure and spa-
tial patterns of aggregation. Using model selection, we investigated the
consistency of the available time series data on local wolf and lynx sub-
populations with a number of single-species population growth models that per-
tain to two types of intraspecific competition, namely, scramble (SC) and contest
competition (CC), and reflect random (R) or aggregated (A) distribution of indi-
viduals. The applied models of population growth—the Ricker (SCR), Skellam
(CCR), Hassell (SCA), and Beverton–Holt (CCA) models—were all parameterized
in terms of intrinsic growth rate and carrying capacity with unified definitions.
The projected carrying capacity was allowed to show a temporal trend, which was
justified by an observed increase in prey abundance in recent decades. For both
species, the models pertaining to contest competition outperformed the scramble
competition models, and the Beverton–Holt model had the greatest weight. How-
ever, for the lynx, the difference of performance between the scramble and contest
competition models was considerably smaller than that for the wolves. In most of
the models, when it was meaningful, an optional time lag operator was added to
account for a delay in individual maturity and reproduction. However, the models
with a time lag had a worse fit than the models without it. This study promotes
the application of population models that reflect intraspecific competition for
modeling population dynamics in a single- or multi-species framework.
KEYWORD S
Canis lupus, contest competition, Lynx lynx, population dynamics, scramble competition
1 | INTRODUCTION
It is common practice to adopt single-species models
when investigating certain aspects in population dynam-
ics, such as time lags (Ruan, 2006), density dependence
(Boukal & Berec, 2002), and environmental conditions
(Liu & Wang, 2009). The important intraspecific factors
affecting population growth also include aspects of the
nature of competition (Brännström & Sumpter, 2005). As
indicated by simulation experiments in previous studies
Received: 18 January 2023 Revised: 23 July 2024 Accepted: 2 August 2024
DOI: 10.1002/1438-390X.12197
This is an open access article under the terms of the Creative Commons Attribution-NonCommercial License, which permits use, distribution and reproduction in any
medium, provided the original work is properly cited and is not used for commercial purposes.
© 2024 The Author(s). Population Ecology published by John Wiley & Sons Australia, Ltd on behalf of The Ecological Society of Japan.
Population Ecology. 2024;1–18. wileyonlinelibrary.com/journal/pope 1
(e.g., Johst et al., 2008), population dynamics may reflect
intraspecific competition patterns and resource partition-
ing. By recognizing these factors and selecting an appro-
priate population model, one may improve the projection
of population dynamics for more appropriate manage-
ment and effective conservation.
After previous persecution and deliberate eradication,
populations of large carnivores have been recovering and
locally re-emerging due to internationally binding legisla-
tion favoring conservation and reintroduction efforts
(Boitani et al., 2015; Chapron et al., 2014). As the avail-
ability of remote wilderness areas is limited, the future of
large carnivores depends on sustainable human–wildlife
coexistence in the human-dominated landscape
(Chapron et al., 2014; Everatt et al., 2019; Lennox
et al., 2018; Linnell et al., 2001; L
opez-Bao et al., 2017).
Under such circumstances, viable population manage-
ment should be based on understanding the characteris-
tics of carnivore population dynamics and the functional
form of population growth.
Two extremes of intraspecific competition, namely,
scramble and contest competition, are generally distin-
guished (Getz & Owen-Smith, 2011; Hassell, 1975;
Krebs, 2014; Nicholson, 1954). Scramble competition
involves equal or random resource partitioning among all
of the individuals in the population. In contrast, contest
competition implies monopolistic resource utilization by
the successful competitors among the local population; the
remaining resources (i.e., leftovers) are insufficient to sus-
tain more individuals, which ultimately fail to reproduce
and are forced to leave the area, live as floaters, or die
from the lack of necessary resources. Although some
authors (e.g., Rockwood, 2015) have suggested that divid-
ing intraspecific relationships into either scramble or
contest competition may not be realistic, the assumptions
of scramble and contest competition have been applied
in studies of a broad range of organisms such as
fish (Petersson & Järvi, 2000), birds (Rieucau &
Giraldeau, 2008), and mammals (Dammhahn & Kappeler,
2010; Luo et al., 2020). Sibly et al. (2007) examined
634 populations of mammals, birds, fish, and insects. Their
results indicate that contest competition, rather than
scramble competition, reflects the status of the species
examined. However, some of the existing studies, which
consider, for example, birds (Black-throated blue warblers;
Rodenhouse et al., 2003) and ungulates (Tibetan antelopes;
Luo et al., 2020), conclude that certain species exhibit
scramble competition. There is a large amount of accumu-
lated empirical and theoretical knowledge regarding cases
in which primates exhibit both types of intraspecific com-
petition (e.g., Howlett & Wheeler, 2021; Isbell, 1991; Jan-
son & Van Schaik, 1988; Kappeler & van Schaik, 2002;
Knott et al., 2008; van Schaik, 1989). To the best of the
authors' knowledge, the applications of intraspecific com-
petition models for large carnivores are limited.
Phenomenological discrete time models Nt+1 = f(Nt),
which predict population size at time step t + 1
(e.g., reproductive season or year) as a function of the pop-
ulation size at the previous time step t according to the
scramble or contest competition scenarios, can be derived
from Verhulst's logistic equation (Johst et al., 2008;
Royama, 1992) or by other means (e.g., Hassell, 1975).
Meaningful derivations are obtained according to pre-
defined fundamental processes, resource partitioning, and
individual interaction (e.g., Anazawa, 2019; Eskola &
Geritz, 2007; Geritz & Kisdi, 2004). The site-based frame-
work provides general underlying principles for develop-
ing population models that predict population size within
a network of sites based on initial number of individuals
and the expected number of remaining individuals and
offspring, which is determined according to explicit
assumptions about intraspecific competition and individ-
ual clustering (Anazawa, 2009; Brännström &
Sumpter, 2005; Johansson & Sumpter, 2003). Both the dis-
tribution of individuals in habitat sites (random or aggre-
gated) and their interactions (scramble or contest) are
incorporated in the site-based framework. Four such
models that employ only two parameters, which can be
expressed as functions of the intrinsic growth rate and car-
rying capacity (i.e., equilibrium density), are listed in
Table 1. These models pertain to scramble or contest com-
petition and assume Poisson or geometric probability dis-
tribution for the number of individuals at the sites.
Regarding intraspecific competition, the Ricker (1954) and
Hassell (1975) models reflect scramble competition while
the Skellam (1951) and Beverton–Holt (Beverton &
Holt, 1957) models reflect contest competition. Regarding
spatial patterns, the Ricker and Skellam models are associ-
ated with random (i.e., Poissonian) individual distribution
while the Hassell and Beverton–Holt models assume the
aggregation of individuals according to the geometric dis-
tribution and a considerable proportion of the unoccupied
sites within the area.
This study aims to examine the consistency of intra-
specific carnivore population dynamics using scramble
and contest population models. This knowledge may be
relevant when justifying the choice of a particular popu-
lation model to predict the population dynamics of a cer-
tain species. For example, the Ricker model has typically
been applied in various cases (Hall, 1988; Hamada, 2024;
Sˇuba et al., 2023), but territorial animals, for example,
gray wolves Canis lupus, are generally supposed to reflect
contest competition (Mech & Boitani, 2003). Such an
examination is of particular interest as both carnivore
species reproduce seasonally but differ in their social
structure. Therefore, we used model selection and multi-
2 SˇUBA ET AL.
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
model inference to investigate whether the population
dynamics of carnivore species are more consistent with
contest competition. We used available census and cul-
ling (i.e., hunting) data on the gray wolf and Eurasian
lynx Lynx lynx populations in Latvia obtained over the
last six decades. Additional assumptions were also con-
sidered. As the prey populations for both species in
Europe have been increasing in recent decades
(Apollonio et al., 2010; Carpio et al., 2021), we investi-
gated whether an increase in carrying capacity would
contribute to the growth of the local wolf and lynx popu-
lations independently of their intrinsic growth rates. Pre-
vious studies on wolf and lynx dynamics in Latvia
(Kawata, 2008), and on wolf dynamics in neighboring
Lithuania (Balcˇiauskas & Kawata, 2009), have found that
a time lag due to delayed maturity was significant for the
wolf population growth. Therefore, we also considered
model versions with different values for the time lag
operator. Lastly, it is known that observation noise and
bias in time series data can affect estimates of density
dependence and, consequently, carrying capacity if not
explicitly accounted for (Knape, 2008). As the game cen-
sus data are often subject to overestimation and contain
process errors of unknown magnitude, we also investi-
gated the effect of abundance estimation bias and the
additional noise of known magnitude on parameter
estimates.
2 | MATERIALS AND METHODS
2.1 | Study species
In Europe, the distribution ranges of the gray wolf and Eur-
asian lynx populations are currently largely overlapping
(Linnell et al., 2008). Wolves have been capable of recover-
ing after drastic reductions by natural recolonization
(Hayes & Harestad, 2000; Linnell et al., 2005; Pletscher
et al., 1997; Reinhardt et al., 2019; Szewczyk et al., 2019;
Vilà et al., 2003) or artificially facilitated acclimatization
(Smith et al., 2003). Reintroduction efforts have played a
more significant part in the lynx recolonization of several
areas in central and western Europe due to the inferior dis-
persal capabilities of lynx (Breitenmoser, 1998;
Breitenmoser et al., 1998; von Arx et al., 2004). The social
structures of gray wolves and Eurasian lynx are different.
The wolves form family packs which consist of a monoga-
mous pair, their offspring before reproductive age, and new-
born pups (Mech, 1970). The lynx males and
nonreproductive females are solitary, while reproductive
females accompany and nurse their cubs until they become
independent (Breitenmoser et al., 2000). Wolf packs defend
their territories from intruding conspecifics (Jędrzejewska &
Jędrzejewski, 1998; Mech, 1970), whereas the home ranges
of neighboring lynx males or members of opposite sexes
may overlap to some extent (Schmidt et al., 1997).
The so-called Baltic wolf and lynx populations, to
which the Latvian sub-populations belong, are thriving
wolf and lynx populations in Europe (Chapron
et al., 2014). In Latvia, wolves and lynx have never been
totally exterminated (Bagrade et al., 2016; Ozolin¸sˇ, Bagrade,
et al., 2017; Ozolin¸sˇ, Zˇunna, et al., 2017; Sˇuba et al., 2021).
The implemented hunting praxis for wolves and lynx has
varied historically (Ozolin¸sˇ et al., 2008, 2016). Before Latvia
became a member state of the European Union in 2004,
wolves were legally hunted all year round (Ozolin¸sˇ,
Bagrade, et al., 2017). After joining the EU, a closed season
and an annual quota were introduced in compliance with
the Habitat Directive (European Council, 1992). Lynx have
always been a protected species in Latvia, and illegal kill-
ings have been penalized. Moderate hunting and lethal
control of this species have been conducted according to
strict supervision by game surveillance authorities
(Bagrade et al., 2016; Ozolin¸sˇ, Zˇunna, et al., 2017). In the
hunting season of 2021/2022, the hunting of lynx was no
longer permitted, and a permanent hunting ban on this
species was introduced in 2022 (Ozolin¸sˇ et al., 2022). In our
analysis, we focused on annual transitions, which are not
affected by seasonal variations in hunting activity. The
poaching of wolves and lynx in Latvia is rare, and its share
within the total number of intentionally removed wolf and
lynx individuals can be considered to be negligible as the
majority of hunters are motivated to report the killings in
compliance with approved legal culling (J. Ozolin¸sˇ, per-
sonal communication).
2.2 | Abundance estimates and
culling data
The datasets used in this study contained the estimated
population (i.e., number of individuals) of wolves and
TABLE 1 Two-parameter (αi and βi) single-species population
models derived according to site-based framework.
Model
name Expression
Competition
type
Distribution
of
individuals
at sites
Ricker αRNtexp(
βRNt) Scramble Random
Skellam αS(1 
 exp[–βSNt]) Contest Random
Hassell αHNt/(1 + βHNt)
2 Scramble Aggregated
Beverton–
Holt
αBHNt/(1 + βBHNt) Contest Aggregated
SˇUBA ET AL. 3
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
lynx and the numbers of culled individuals in Latvia from
1958 to 2022. Such data are available on the Latvian State
Forest Service website (www.vmd.gov.lv) and in previous
publications (e.g., Andersone-Lilley & Ozolins, 2005;
Kawata, 2008; Vanags, 2010). The numbers of wolves and
lynx are estimated for every regional game management
district at the closure of the hunting season, that is,
31 March or 30 April, and added together to calculate the
total abundance over larger game management territories
and the whole country. This estimate can be understood
as the abundance before reproduction. However, it is
prone to overestimation as the home ranges of large car-
nivores are likely to encompass several adjacent hunting
districts (e.g., Herfindal et al., 2005), and the individual
animals may have been accounted for more than once.
The extent of this overestimation bias is unknown. Other
studies (Ozolin¸sˇ, Bagrade, et al., 2017; Ozolin¸sˇ, Zˇunna,
et al., 2017; Sˇuba et al., 2021) suggest that the actual
abundance differs from the reported number by less than
50%. To examine the potential effects on model perfor-
mance caused by biased population estimates, we investi-
gated the tested models with adjusted abundance
estimates using Xt = cNt, where the reported
abundance estimates Nt were multiplied by a factor of
c = 0.7, 0.8, 0.9, or 1.1 to account for potential bias, or by
setting c = 1 if no bias was assumed.
The effects of potential observation errors in the
abundance data on the model statistics and parameter
values were investigated using 5000 simulated datasets
per species. These datasets were created by multiplying
the original wolf and lynx data by a set of independent
log-normally distributed annual random numbers, with
1 as the expected value (parameters μ = 
½σ2noise and
σ2 = σ2noise), and by setting σnoise to 0.1. While we do not
know the true level of noise in the data, this reflects a
realistic level, which at least helps to identify the occur-
rence and direction of bias due to observation errors and
also gives a rough clue about the level of bias in the dif-
ferent parameters. The parameter estimates were also
bias-corrected by subtracting the levels of bias (estimated
in the simulations) from the originally estimated parame-
ter values. The approach used here for exploring the
effect of simulated additional observation error, and
using that for bias correction, applies the idea of simula-
tion extrapolation (i.e., SIMEX), which is thoroughly
explicated in Cook and Stefanski (1994) and Stefanski
and Cook (1995).
2.3 | Models and parameters
In this study, we applied phenomenological single-
species discrete time population models that use
parameters related to intrinsic growth rate, carrying
capacity, and residual variance. The analyses were
carried out separately for each species. The intrinsic
growth rate was parameterized as an instantaneous
(i.e., logarithmic) rate r0 instead of a finite (i.e., geomet-
ric) rate λ0 (Skalski et al., 2005). The effects on population
growth following the removal of known individuals due
to culling, in relation to the abundance estimate
before breeding and culling Xt, were accounted for by
adding the total number of culled individuals ht to the
estimated abundance at the closure of the hunting season
Xt+1, obtained in the spring of the following calendar
year, that is, Xt+1 + ht = f(Xt) or Xt+1 = f(Xt) 
 ht
(Eberhardt, 1987). The annual logarithmic per capita
growth rate was calculated as rt = ln(Xt+1 + ht) 
 ln(Xt).
A monotonic multiplicative change in the projected car-
rying capacity over the study period was allowed by
expressing Kt as a function of year t and using the two
parameters a and b; the latter expressed the temporal
trend as
Kt ¼ exp aþbtð Þ: ð1Þ
Hence, the abundance of wolves or lynx during the
temporal transition from year t to year t + 1 was modeled
according to the focal species' abundance before breeding
Xt, the number of culled individuals ht, the intrinsic
growth rate r0, the parameters a and b (which defines the
carrying capacity at year t), and the stochastic parameter
εt, quantifying the residual variation due to errors,
uncontrolled environmental factors, and demographic
stochasticity:
Xtþ1¼ f Xt, tjr0,a,bð Þeεt 
ht: ð2Þ
The error term εt was treated as a normally distrib-
uted random variable (see Ruokolainen et al., 2009). Set-
ting its mean and variance as 
½σε2 and σε2,
respectively, ensures that the expected value of a log-
normally distributed ratio Xtþ1þhtf Xt , tjr0,a,bð Þ equals unity and
therefore provides unbiased parameter estimates. Fur-
thermore, although the model itself is a first-order auto-
regression, we allowed the residuals to be serially
autocorrelated according to a first-order linear autore-
gressive process with the autocorrelation coefficient φ.
A significant proportion of the wolf and lynx popula-
tions in Latvia consists of juveniles and subadults
(Bagrade et al., 2016; Sˇuba et al., 2021). Hence, a time lag
operator τ, represented by a non-negative integer, was
applied to account for the effect of delayed reproduction
due to maturity. In previous studies, 0 < τ ≤ 3 improved
the model fitting to the observed wolf dynamics in Latvia
and Lithuania (Balcˇiauskas & Kawata, 2009;
4 SˇUBA ET AL.
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Kawata, 2008). A time lag of τ = 0, 1, 2, 3, 4 was intro-
duced into the Ricker, Hassell, and Beverton–Holt
models (see Equations 3, 5, and 6); this introduction was
similar to that carried out in Kawata (2008) and
Balcˇiauskas and Kawata (2009). The Skellam model
requires only a single term with the explanatory variable
(Table 1). Hence, the introduction of τ > 0 into this
model was incompatible with a first-order autoregressive
process. Rearrangements of the Ricker, Skellam, Hassell,
and Beverton–Holt equations are given in
Equations (3–6), and a description of their derivation is
given in Appendix 1 (a more detailed description of the
modified Skellam model is also provided by Sˇuba
et al., 2023).
f R Xt,Xt
τ, tjr0,a,b,τð Þ¼Xt exp r0 1

Xt
τ
exp aþ t
 τ½ bð Þ

  
:
ð3Þ
f S Xt, tjr0,a,bð Þ¼ exp aþbtð Þ
1
 exp 
Xt 1þe

r0W 
exp r0
er0f gð Þ
eaþbt
r0
 
1þ e
r0W 
exp r0
 er0f gð Þ :
ð4Þ
fH Xt,Xt
τ, tjr0,a,b,τð Þ¼
Xt exp r0ð Þ
1þ exp r02
 	
1
  Xt
τexp aþ t
τ½ bð Þ
 2 :
ð5Þ
f BH Xt,Xt
τ, tjr0,a,b,τð Þ¼
Xt exp r0ð Þ
1þ exp r0ð Þ
1½  Xt
τexp aþ t
τ½ bð Þ
:
ð6Þ
2.4 | Fitting procedure and model
selection
To fit the investigated models to the data and to find the
maximum likelihood estimates for r0, a, b, σε
2, and φ,
numerical optimization and generalized least squares
(GLS) techniques were applied. Nonlinear numerical
optimization using the Nelder–Mead simplex algorithm
was used to estimate the parameters φ, r0, or b, which
depended on the particular model. The remaining param-
eters could be estimated analytically—during each itera-
tion of the optimization algorithm—using the GLS
techniques. The discrete parameter τ was given by the
model setting, that is, its value was assessed in terms of
model uncertainty.
To assess the consistency of the investigated models
with the data, we applied model selection and assessed
individual model weights (Aho et al., 2014; Symonds &
Moussalli, 2011; Wagenmakers & Farrell, 2004). For each
species, we evaluated a set of 16 candidate models with
optimized parameter values. The set included the modifi-
cations of the Ricker, Hassell, and Beverton–Holt models
described above, each with the time lag τ of 0, 1, 2, 3, or
4, as well as the Skellam model. In our case, as the num-
ber of parameters was the same for every model (k = 5),
the penalty terms of the model complexity would be con-
stant in information criteria such as AIC(c) and BIC, and
the differences in the performances of the models would
in all cases be determined solely by differences in 
2ln
(L) = ln(L
2), where L is the maximized likelihood func-
tion value. Hence, model selection results that are deter-
mined according to the Bayesian or Akaike information
criterion provide essentially the same differences between
the models, that is, ΔBIC = ΔAIC(c) = Δln(L
2).
The standard errors (SEs) for the estimated model
parameters (r 0^, â, b ,^ σ ε^
2, and φ )^ were calculated by
presuming that the inverse Hessian matrices of the nega-
tive log-likelihood function provided the variance–
covariance matrix of the parameters. Confidence inter-
vals (95%) were constructed by adding 1.96  SE to or
subtracting it from the parameter estimates. All the cal-
culations were conducted using the programming envi-
ronment R (R Core Team, 2017, see Supporting
Information for the code).
3 | RESULTS
The Latvian wolf and lynx abundance data exhibited an
overall increasing trend with considerable fluctuations
over the last six decades (Figure 1a,c), but with a decrease
in the calculated annual growth rate over time
(Figure 1b,d). The mean observed annual per capita
growth rates r¯t for the wolf and lynx populations in
Latvia were 0.577 (SD = 0.425, n = 58) and 0.225
(SD = 0.260, n = 58), respectively.
The estimates of the intrinsic growth rate, the param-
eters a and b that determined the projected carrying
capacity Kt = exp(a + bt), the residual variance σε
2, and
the autocorrelation coefficient φ according to the investi-
gated population models are given in Tables 2 and 3. The
unconditional parameter estimates according to full-
model averaging by the model weights w are given in
Table 4. The model-averaged estimates of the residual
parameters were σ ̄ε2 = 0.0387 and φ̄ = 0.5939 for the
wolf population models and σ̄ε
2 = 0.0386 and φ̄ = 0.2671
for the lynx population models, respectively. The
Beverton–Holt model with τ = 0 was found to be the best
performing model, with parameter values of r 0^ = 1.3903
(SE = 0.1509), â = 5.9272 (SE = 0.3458), and b^= 0.0245
(SE = 0.0072) for the wolves and r 0^ = 0.7581
SˇUBA ET AL. 5
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
(SE = 0.1234), â = 5.3668 (SE = 0.2844), and b^= 0.0381
(SE = 0.0062) for the lynx. The most obvious difference
between the parameter estimates for the two species was
a higher intrinsic growth rate for the wolves compared to
that of the lynx. The model fit with the data and the trend
in the projected carrying capacity of studied wolf and
lynx populations are illustrated in Figures 2a–h and
3a–h, respectively.
For the populations of both carnivore species, the
models that implied contest competition (i.e., the
FIGURE 1 Population
dynamics of wolves (top panels
a, b) and lynx (bottom panels c,
d) in Latvia during the last six
decades in terms of pre-breeding
and pre-culling abundance Nt
(panels a, c) and annual per
capita growth rate rt =
ln(Nt+1 + ht) 
 ln(Nt) (panels
b, d).
TABLE 2 Summary of estimated intrinsic growth rate r0 and parameters determining the projected carrying capacity Kt = exp(a + bt),
residual variance, autocorrelation coefficient, and information-theoretic statistics of tested models describing wolf population dynamics in
Latvia from 1963 to 2021.
Model τ r 0^ â b^ σ ε^
2 φ^ ln(L) Δln(L
2) w
Ricker 0 1.092 6.169 0.019 0.053 0.693 21.1 8.1 0.013
Ricker 1 1.068 6.174 0.019 0.051 0.651 19.3 11.6 0.002
Ricker 2 0.970 6.587 0.013 0.063 0.689 16.0 18.3 <0.001
Ricker 3 1.027 6.143 0.021 0.060 0.697 17.9 14.6 0.001
Ricker 4 0.890 7.012 0.007 0.082 0.768 15.3 19.6 <0.001
Skellam 0 1.249 5.946 0.023 0.044 0.640 22.9 4.5 0.078
Hassell 0 1.228 5.997 0.023 0.045 0.642 22.8 4.6 0.074
Hassell 1 1.190 6.057 0.021 0.042 0.574 20.8 8.6 0.010
Hassell 2 1.077 6.358 0.017 0.050 0.614 16.8 16.7 <0.001
Hassell 3 1.129 6.129 0.021 0.052 0.645 18.8 12.6 0.001
Hassell 4 0.942 7.005 0.007 0.075 0.742 15.6 19.0 <0.001
Beverton–Holt 0 1.390 5.927 0.025 0.037 0.593 25.1 0.0 0.743
Beverton–Holt 1 1.321 6.042 0.022 0.035 0.489 22.8 4.7 0.072
Beverton–Holt 2 1.244 6.203 0.019 0.041 0.488 18.3 13.7 0.001
Beverton–Holt 3 1.222 6.221 0.019 0.045 0.605 20.0 10.3 0.004
Beverton–Holt 4 1.000 6.991 0.007 0.068 0.727 16.0 18.3 <0.001
6 SˇUBA ET AL.
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Beverton–Holt and Skellam models) had a greater cumula-
tive model weight (0.899 for the wolf population and 0.679
for the lynx population) compared to the models of the
scramble competition (i.e., the Hassell and Ricker models).
The 95% confidence set of models describing wolf popula-
tion dynamics included the Beverton–Holt model with
τ = 0 (w = 0.743), the Skellam model (w = 0.078), the
Hassell model with τ = 0 (w = 0.074), and the Beverton–
Holt model with τ = 1 (w = 0.072). For the lynx popula-
tion, the 95% confidence set included the Beverton–Holt
model with τ = 0 (w = 0.471), the Hassell model with
τ = 0 (w = 0.187), the Skellam model (w = 0.163), the
Ricker model with τ = 0 (w = 0.088), the Beverton–Holt
model with τ = 1 (w = 0.026), the Hassell model with
τ = 1 (w = 0.013), and the Beverton–Holt model with
τ = 3 (w = 0.012). The models that assumed a random
distribution of individuals at sites (i.e., the Ricker and
Skellam models) had a lower cumulative model weight
(0.093 for the wolf population and 0.269 for the lynx popu-
lation) than the models involving non-random distribution
and individual clustering (i.e., the Hassell and Beverton–
Holt models). The models that contained a time lag opera-
tor τ (i.e., the models with τ > 0) had a low cumulative
model weight (0.092 for the wolf population and 0.091 for
the lynx population). Hence, in further analyses, we
focused on models with τ = 0.
Multiplying the abundance estimates by a bias factor
c = {0.7, 0.8, 0.9, 1, 1.1} indicated that overestimation
would have influenced the relative performances of the
candidate models (Figure 4a–f). When overestimation
was assumed for the lynx population, the model weights
of the Beverton–Holt model with τ = 0 were greater com-
pared to the assumption of no bias or when assuming
underestimation of the abundance data (Figure 4d); these
results therefore provided robust support in favor of this
model. The opposite pattern was observed for the worst
performing models (Figure 4e,f). For the wolf population,
the weights of the models with assumed bias in the abun-
dance data changed to a lesser extent (Figure 4a–c), again
suggesting that the results from the model selection
would, in practice, be the same regardless of the bias in
the abundance estimates.
The potential effect of observation errors was investi-
gated by testing the models' performance in 5000
TABLE 3 Summary of estimated intrinsic growth rate r0 and parameters determining the projected carrying capacity Kt = exp(a + bt),
residual variance, autocorrelation coefficient, and information-theoretic statistics of tested models describing lynx population dynamics in
Latvia from 1963 to 2021.
Model τ r 0^ â b^ σ ε^
2 φ^ ln(L) Δln(L
2) w
Ricker 0 0.645 5.422 0.037 0.040 0.294 13.3 3.4 0.088
Ricker 1 0.581 5.562 0.035 0.042 0.207 10.9 8.3 0.007
Ricker 2 0.555 5.603 0.034 0.044 0.230 10.0 10.0 0.003
Ricker 3 0.571 5.487 0.036 0.044 0.296 10.9 8.3 0.007
Ricker 4 0.464 6.037 0.026 0.049 0.313 8.1 13.8 <0.001
Skellam 0 0.687 5.405 0.037 0.039 0.280 13.9 2.1 0.163
Hassell 0 0.693 5.403 0.037 0.039 0.277 14.1 1.9 0.187
Hassell 1 0.619 5.539 0.035 0.041 0.184 11.4 7.2 0.013
Hassell 2 0.583 5.601 0.034 0.043 0.214 10.3 9.4 0.004
Hassell 3 0.597 5.496 0.036 0.043 0.288 11.1 7.8 0.010
Hassell 4 0.473 6.089 0.026 0.049 0.309 8.2 13.6 0.001
Beverton–Holt 0 0.758 5.367 0.038 0.037 0.263 15.0 0.0 0.471
Beverton–Holt 1 0.666 5.510 0.036 0.040 0.159 12.1 5.8 0.026
Beverton–Holt 2 0.614 5.602 0.035 0.042 0.197 10.7 8.6 0.006
Beverton–Holt 3 0.625 5.511 0.036 0.043 0.281 11.4 7.3 0.012
Beverton–Holt 4 0.482 6.152 0.025 0.048 0.305 8.3 13.5 0.001
TABLE 4 Unconditional estimates of intrinsic growth rate r0
and parameters determining the projected carrying capacity
Kt = exp(a + bt) for wolf and lynx population dynamics in Latvia
from 1963 to 2021 (the estimates are model-averaged values of the
Ricker, Skellam, Hassell, and Beverton–Holt models according to
model weights; numbers in brackets are parameter standard
errors).
Population r¯0 a¯ b¯
Wolf 1.354 (0.160) 5.949 (0.357) 0.024 (0.007)
Lynx 0.711 (0.123) 5.400 (0.289) 0.037 (0.006)
SˇUBA ET AL. 7
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
simulated datasets, where a log-normally distributed
noise was added to the original data. The performance of
the Beverton–Holt model for the wolf time series
remained superior according to the mean values of the
model weights (Figure 5a–c). In the lynx time series,
however, the mean value of the model weights consider-
ably decreased for the Beverton–Holt and Hassell models
and increased for the Ricker and Skellam models
(Figure 5d–f).
The assumptions about bias in the abundance data
affected the estimates of the intrinsic growth rate and the
constant parameter defining the carrying capacity
(Figure 6a–f). As the bias factor was increased from 0.7 to
1.1, the intrinsic growth rate decreased for both wolves
and lynx. In contrast, the constant parameter â (level
of Kt) increased for lynx, but less so for wolves. The trend
parameter b^(a logarithmic trend in Kt) remained largely
unaffected by bias in the abundance estimates, particu-
larly when related to the parameter uncertainty that was
measured as the confidence interval range.
The introduction of noise into the original data
mainly increased the error variance and decreased the
autocorrelation (Figure 7a–j). It also had a considerable
effect on the parameter estimates that determined the
growth rate and the carrying capacity of the wolf popula-
tion. In particular, the added noise increased the esti-
mates of parameters r 0^ and b^and decreased that of the
parameter â. However, the additional noise had no signif-
icant effect on the parameter estimates of the lynx
population.
FIGURE 2 Fit of the tested
scramble and contest
competition models (solid line)
with temporally increasing
carrying capacity (dashed line)
to the pre-culling abundance
estimates (left panels) and
corresponding observations and
estimates of instantaneous per
capita growth rate (right panels)
of the wolf population in Latvia
during the last six decades:
Ricker model (panels a, b),
Skellam model (panels c, d),
Hassell model (panels e, f), and
Beverton–Holt model (panels
g, h).
8 SˇUBA ET AL.
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
4 | DISCUSSION
Gray wolves and Eurasian lynx, as with large carnivores
in general, are territorial animals, albeit with different
social structures. The differences in individual interaction
are expected to conform to the different functional forms
of the models describing population dynamics. The main
difference between the scramble and contest competition
models, as determined by their functional forms is related
to the predicted response when the expected abundance
Nt+1 = f(Nt) exceeds the carrying capacity of the initial
abundance Nt due to an influx of individuals from else-
where. In scramble competition, under such circum-
stances, the resulting abundance decreases due to density
dependence and overcompensation, whereas in contest
competition, the abundance asymptotically approaches a
certain level as individuals attack each other and the
excess individuals are killed or forced to leave the sites
(e.g., Brännström & Sumpter, 2005; Eskola &
Geritz, 2007).
In this study, we investigated the consistency of the
time series data with a range of phenomenological
single-species population models that are derivable from
first principles (i.e., a mechanistic bottom-up manner of
deriving population models based on the process at
an individual level) according to site-based frame-
works (Anazawa, 2009; Brännström & Sumpter, 2005;
Johansson & Sumpter, 2003). The models were rear-
ranged to contain the same set of parameters, namely,
the intrinsic growth rate and the carrying capacity. The
Ricker and Beverton–Holt population models assume
scramble competition and a random distribution of
FIGURE 3 Fit of the tested
scramble and contest
competition models (solid line)
with temporally increasing
carrying capacity (dashed line)
to the pre-culling abundance
estimates (left panels) and
corresponding observations and
estimates of instantaneous per
capita growth rate (right panels)
of the lynx population in Latvia
during the last six decades:
Ricker model (panels a, b),
Skellam model (panels c, d),
Hassell model (panels e, f), and
Beverton–Holt model (panels
g, h).
SˇUBA ET AL. 9
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
FIGURE 4 Effect of overestimation (c < 1) and underestimation bias (c > 1) in abundance estimates Xt = cNt on the weights of the
candidate models describing wolf (top panels a–c) and lynx (bottom panels d–f) population dynamics: Beverton–Holt model (panels a, d),
Hassell and Skellam models (panels b, e), Ricker model (panels c, f); τ values correspond to introduced time lags from 0 to 4 years to account
for the effect of delayed reproduction due to maturity.
FIGURE 5 Effect of
introduced noise in the data in
5000 simulations, multiplying
original values by a log-
normally distributed variable
with μ = 
½σ2noise and
σ2 = σ2noise set to 0.1, on the
weights of the candidate models
describing wolf (top panels a–c)
and lynx (bottom panels d–f)
population dynamics: Beverton–
Holt model (panels a, d), Hassell
and Skellam models (panels
b, e), Ricker model (panels c, f);
τ values correspond to
introduced time lags from 0 to
4 years to account for the effect
of delayed reproduction due to
maturity.
10 SˇUBA ET AL.
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
individuals at habitat sites and contest competition and
individual aggregation, respectively. The Hassell and
Skellam models assume other combinations of underly-
ing assumptions, namely scramble competition and indi-
vidual aggregation and contest competition and a
random distribution of individuals (Table 1). Our study
confirmed the different performances of the investigated
models in describing the population dynamics of the wolf
and lynx population in Latvia between 1963 and 2021, as
well as the superior performance of the other models
(Tables 2 and 3) when compared to the frequently
applied Ricker model (Hamada, 2024). Specifically, the
changes in wolf abundance were more consistent with
the models that assumed contest competition and spatial
aggregation and matched the expectations related to their
territorial behavior and assemblage into packs. For lynx,
the Beverton–Holt model with τ = 0 also had the highest
model weights, but differed from that of the other models
to a lesser extent. Moreover, the differences between the
contest and scramble competition scenarios were less
clear when the other models were considered. This may
be explained by the fact that the territorial behavior of
lynx is not as extensive and competitive as it is for wolves.
Although individuals in neighboring territories tend to
avoid each other, lynx males may share 30%–75% of their
home ranges with other individuals (Schmidt
et al., 1997), whereas wolf packs violently attack and pur-
sue intruding conspecifics, and the overlapping of neigh-
boring territories is rare (Mech, 1970).
The investigated population models pertain to some-
what extreme cases of individual spatial distribution
(i.e., random or extremely clustered cases). The inclusion
of a parameter that quantifies spatial clustering
(e.g., Anazawa, 2009; Brännström & Sumpter, 2005) was
deliberately avoided due to the general scope of this study
and the focus on parameters associated with intrinsic
growth rate and carrying capacity. However, the differ-
ences in individual aggregation are undoubtedly signifi-
cant enough to require a quantitative scaling. Moreover,
ideal scramble and ideal contest competition can be
regarded as extreme cases in intraspecific competition.
Anazawa (2009) has developed a general population
model with a parameter β that quantifies deviation from
the ideal contest competition scenario and ranges
between zero and one (i.e., β = 0 for the ideal contest
competition and β = 1 for the ideal scramble
FIGURE 6 Effect of overestimation (c < 1) and underestimation bias (c > 1) in abundance estimates Xt = cNt on estimated model
parameters describing wolf (top panels a–c) and lynx (bottom panels d–f) population dynamics: Intrinsic growth rate r 0^ (panels a, d) and
parameters â (panels b, e) and b^(panels c, f) determining the projected carrying capacity Kt = exp(a + bt) (whiskers indicate 95%
confidence intervals).
SˇUBA ET AL. 11
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
competition). A general population model with addi-
tional parameters may be worth exploring in further
studies. This may be relevant when considering the vari-
ous possibilities of social aggregation and territorial
behavior since the lynx case in this study demonstrated
less conformity with the ideal contest competition
scenario.
A viable population is expected to grow at densities
below the carrying capacity, which is frequently treated
as a constant. This was an underlying assumption in pre-
vious investigations of the wolf and lynx populations in
Latvia (Kawata, 2008) and of the wolf population
in Lithuania (Balcˇiauskas & Kawata, 2009). The previous
estimates of the carrying capacities for wolves and lynx in
Latvia, which used the same dataset until the year 2005
but applied a quadratic population growth model, were
1066–1092 and 971–1188 individuals, respectively
(Kawata, 2008). In this study, however, the projected car-
rying capacity or equilibrium density Kt was allowed to
increase, as suggested by the increase in prey abundance
observed locally (Ozolin¸sˇ et al., 2016) and regionally
(Apollonio et al., 2010; Carpio et al., 2021). Prey biomass
has been called the most important factor affecting wolf
population dynamics (Fuller et al., 2003). For lynx, an
increasing prey density would allow smaller home ranges
(Herfindal et al., 2005). We assumed an exponential func-
tion according to the growth process and defined the car-
rying capacity as Kt = exp(a + bt). In the majority of
models, zero was beyond the 95% confidence limits of the
parameter b, supporting the validity of the assumption
regarding the increase in the carrying capacity
(Figure 6c,f). Consequently, the estimates of intrinsic
growth rate r0 were higher than the estimates provided
by Kawata (2008), which assumed a constant carrying
capacity (1.354 vs. 0.962 for the wolf population and
0.711 vs. 0.399 for the lynx population).
Apart from the increasing trend, the population
dynamics of both species also exhibited a fluctuating pat-
tern (Figure 1a,c), which may have been caused by
changes in the carrying capacity (e.g., fluctuations in prey
availability) or stochastic factors. Moreover, in the over-
compensatory models, oscillations between stable states
and chaotic dynamics may occur with elevated growth
rates (e.g., May, 1974, 1976). As this is not the case for con-
test competition, which limits the population density, the
fluctuations in population density under the assumption
of contest competition are expected according to changes
in the limiting factors associated with the carrying capac-
ity. Although a temporally increasing and oscillating carry-
ing capacity may be introduced into the models by
applying a periodical function and adding parameters that
determine the period and range of the oscillations, we
treated the fluctuations in population density as autocorre-
lated noise to avoid overparameterization of the models.
The model-averaged intrinsic growth rates r¯0 for the
wolf and lynx population in Latvia were 1.354 and 0.711
FIGURE 7 Effect of introduced noise in the data in 5000 simulations, multiplying original values by a log-normally distributed variable
with μ = 
½σ2noise and σ2 = σ2noise set to 0.1, on estimated model parameters describing wolf (top panels a–e) and lynx (bottom panels f–j)
population dynamics: Intrinsic growth rate r 0^ (panels a, f) and parameters â (panels b, g) and b^(panels c, h) determining the projected
carrying capacity Kt = exp(a + bt) as well as error variance σ ε^
2 (panels d, i) and autocorrelation coefficient φ^(panels e, j).
12 SˇUBA ET AL.
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
(Table 4), respectively. These rates are associated with
stable population dynamics and would correspond to the
geometric growth rates (λ0) of 3.87 and 2.04, respectively.
The intrinsic growth rates of both species are expected to
differ due to different fecundity. The mean litter sizes of
gray wolves and Eurasian lynx in Europe vary from 4.4 to
7.7 (on average, 5.9, Jędrzejewska et al., 1996) and from
2.08 to 2.12 (Nilsen et al., 2012), respectively. In Latvia,
prenatal fecundity for wolves and lynx during the last
two decades has been 6.1 (Sˇuba et al., 2021) and 3.2
(Bagrade et al., 2016), respectively. This matches well
with the difference in growth rates observed between the
two species. Assuming mean wolf and lynx litter sizes of
5.9 and 2.1, which may account for postnatal mortality,
equal sex ratio, and the complete survival of the recruited
individuals, the annual geometric growth rates of 3.87
and 2.04 could be achieved if 97.4% of all wolf females
and 98.7% of all lynx females reproduced. This estimate is
based on the assumption that the proportion of reprodu-
cing females in the population under the stated circum-
stances equals 2(λ0 
 1)/y, where y corresponds to the
mean number of offspring per female (see Kelker, 1947).
In this study, the mean observed annual growth rates r¯t
were 0.577 for wolves and 0.225 for lynx. These rates cor-
respond to geometric growth rates λt of 1.78 and 1.25,
respectively. As reproductive output decreases due to
density dependence, the litter size or the proportion of
reproductive females may remain high to counterbalance
additional mortality. During the last two decades, 63.1%
and 87% of the wolf and lynx adult females culled in
Latvia contained traces of reproduction (Bagrade
et al., 2016; Sˇuba et al., 2021). Such proportions were
lower than expected in the absence of competition, but
comparable to other areas with considerable hunting
pressure (Frank & Woodroffe, 2001; Fuller et al., 2003;
L
opez-Bao et al., 2019; Mech et al., 2016).
Wolf and lynx individuals usually become mature
later than their first or second year of life (Fuller
et al., 2003; Nilsen et al., 2012). Therefore, the effect of
delayed reproduction on the population growth rate has
been considered in previous studies on wolf and lynx
population dynamics by introducing a time lag operator τ
(Balcˇiauskas & Kawata, 2009; Kawata, 2008). This
approach estimates the growth rate at year t by consider-
ing only mature individuals, that is, all the individuals at
year t 
 τ would have become mature by the year t and
capable of producing offspring. In this study, we applied
a time lag operator for the models where meaningful
introduction of τ was possible (i.e., for all the models but
the Skellam model). Unlike in previous research, the
models with τ > 0 had significantly lower model weights
than the models with τ = 0. Conversely, in previous stud-
ies, the quadratic model describing wolf dynamics was
the most consistent with the data in the case of τ = 3
(Balcˇiauskas & Kawata, 2009, Kawata, 2008). As men-
tioned above, the previous studies also differed from the
present one in that they assumed a constant carrying
capacity. As different results may arise from different
underlying assumptions, our findings do not dismiss the
potential significance of delayed maturity in modeling
wolf dynamics altogether. Rather they demonstrate that
it had an insignificant impact on determining wolf popu-
lation growth under the modeling circumstances applied
here, which included an unstructured model, a carrying
capacity with a temporal trend, and autocorrelated
residuals.
As the abundance data used in this study likely
included a considerable overestimation bias, we investi-
gated whether such bias could have influenced the model
performance. Assuming overestimation up to 42.9% and
an unlikely underestimation down to 
9.1%, the abun-
dance estimates for both species were multiplied by a bias
factor c, which varied between 0.7 and 1.1 (i.e., c = 1
assumed no bias, c < 1 assumed overestimation in the
abundance data, and c > 1 assumed underestimation).
We found that such bias affected not only the parameter
values but also the model performance in terms of model
weights. In both species, the greatest weight was
observed for the Beverton–Holt model with τ = 0.
Changes in the bias factor did not significantly influence
the weight of this model for wolves, but for lynx, the
weight of the Beverton–Holt model increased as the
assumed overestimation was increased. However, the
weights of the other models displayed an inverse trend
according to the changes in the bias factor. Therefore, the
likely overestimation in the wolf and lynx abundance
data indicated a greater correspondence with the
Beverton–Holt model, and the model selection results
were robust regardless of the existence of estimation bias.
The effect of unknown observation errors was investi-
gated by examining 5000 simulated datasets, in which a
log-normally distributed noise was added to the original
data. This treatment had almost no effect on the average
performance of the Beverton–Holt model with τ = 0 for
the wolf population as the mean value of the model
weight was similar to the original value. For the lynx
population, this was not the case as the mean value of
the model weight in the simulations increased for the
models of spatially random individual distribution and
scramble competition. As the bias correction is the oppo-
site of adding more noise, when a parameter value
decreases with an added observation error, then the true
value is on average larger than the original estimate (and
vice versa). Therefore, when the model weights decreased
with added noise, the contest competition models would
be supported even more strongly than was observed.
SˇUBA ET AL. 13
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
This study examined the consistency of wolf and lynx
population dynamics with the models based on scramble
or contest competition. The differences found in model
performance may be relevant when considering the
application of models other than the most popular ones,
such as the Ricker model, to study and predict changes in
population dynamics in single- or multi-species frame-
works. In our case, the contest competition models, such
as the Beverton–Holt model, showed better fit. Hence,
these models may be better suited for predicting carni-
vore population dynamics. The time series used in our
study, which spanned more than six decades, contained
no additional information on age structure or on the pro-
portion of males and females within the studied popula-
tions. As more detailed demographic data become
available, we believe that age- or stage-structured popula-
tion models will be more appropriate for species like
wolves and lynx. Matrix models are adaptable and can
include density dependence (e.g., Miller et al., 2002), or
multiple age or stage classes can be introduced into popu-
lation models such as the Ricker model (Shelton
et al., 2012). Whether and how such models reflect the
intraspecific competition remains an open question for
further studies.
ACKNOWLEDGMENTS
This study was conducted as part of the postdoctoral
research project “Sub-population dynamics of grey wolf
Canis lupus and Eurasian lynx Lynx lynx in Latvia and
identification of depredation risk on livestock”
(No.1.1.1.2/VIAA/3/19/511), funded by the European
Regional Development Fund (Agreement No. 1.1.1.2./
16/I/001). We thank our colleague Dr. Janis Ozolin¸sˇ for
reading the manuscript and providing valuable sugges-
tions and Ginta Sˇuba for assistance in image
preparation.
CONFLICT OF INTEREST STATEMENT
The authors declare no conflicts of interest.
ORCID
Jurģis Sˇuba https://orcid.org/0000-0001-5639-5863
Yukichika Kawata https://orcid.org/0000-0001-9511-
6863
Andreas Lindén https://orcid.org/0000-0002-5548-2671
REFERENCES
Aho, K., Derryberry, D. W., & Peterson, T. (2014). Model selection
for ecologists: The worldviews of AIC and BIC. Ecology, 95,
631–636. https://doi.org/10.1890/13-1452.1
Anazawa, M. (2009). The mechanistic basis of population models
with various types of competition (Theory of Biomathematics
and its Applications V). RIMS Koky uroku, 1663, 176–181.
Anazawa, M. (2019). Inequality in resource allocation and popula-
tion dynamics models. Royal Society Open Science, 6, 182178.
https://doi.org/10.1098/rsos.182178
Andersone-Lilley, Z., & Ozolins, J. (2005). Game mammals in
Latvia: Present status and future prospects. Scottish Forestry, 59,
13–18.
Apollonio, M., Andersen, R., & Putman, R. (Eds.). (2010). European
ungulates and their management in the 21st century. Cambridge
University Press.
Bagrade, G., Run¸ģis, D. E., Ornicans, A., Sˇuba, J., Zˇunna, A.,
Howlett, S. J., Lukins, M., Gailīte, A., Stepanova, A., Done, G.,
Gaile, A., Bitenieks, K., Mihailova, L., Baumanis, J., &
Ozolin¸sˇ, J. (2016). Status assessment of Eurasian lynx in Latvia
linking genetics and demography – A growing population or a
source-sink process? Mammal Research, 61, 337–352. https://
doi.org/10.1007/s13364-016-0279-8
Balcˇiauskas, L., & Kawata, Y. (2009). Estimation of carrying capac-
ity and growth rate of wolf in Lithuania. Acta Zoologica Litua-
nica, 19, 79–84. https://doi.org/10.2478/v10043-009-0018-3
Beverton, R. J. H., & Holt, S. J. (1957). On the dynamics of exploited
fish populations. HM Stationary Office.
Boitani, L., Alvarez, F., Anders, O., Andren, H., Avanzinelli, E.,
Balys, V., Blanco, J. C., Breitenmoser, U., Chapron, G.,
Ciucci, P., Dutsov, A., Groff, C., Huber, D., Ionescu, O.,
Knauer, F., Kojola, I., Kubala, J., Kutal, M., Linnell, J., …
Zlatanova, D. (2015). Key actions for large carnivore populations
in Europe. Report to DG Environment, European Commission,
Bruxelles. Contract no. 07.0307/2013/654446/SER/B3. Institute
of Applied Ecology.
Boukal, D. S., & Berec, L. (2002). Single-species models of the Allee
effect: Extinction boundaries, sex ratios and mate encounters.
Journal of Theoretical Biology, 218, 375–394. https://doi.org/10.
1006/jtbi.2002.3084
Brännström, Å., & Sumpter, D. J. T. (2005). The role of competition
and clustering in population dynamics. Proceedings of the Royal
Society B, 272, 2065–2072. https://doi.org/10.1098/rspb.2005.
3185
Breitenmoser, U. (1998). Large predators in the Alps: The fall and
rise of man's competitors. Biological Conservation, 83, 279–289.
https://doi.org/10.1016/S0006-3207(97)00084-0
Breitenmoser, U., Breitenmoser-Wursten, C., & Capt, S. (1998). Re-
introduction and present status of the lynx (Lynx lynx) in
Switzerland. Hystrix, 10(1), 17–30. https://doi.org/10.4404/
hystrix-10.1-4118
Breitenmoser, U., Breitenmoser-Wursten, C., Okarma, H.,
Kaphegyi, T., Kaphegyi-Wallmann, U., & Muller, U. M. (2000).
Action plan for the conservation of the Eurasian lynx in Europe
(Lynx lynx). Nature and Environment, 112, 1–69.
Carpio, A. J., Apollonio, M., & Acevedo, P. (2021). Wild ungulate
overabundance in Europe: Contexts, causes, monitoring and
management recommendations. Mammal Review, 51, 95–108.
https://doi.org/10.1111/mam.12221
Chapron, G., Kaczensky, P., Linnell, J. D. C., von Arx, M.,
Huber, D., Andrén, H., L
opez-Bao, J. V., Adamec, M.,
A´lvares, F., Anders, O., Balcˇiauskas, L., Balys, V., Bed}o, P.,
Bego, F., Blanco, J. C., Breitenmoser, U., Brøseth, H., Bufka, L.,
Bunikyte, R., … Boitani, L. (2014). Recovery of large carnivores
in Europe's modern human-dominated landscapes. Science,
346, 1517–1519. https://doi.org/10.1126/science.1257553
14 SˇUBA ET AL.
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Cook, J. R., & Stefanski, L. A. (1994). Simulation-extrapolation esti-
mation in parametric measurement error models. Journal of the
American Statistical Association, 89, 1314–1328. https://doi.org/
10.2307/2290994
Dammhahn, M., & Kappeler, P. M. (2010). Scramble or contest
competition over food in solitarily foraging mouse lemurs
(Microcebus spp.): New insights from stable isotopes. American
Journal of Physical Anthropology, 141, 181–189. https://doi.org/
10.1002/ajpa.21129
Eberhardt, L. L. (1987). Population projections from simple models.
Journal of Applied Ecology, 24, 103–118. https://doi.org/10.
2307/2403790
Eskola, H. T. M., & Geritz, S. A. H. (2007). On the mechanistic deri-
vation of various discrete-time population models. Bulletin of
Mathematical Biology, 69, 329–346. https://doi.org/10.1007/
s11538-006-9126-4
Euorpean Council. (1992). Council Directive 92/43/EEC of 21 May
1992 on the conservation of natural habitats and of wild fauna
and flora. https://eur-lex.europa.eu/legal-content/EN/TXT/?
uri=CELEX:31992L0043
Everatt, K. T., Moore, J. F., & Kerley, G. I. H. (2019). Africa's apex
predator, the lion, is limited by interference and exploitative
competition with humans. Global Ecology and Conservation, 20,
e00758. https://doi.org/10.1016/j.gecco.2019.e00758
Frank, L. G., & Woodroffe, R. (2001). Behaviour of carnivores in
exploited and controlled populations. In J. L. Gittleman, S. M.
Funk, D. W. Macdonald, & R. K. Wayne (Eds.), Carnivore con-
servation (pp. 419–442). Cambridge University Press.
Fuller, T., Mech, L. D., & Cochrane, J. F. (2003). Wolf population
dynamics. In D. L. Mech & L. Boitani (Eds.), Wolves: Behavior,
ecology, and conservation (pp. 161–191). University of Chicago
Press.
Geritz, S. A. H., & Kisdi, E. (2004). On the mechanistic
underpinning of discrete-time population models with complex
dynamics. Journal of Theoretical Biology, 228, 261–269. https://
doi.org/10.1016/j.jtbi.2004.01.003
Getz, W. M., & Owen-Smith, N. (2011). Consumer-resource dynam-
ics: Quantity, quality, and allocation. PLoS One, 6, e14539.
https://doi.org/10.1371/journal.pone.0014539
Hall, C. A. S. (1988). An assessment of several of the historically
most influential theoretical models used in ecology and of the
data provided in their support. Ecological Modelling, 43, 5–31.
https://doi.org/10.1016/0304-3800(88)90070-1
Hamada, M. Y. (2024). Dynamics of a ricker type predator–prey
model. Rendiconti del Circolo Matematico di Palermo Series 2.
https://doi.org/10.1007/s12215-024-01062-y
Hassell, M. P. (1975). Density-dependence in single-species popula-
tions. Journal of Animal Ecology, 44, 283–295.
Hayes, R. D., & Harestad, A. S. (2000). Demography of a recovering
wolf population in the Yukon. Canadian Journal of Zoology, 78,
36–48. https://doi.org/10.1139/z99-186
Herfindal, I., Linnell, J. D. C., Odden, J., Nilsen, E. B., &
Andersen, R. (2005). Prey density, environmental productivity
and home-range size in the Eurasian lynx (Lynx lynx). Journal
of Zoology (London), 265, 63–71. https://doi.org/10.1017/
S0952836904006053
Howlett, C., & Wheeler, B. C. (2021). Prenatal androgen effects as a
proximate mechanism underpinning variation in social behav-
ior among female nonhuman primates. International Journal of
Primatology, 42, 301–329. https://doi.org/10.1007/s10764-021-
00204-8
Isbell, L. A. (1991). Contest and scramble competition: Patterns of
female aggression and ranging behavior among primates.
Behavioral Ecology, 2, 143–155. https://doi.org/10.1093/beheco/
2.2.143
Janson, C. H., & Van Schaik, C. P. (1988). Recognizing the
many faces of primate food competition: Methods. Behav-
iour, 105, 165–186. https://doi.org/10.1163/156853988X0
0502
Jędrzejewska, B., & Jędrzejewski, W. (1998). Predation in vertebrate
communities. The Białowieża Primeval Forest as a case study.
Springer Verlag.
Jędrzejewska, B., Jędrzejewski, W., Bunevich, A. N.,
Minkowski, L., & Okarma, H. (1996). Population dynamics of
wolves Canis lupus in Białowieża Primeval Forest (Poland and
Belarus) in relation to hunting by humans, 1847–1993. Mam-
mal Review, 26, 103–126. https://doi.org/10.1111/j.1365-2907.
1996.tb00149.x
Johansson, A., & Sumpter, D. J. T. (2003). From local interactions
to population dynamics in site-based models of ecology. Theo-
retical Population Biology, 64, 497–517. https://doi.org/10.1016/
S0040-5809(03)00076-5
Johst, K., Berryman, A., & Lima, M. (2008). From individual inter-
actions to population dynamics: Individual resource partition-
ing simulation exposes the causes of nonlinear intra-specific
competition. Population Ecology, 50, 79–90. https://doi.org/10.
1007/s10144-007-0061-5
Kappeler, P. M., & van Schaik, C. P. (2002). Evolution of primate
social systems. International Journal of Primatology, 23, 707–
740. https://doi.org/10.1023/A:1015520830318
Kawata, Y. (2008). Estimation of carrying capacities of large carni-
vores in Latvia. Acta Zoologica Lituanica, 18, 3–9. https://doi.
org/10.2478/v10043-008-0001-4
Kelker, G. H. (1947). Computing the rate of increase for deer. The
Journal of Wildlife Management, 11, 177–183. https://doi.org/
10.2307/3795562
Knape, J. (2008). Estimability of density dependence in models of
time series data. Ecology, 89, 2994–3000. https://doi.org/10.
1890/08-0071.1
Knott, C., Beaudrot, L., Snaith, T., White, S., Tschauner, H., &
Planansky, G. (2008). Female-female competition in Bornean
orangutans. International Journal of Primatology, 29, 975–997.
https://doi.org/10.1007/s10764-008-9278-1
Krebs, C. J. (2014). Ecology: The experimental analysis of distribution
and abundance (6th ed.). Pearson Education Limited.
Lennox, R. J., Gallagher, A. J., Ritchie, E. G., & Cooke, S. J. (2018).
Evaluating the efficacy of predator removal in a conflict-prone
world. Biological Conservation, 224, 277–289. https://doi.org/10.
1016/j.biocon.2018.05.003
Linnell, J., Salvatori, V., & Boitani, L. (2008). Guidelines for popula-
tion level management plans for large carnivores in Europe.
A large carnivore initiative for Europe report prepared for
the European Commission (contract 070501/2005/
424162/MAR/B2).
Linnell, J. D. C., Brøseth, H., Solberg, E. J., & Brainerd, S. M.
(2005). The origins of the southern Scandinavian wolf popula-
tion: Potential for natural immigration in relation to dispersal
distances, geography and Baltic ice. Wildlife Biology, 11,
SˇUBA ET AL. 15
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
383–391. https://doi.org/10.2981/0909-6396(2005)11[383:
TOOTSS]2.0.CO;2
Linnell, J. D. C., Swenson, J. E., & Andersen, R. (2001). Predators
and people: Conservation of large carnivores is possible at high
human densities if management policy is favourable. Animal
Conservation, 4, 345–349. https://doi.org/10.1017/
S1367943001001408
Liu, M., & Wang, K. (2009). Survival analysis of stochastic single-
species population models in polluted environments. Ecological
Modelling, 220, 1347–1357. https://doi.org/10.1016/j.ecolmodel.
2009.03.001
L
opez-Bao, J. V., Aronsson, M., Linnell, J. D. C., Odden, J.,
Persson, J., & Andrén, H. (2019). Eurasian lynx fitness shows
little variation across Scandinavian human-dominated land-
scapes. Scientific Reports, 9, 8903. https://doi.org/10.1038/
s41598-019-45569-2
L
opez-Bao, J. V., Bruskotter, J., & Chapron, G. (2017). Finding
space for large carnivores. Nature Ecology & Evolution, 1, 0140.
https://doi.org/10.1038/s41559-017-0140
Luo, Y., Wang, L., Yang, L., Wang, X., Tan, M., & Li, Z. (2020).
Inter- and intraspecific vigilance patterns of two sympatric
Tibetan ungulates. Journal of Mammalogy, 101, 498–506.
https://doi.org/10.1093/jmammal/gyz175
May, R. M. (1974). Biological populations with nonoverlapping gen-
erations: Stable points, stable cycles, and chaos. Science, 186,
645–647. https://doi.org/10.1126/science.186.4164.645
May, R. M. (1976). Simple mathematical models with very compli-
cated dynamics. Nature, 261, 459–467. https://doi.org/10.1038/
261459a0
Mech, L. D. (1970). The wolf: The ecology and behaviour of an
endangered species. University of Minnesota Press.
Mech, L. D., Barber-Meyer, S. M., & Erb, J. (2016). Wolf (Canis
lupus) generation time and proportion of current breeding
females by age. PLoS One, 11, e0156682. https://doi.org/10.
1371/journal.pone.0156682
Mech, L. D., & Boitani, L. (2003). Wolf social ecology. In D. L.
Mech & L. Boitani (Eds.), Wolves: Behavior, ecology, and conser-
vation (pp. 1–34). University of Chicago Press.
Miller, D. H., Jensen, A. L., & Hammill, J. H. (2002). Density depen-
dent matrix model for gray wolf population projection. Ecologi-
cal Modelling, 151, 271–278. https://doi.org/10.1016/S0304-3800
(01)00493-8
Nicholson, A. J. (1954). An outline of the dynamics of animal popu-
lations. Australian Journal of Zoology, 2, 9–65. https://doi.org/
10.1071/ZO9540009
Nilsen, E. B., Linnell, J. D. C., Odden, J., Samelius, G., &
Andrén, H. (2012). Patterns of variation in reproductive param-
eters in Eurasian lynx (Lynx lynx). Acta Theriologica, 57, 217–
223. https://doi.org/10.1007/s13364-011-0066-5
Ozolin¸sˇ, J., Bagrade, G., Männil, P., & Balcˇiauskas, L. (2022). Eur-
asian lynx in Latvia: Experience gained and future challenges
at a population level. Cat News, 75, 37–41.
Ozolin¸sˇ, J., Bagrade, G., Ornicans, A., Zˇunna, A., Done, G.,
Stepanova, A., Pilate, D., Sˇuba, J., Lukins, M., & Howlett, S. J.
(2017). Action plan for Eurasian lynx Lynx lynx conservation
and management. Latvian State Forest Research Institute
Silava.
Ozolin¸sˇ, J., Pupila, A., Ornicans, A., & Bagrade, G. (2008). Lynx
management in Latvia: Population control or sport hunting? In
O. Opermanis & G. Whitelaw (Eds.), Economic, social and cul-
tural aspects in biodiversity conservation (pp. 59–72). Press of
the University of Latvia.
Ozolin¸sˇ, J., Zˇunna, A., Howlett, S. J., Bagrade, G., Pilate, D.,
Ornicans, A., & Peterhofs, E. (2016). Population dynamics of
large mammals in Latvia with an emphasis on prey-predator
interactions. In M. Stubbe (Ed.), Beiträge zur Jagd- und Wild-
forschung, Bd 41 (pp. 59–73). Gesellschaft für Wildtier- und
Jagdforschung e.V.
Ozolin¸sˇ, J., Zˇunna, A., Ornicans, A., Done, G., Stepanova, A.,
Pilate, D., Sˇuba, J., Lukins, M., Howlett, S. J., & Bagrade, G.
(2017). Action plan for Grey wolf Canis lupus conservation and
management. Latvian State Forest Research Institute Silava.
Petersson, E., & Järvi, T. (2000). Both contest and scramble compe-
tition affect the growth performance of brown trout, Salmo
Trutta, parr of wild and of sea-ranched origins. Environmental
Biology of Fishes, 59, 211–218. https://doi.org/10.1023/A:
1007645411586
Pletscher, D. H., Ream, R. R., Boyd, D. K., Fairchild, M. W., &
Kunkel, K. E. (1997). Population dynamics of a recolonizing
wolf population. Journal of Wildlife Management, 61, 459–465.
https://doi.org/10.2307/3802604
R Core Team. (2017). R: A language and environment for statistical
computing. R Foundation for Statistical Computing. https://
www.R-project.org/
Reinhardt, I., Kluth, G., Nowak, C., Szentiks, C. A., Krone, O.,
Ansorge, H., & Mueller, T. (2019). Military training areas facili-
tate the recolonization of wolves in Germany. Conservation Let-
ters, 12, e12635. https://doi.org/10.1111/conl.12635
Ricker, W. E. (1954). Stock and recruitment. Journal of the Fisheries
Board of Canada, 11, 559–623. https://doi.org/10.1139/f54-039
Rieucau, G., & Giraldeau, L.-A. (2008). Group size effect caused by
food competition in nutmeg mannikins (Lonchura punctulata).
Behavioral Ecology, 20, 421–425. https://doi.org/10.1093/
beheco/arn144
Rockwood, L. L. (2015). Introduction to population ecology (2nd
ed.). Wiley-Blackwell.
Rodenhouse, N. L., Sillett, T. S., Doran, P. J., & Holmes, R. T.
(2003). Multiple density-dependence mechanisms regulate a
migratory bird population during the breeding season. Proceed-
ings of the Royal Society B, 270, 2105–2110. https://doi.org/10.
1098/rspb.2003.2438
Royama, T. (1992). Analytical population dynamics. Chapman &
Hall.
Ruan, S. (2006). Delay differential equations in single species
dynamics. In O. Arino, M. L. Hbid, & E. A. Dads (Eds.), Delay
differential equations and applications. NATO science series
(Vol. 205, pp. 477–517). Springer. https://doi.org/10.1007/1-
4020-3647-7_11
Ruokolainen, L., Lindén, A., Kaitala, V., & Fowler, M. S. (2009).
Ecological and evolutionary dynamics under coloured environ-
mental variation. Trends in Ecology and Evolution, 24, 555–563.
https://doi.org/10.1016/j.tree.2009.04.009
Schmidt, K., Jędrzejewski, W., & Okarma, H. (1997). Spatial organi-
zation and social relations in the Eurasian lynx population in
Białowieża Primeval Forest, Poland. Acta Theriologica, 42, 289–
312. https://doi.org/10.4098/AT.arch.97-30
Shelton, A. O., Munch, S. B., Keith, D., & Mangel, M. (2012). Mater-
nal age, fecundity, egg quality, and recruitment: Linking stock
16 SˇUBA ET AL.
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
structure to recruitment using an age-structured Ricker model.
Canadian Journal of Fisheries and Aquatic Sciences, 69, 1631–
1641. https://doi.org/10.1139/F2012-082
Sibly, R. M., Barker, D., Hone, J., & Pagel, M. (2007). On the stabil-
ity of populations of mammals, birds, fish and insects. Ecology
Letters, 10, 970–976. https://doi.org/10.1111/j.1461-0248.2007.
01092.x
Skalski, J. R., Ryding, K. E., & Millspaugh, J. (2005). Wildlife demog-
raphy: Analysis of sex, age, and count data. Elsevier Academic
Presss.
Skellam, J. G. (1951). Random dispersal in theoretical populations.
Biometrika, 38, 196–218. https://doi.org/10.2307/2332328
Smith, D. W., Peterson, R. O., & Houston, D. B. (2003). Yellowstone
after wolves. Bioscience, 53, 330–340. https://doi.org/10.1641/
0006-3568(2003)053[0330:YAW]2.0.CO;2
Stefanski, L. A., & Cook, J. R. (1995). Simulation-extrapolation: The
measurement error jackknife. Journal of the American Statisti-
cal Association, 90, 1247–1256. https://doi.org/10.2307/2291515
Sˇuba, J., Kawata, Y., & Lindén, A. (2023). Properties and interpreta-
tion of the Skellam model—A discrete-time contest competition
population model. Population Ecology, 66, 115–122. https://doi.
org/10.1002/1438-390X.12169
Sˇuba, J., Zˇunna, A., Bagrade, G., Done, G., Lukins, M.,
Ornicans, A., Pilate, D., Stepanova, A., & Ozolin¸sˇ, J. (2021).
Closer to carrying capacity: Analysis of the internal demo-
graphic structure associated with the management and density
dependence of a controlled wolf population in Latvia. Sustain-
ability, 13, 9783. https://doi.org/10.3390/su13179783
Symonds, M. R. E., & Moussalli, A. (2011). A brief guide to model
selection, multimodel inference and model averaging in beha-
vioural ecology using Akaike's information criterion. Behavioral
Ecology and Sociobiology, 65, 13–21. https://doi.org/10.1007/
s00265-010-1037-6
Szewczyk, M., Nowak, S., Niedz´wiecka, N., Hulva, P., Sˇpinkytė-
Bacˇkaitienė, R., Demjanovicˇov
a, K., Cˇern
a Bolfíkov
a, B.,
Antal, V., Fenchuk, V., Figura, M., Tomczak, P., Stachyra, P.,
Stępniak, K., M., Zwijacz-Kozica, T., & Mysłajek, R. W. (2019).
Dynamic range expansion leads to establishment of a new,
genetically distinct wolf population in Central Europe. Scientific
Reports, 9, 19003. https://doi.org/10.1038/s41598-019-55273-w
van Schaik, C. P. (1989). The ecology of social relationships
amongst female primates. In V. Standen & R. A. Foley (Eds.),
Comparative Socioecology (pp. 195–218). Blackwell.
Vanags, J. (2010). Medības: Atzin¸as un Patiesības [Hunting: Opin-
ions and Truth]. Author's edition (in Latvian).
Vilà, C., Sundqvist, A.–. K., Flagstad, Ø., Seddon, J., Björnerfeldt, S.,
Kojola, I., Casulli, A., Sand, H., Wabakken, P., & Ellegren, H.
(2003). Rescue of a severely bottlenecked wolf (Canis lupus)
population by a single immigrant. Proceedings of the Royal Soci-
ety B, 270, 91–97. https://doi.org/10.1098/rspb.2002.2184
von Arx, M., Breitenmoser-Würsten, C., Zimmermann, F., &
Breitenmoser, U. (Eds.). (2004). Status and conservation of the
Eurasian lynx (Lynx lynx) in Europe in 2001. KORA Bericht
No 19.
Wagenmakers, E.-J., & Farrell, S. (2004). AIC model selection using
Akaike weights. Psychonomic Bulletin & Review, 11, 192–196.
https://doi.org/10.3758/BF03206482
SUPPORTING INFORMATION
Additional supporting information can be found online
in the Supporting Information section at the end of this
article.
How to cite this article: Sˇuba, J., Kawata, Y., &
Lindén, A. (2024). Local population dynamics of
gray wolves Canis lupus and Eurasian lynx Lynx
lynx exhibit consistency with intraspecific contest
competition models. Population Ecology, 1–18.
https://doi.org/10.1002/1438-390X.12197
SˇUBA ET AL. 17
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
APPENDIX 1: EXPRESSION OF TWO-
PARAMETER POPULATION MODELS IN TERMS
OF INTRINSIC GROWTH RATE AND CARRYING
CAPACITY
To estimate the carrying capacity and intrinsic growth
rate of wolf and lynx populations by fitting the selected
models to the data, the parameters αi and βi of the Ricker
(i = R), Skellam (i = S), Hassell (i = H), and Beverton–
Holt (i = BH) model equations (Table 1) were expressed
in terms of carrying capacity and intrinsic growth rate.
Intrinsic growth rate is the maximum population growth
rate in the absence of competition. It can be formulated
as a finite rate λ0 or an instantaneous rate r0. For species
such as the gray wolf and Eurasian lynx, which repro-
duce once per year during a breeding season, λ0 may be
preferred as it is analogous to the annual or geometric
growth rate λt = Nt+1/Nt. The analytical expression of λ0
in terms of the given model parameters was derived as:
λ0¼ lim
Nt!0
f Ntjαi,βið Þ
Nt
: ðA1Þ
Parameter K pertaining to carrying capacity was
derived from the expected equilibrium density:
f Ntjαi,βið Þ
Nt

Nt¼K
¼ 1: ðA2Þ
Solving parameters αi and βi for the Ricker, Hassel,
and Beverton–Holt models according to Equations (A1)
and (A2) provided αR = αH = αBH = λ0 and βR = ln(λ0)/
K, βH = (
ffiffiffiffi
λ0
p 
1)/K, and βBH= (λ0
 1)/K, respectively.
For the Skellam model, Equation (A1) yielded:
αS¼ λ0
βS
: ðA3Þ
Combining Equations (A2) and (A3) and conducting
algebraic rearrangements provided a nonlinear equation:
βSK
λ0ð ÞeβSK
λ0 ¼

λ0
eλ0
: ðA4Þ
The parameter βS was found by applying Lambert's
W function (also called the omega function or product
logarithm) to Equation (A4):
βS¼
λ0þW 
 λ0eλ0
 	
K
: ðA5Þ
Rearranged equations of the applied population
models in terms of λ0 and K are given in Table A1. Cer-
tain calculations are facilitated by applying instantaneous
growth rates. Ultimately, using r0 = ln(λ0) was found to
be more convenient for our study.
TABLE A1 Expression of Ricker, Skellam, Hassell, and
Beverton–Holt population models in terms of intrinsic geometric
growth rate λ0 and carrying capacity K.
Model Expression f(Ntjλ0, K)
Ricker Ntλ01

Nt
K
Skellam
K
1
exp 
 λ0þW 
 λ0
eλ0
 h i
Nt
K
 
1þλ0
1W 
 λ0
eλ0
 
Hassell Ntλ0
1þ ffiffiffiλ0p 
1½ NtKð Þ2
Beverton–Holt Ntλ0
1þ λ0
1ð ÞNtK
18 SˇUBA ET AL.
 1438390x, 0, D
ow
nloaded from
 https://esj-journals.onlinelibrary.wiley.com/doi/10.1002/1438-390X.12197 by Luonnonvarakeskus, W
iley Online Library on [10/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on W
iley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License