Spatial distribution and optimal harvesting of an age-structured population in a fluctuating environment

Mathematical Biosciences 296 (2018) 36–44

Contents lists available at ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Spatial distribution and optimal harvesting of an age-structured population in a fluctuating environment

T

⁎

Steinar Engen ,a, Aline Magdalena Leeb, Bernt-Erik Sætherb a b

Centre for Biodiversity Dynamics, Department of Mathematical Sciences, Norwegian University of Science and Technology, Trondheim N-7491, Norway Centre for Biodiversity Dynamics, Department of Biology, Norwegian University of Science and Technology, Trondheim N-7491, Norway

A R T I C L E I N F O

A B S T R A C T

Keywords: Age-specific harvest model Environmental stochasticity Optimal harvesting Spatial synchrony Population fluctuations Stochastic population dynamics

We analyze a spatial age-structured model with density regulation, age specific dispersal, stochasticity in vital rates and proportional harvesting. We include two age classes, juveniles and adults, where juveniles are subject to logistic density dependence. There are environmental stochastic effects with arbitrary spatial scales on all birth and death rates, and individuals of both age classes are subject to density independent dispersal with given rates and specified distributions of dispersal distances. We show how to simulate the joint density fields of the age classes and derive results for the spatial scales of all spatial autocovariance functions for densities. A general result is that the squared scale has an additive term equal to the squared scale of the environmental noise, corresponding to the Moran effect, as well as additive terms proportional to the dispersal rate and variance of dispersal distance for the age classes and approximately inversely proportional to the strength of density regulation. We show that the optimal harvesting strategy in the deterministic case is to harvest only juveniles when their relative value (e.g. financial) is large, and otherwise only adults. With increasing environmental stochasticity there is an interval of increasing length of values of juveniles relative to adults where both age classes should be harvested. Harvesting generally tends to increase all spatial scales of the autocovariances of densities.

1. Introduction Predicting consequences of harvesting natural populations is challenging because exploitation is likely to interfere with natural processes that affect patterns of variation in population size [2,17,21,23,39]. Such harvest-induced modifications of the population dynamics may strongly affect the sustainability of long-term harvesting. An important achievement from analyses of population models has been the derivation of general principles for choice of harvest strategies that can act as guidelines for management of harvested populations. For example, a density-regulated population in a stable environment should according to the principle of maximum sustainable yield be maintained at a population size that maximizes the population growth [9]. Similarly, using a simple age-structured model, Beverton and Holt [5] introduced the yield per recruit as a key parameter, which can be used to find the optimal fishing rate that maximizes the annual yield. More complex harvesting models with interactions among two [8] or more species [29], and possibly harvest of more than one species have also been analyzed. These initial models ignored two important characteristics of most exploited populations. First, the analyses were based on purely deterministic models. However, the dynamics of most natural

⁎

populations are strongly affected by random variation in the environment [37], which requires models with environmental stochasticity. Such models generally imply modifications toward more conservative harvesting. Large mean annual yield is then obtained by different types of threshold strategies that reduce the probability of an unexpected population crash or even extinction [24,25,36]. Second, the early harvest models all assumed that individuals are homogeneously distributed in space. This contradicts recently accumulated evidence that most harvested populations show distinct spatial structures, strongly affected by the choice of harvest strategy [7,15,16,22]. In practice, this often means that individuals far apart in space may experience quite different stochastic environmental effects (e.g. [18]), but may still interact through dispersal. Distribution of individuals in space is commonly described through the spatial autocovariance function for population density [6,13,20], that is, the temporal covariance of densities at two locations expressed as a function of the distance between them. Lande et al. [26] scaled the autocovariance function to become a distribution and used its standard deviation as a measure of spatial scale and derived the remarkably simple formula l 2 = le2 + mlg2/ γ for small fluctuations. Here l, le and lg are the scales of the density, the autocovariance of environmental

Corresponding author. E-mail addresses: [email protected] (S. Engen), [email protected] (A.M. Lee), [email protected] (B.-E. Sæther).

https://doi.org/10.1016/j.mbs.2017.12.003 Received 22 June 2017; Received in revised form 20 November 2017; Accepted 9 December 2017 Available online 11 December 2017 0025-5564/ © 2017 Published by Elsevier Inc.

Mathematical Biosciences 296 (2018) 36–44

S. Engen et al.

fluctuations, and the dispersal distance, while m is dispersal rate and γ is the strength of local density regulation. If there is no dispersal, then l = le , known as the Moran effect [32]. An interesting consequence of the above result is that the effect of small dispersal rates on the scale of population densities may be large if the local density regulation is weak (small γ). Here we analyze a simplified age-structured model in continuous space and time with two age classes (or stages), juveniles and adults. To simplify the analyses, we assume no production of offspring by juveniles and that density dependence acts only through the juvenile survival rate, dependent on the density of juveniles, whereas no density regulation occurs in the adult segment. Individuals of the two classes may have different dispersal rates and distributions of dispersal distances. We derive some general spatial scaling results for this model and also study optimal proportional harvesting strategies allowing for different age-specific harvest rates.

dn1 (x ) = n2 (x ) βdt − n1 (x ) dt − n1 (x )[μ + νn1 (x )] dt − m1 n1 (x ) dt + m1 dt

∫ n1 (x − u) g1 (u) du + n2 (x ) σβ dBβ (x )

− n1 (x ) σμ dBμ (x ).

This model formulation implies that there are large population densities so that the effects of local demographic noise [14] can be ignored. The variances σβ2, σδ2 and σμ2 describe noise in rates affecting all individuals in a class and are hence environmental variance components. The noise components can not be chosen freely because any noise covariance matrix derived from them has to be positive definite. In Appendix A we show how the three spatial noise terms can be expressed in a general way using three independent environmental variables to express dBβ(x), dBμ(x), and dBδ(x) in a logically consistent way. We analyze this model in the case that the fields n1(x) and n2(x) fluctuate in a stationary way in space near the deterministic equilibrium values K1 and K2 obtained by ignoring the noise terms, obeying the equations βK2 − K1 − K1 (μ + νK1) = 0 and K1 − δK2 = 0 with solution

2. Model We consider an organism with a juvenile age class with density n1(x) at location x = (x1, x2) in space that does not produce offspring, and an adult class with density n2(x). The juveniles are subject to density-dependent mortality within the class as well as instantaneous densityindependent dispersal at rate m1 with a given two-dimensional distribution of dispersal distance g1(u). The model may be initiated as a discrete time model expressing changes during a time step (usually one year) which under moderate fluctuations can be approximated by a continuous time diffusion model, or we may formulate the model directly in continuous time. We assume that individuals stay in the juvenile state for on average one time unit before entering the class of adults producing offspring. The adults disperse continuously at a different rate m2 with distribution of distance g2(u), but are not subject to density dependence. The temporal change in n1(x) during time dt due to individuals changing age class is then − n1 (x ) dt , and the same amount is added to n2(x). In a large population this is approximately equivalent to assuming that each individual remains juvenile for exactly one time unit. Adults at location x have a temporally fluctuating death rate δ + σδ dBδ (x )/ dt . In general we write dB(x) with appropriate subscript for the temporal increment during time dt of a standard Brownian motion so that E[dB (x )] = 0 and var[dB (x )] = dt . The discrete time analogy would be to express the noise operating during a time step through a substitution of dB(x) with a stochastic variable ΔB(x) with zero mean and unit variance, while replacing dt with Δt = 1. The continuous increments are possibly correlated in space expressed by E[dB (x ) dB (x + y )] = ρ (y ) dt (with appropriate subscript), depending on the distance only. Hence, the average lifetime after reaching the adult state is 1/δ in the average environment. With these assumptions the dynamics of the adult population is given by

K1 = r / ν ,

K2 = K1/ δ = r /(δν ),

− n2 (x ) σδ dBδ (x ).

(3)

where r = R 0 − 1 − μ, and R 0 = β / δ . Since δ is the average lifetime of adults measured from the time individuals reach the adult stage and β is their rate of reproduction, R0 is the mean lifetime reproductive success of individuals that do not die as juveniles. Here (1 + μ) represents the removal rate of juveniles, that is, death rate plus transmission rate to the adult state, at small densities (n1(x) ≈ 0), while since K2 = K1/ δ , the component of the rate of increase in n1(x) at equilibrium due to births is R0. Hence, both carrying capacities are positive if this rate of increase is larger than the rate of removal of individuals at small densities, that is R 0 > 1 + μ, or r > 0. The parameter r is the natural measure of the deterministic population growth rate at small densities, and since the density regulation is of the logistic type this is also the natural measure of the strength of density regulation. 3. Linearization by small noise approximation First, assuming small fluctuations we approximate the ni(x), i = 1, 2, in the noise terms by Ki as in Lande et al. [26]. Second, we introduce the relative deviations from the deterministic equilibrium as z i (x ) = ni (x )/ Ki − 1 that are zero at the deterministic equilibrium. The dynamic equations linearized at the carrying capacities are then

dz1 (x ) = [−α11 z1 (x ) + α12 z2 (x )] dt + m1 dt

∫ z1 (x − u) g1 (u) du

+ σ1 dB1 (x ) dz2 (x ) = [α21 z1 (x ) − α22 z2 (x )] dt + m2 dt

dn2 (x ) = n1 (x ) dt − δn2 (x ) dt − m2 n2 (x ) dt + m2 dt

(2)

∫ n2 (x − u) g2 (u) du

+ σ2 dB2 (x ),

(1)

(4)

∫ z2 (x − u) g2 (u) du (5)

where α11 = 1 + μ + m1 + 2r , α12 = β / δ , α21 = δ, α22 = δ + m2 , σ2 = σδ and and σ12 = (σβ / δ )2 + σμ2 − 2ρβ, μ (0) σβ σμ/ δ dB2 (x ) = −dBδ (x ), dB1 (x ) = [(σβ / δ ) dBβ (x ) − σμ dBμ (x )]/ σ1. Engen [11] has shown, in a spatial model with no age structure, that linearization yields surprisingly accurate results for autocovariances and their spatial scales under moderate and even under fairly large fluctuations in population density. Notice that the dynamics of z1(x) and z2(x) do not depend on the parameter ν used to define the density regulation. Actually, ν is simply a scaling parameter for K1 and K2. Hence, we can later choose ν = 1 without loss of generality. To avoid too many complex parameter combinations describing the noise we shall in our examples consider a simplified model for the noise assuming that all spatial correlations among noise components are proportional to a common environmental correlation ρe(y). Then,

Adult individuals produce offspring at rate β + σβ dBβ (x )/ dt giving an input [βdt + σβ dBβ (x )] n2 (x ) to the juvenile class during dt. Furthermore, juveniles are subject to density-dependent survival by a death rate μ + σμ dBμ (x )/ dt + νn1 (x ) increasing linearly with juvenile density in analogy to standard logistic dynamics of a population with no age structure. Alternatively this logistic formulation of the density dependence may be considered as a first order approximation around the carrying capacity of other types of density dependence, such as the Beverton–Holt [5] or Ricker model [34] commonly used in analyses of fish dynamics. However, these two models express the density regulation through the adult population size, which is approximately proportional to the density of juveniles used in the present model. The dynamics of n1(x) are then given by 37

Mathematical Biosciences 296 (2018) 36–44

S. Engen et al.

δ + h2 μ + h1 replaced by and δ by so that K2 (h1, h2) = K1 (h1, h2)/(δ + h2) . Hence, introducing H2 = h2/(δ + h2) which gives β /(δ + h2) = β (1 − H2) and h2 K2 = K1 h2/(δ + h2) = K1 H2, and using the expressions for K1 and K2 under harvesting then yields

ρβ (y ) = ρδ (y ) = ρμ (y ) = ρe (y ), ρβδ (y ) = ρβδ (0) ρe (y ), ρβμ while (y ) = ρβμ (0) ρe (y ) and ρδμ (y ) = ρδμ (0) ρe (y ) . These assumptions lead to E[dB2 (x ) dB2 (x + y )]/ dt = E[dB1 (x ) dB1 (x + y )]/ dt = ρe (y ) ρ12 and (y ) = E[dB1 (x ) dB2 (x + y )] = [−ρβδ (0) σβ / δ + σμ ρδμ (0)] ρe (y )/σ1. In Appendix A we give a full description of the general case in which the three noise terms for the birth rates and the two death rates have a general environmental noise. In Appendix B we derive a set of balance equations characterizing the stationary two-dimensional spatial process z1(x), z2(x) by the spatial auto-covariance functions cij (y ) = cov[z i (x ), z j (x + y )]. In Appendix C these equations are solved by Fourier transforms. These can be transformed back to the actual covariances, transformations that are simple one dimensional integrations in the case of isotropic models, that is, when all spatial functions depend only on Euclidian distance so that all directions in space appear mathematically equivalent (Eq. (C1) in Appendix C).

νY = (Ph1 + H2)[R 0 (1 − H2) − (1 + μ + h1)]. Searching for solutions h1, H2 ≥ 0, we first find the maximum under the constraint h1 = cH2, c ranging from zero to infinity, giving the equation

νY = (cP + 1) H2 [r − H2 (R 0 + c )], from which we find that the maximum of Y with respect to H2 for a given c is proportional to (cP + 1)/(c + R 0) . Inspection of this function reveals that the maximum is attained for c = 0 if P < 1/R0 and c = ∞ if P > 1/ R 0 = δ / β (and is independent of c if P = 1/ R 0 ). Hence, the optimal solution is always to harvest only one age class. If the first class has low value, P < 1/R0, then the solution is h1 = 0 and

4. Spatial scale of autocovariances

h2 =

We define the spatial scale of covariances as proposed by Lande et al. [26], first scaling the autocovariance functions to become distributions and then use the standard deviation of these two-dimensional distributions in a given direction as a measurement of spatial scale in that direction. Here we only show examples of isotropic models, so the scales are the same in any direction. In Appendix D we derive the spatial scale lij of the covariance functions cij(y) for the simplified noise model given above. Writing le for the spatial scale of the common environmental correlation factor ρe(y), and lgi for the scale of the distribution of the dispersal distance of individuals in age-class i, the squared scales for the covariances take the form

lij2 = le2 + Gij1 m1 lg21 + Gij2 m2 lg22,

β − (1 + μ) δ rδ = . β / δ + (1 + μ) R 0 + (1 + μ)

If the value of the first class is higher, P > δ/β, then only the first class should be harvested, h2 = 0, and

h1 =

1 [R 0 − (1 + μ)] = r /2. 2

It appears that the harvesting rate is half the deterministic population growth rate for small population sizes in the absence of harvesting and that harvesting then reduces the value of n1 by 50%, which is a well known result for the simple deterministic logistic model with no age structure. More interestingly, however, when the solution is to harvest the adults only, then it appears from the above solution that this should also be done in such a way that the juvenile class and the deterministic growth rate are reduced to 50% of their values under no harvesting. To see this we simply replace δ by δ + h2 = 2β /(R 0 + 1 + μ) giving the new growth rate under harvesting as

(6)

where the Gijk are rather complex functions of the population parameters, derived in Appendix D. However, some numerical examples given below indicate that these factors are approximately inversely proportional to the natural measure r = β / δ − 1 − μ of strength of density regulation in the present age-structured model for moderate values of r (0 < r < 0.2), in analogy with the scaling result in the simple model of Lande et al. [26].

β (R 0 + 1 + μ) β − (1 + μ) = − (1 + μ) = r /2. 2β δ + h2 5.2. Solution for the stochastic model For given parameter values we have shown in Appendix C how to compute the Fourier transform of c11 and transform it back to reveal the full autocovariance function c11(y) and in particular var(z1) = c11 (0) and 2 accordingly var(n1) = K12 var(z1) = τ11 . Taking the expectations of the dynamic Eqs. (1) and (2) then yields En1 = (δ + h2)En2 and, using En12 = (En1)2 + τ11 then yields

5. Optimal proportional harvesting Proportional harvesting may easily be incorporated in this model. Suppose the two age classes are harvested in space proportional to their density at each location, that is at rates h1n1(x) and h2n2(x). This does not change the structure of the model, only the density-independent death rates μ and δ changing to μ + h1 and δ + h2 . Then, in the deterministic model, the yields per time unit are K1h1 and K2h2 for the two age classes, where K1 and K2 now are functions of h1 and h2. If the ‘value’ (for example price or average weight) of juveniles relative to the adults is P, there is an optimal solution for h1 and h2 for maximizing the total value PK1 h1 + K2 h2 . In a stochastic environment, however, one often attempts to maximize the mean annual yield [24,25] corresponding to the mean value Ph1 En1 + h2 En2 . In general, the expectations will depend on the stochastic fluctuations, as will be discussed below.

β En2 − (1 + μ + h1)En1 − ν (En1)2 − ντ11 = 0. Inserting En2 = En1/(δ + h2) and using the fact β /(δ + h2) − (1 + μ + h1) = νK1 (h1, h2) the solution for En1 is

1 En1 = K1 ⎡1 − (1 − 2 ⎣

that

2 1 − 4τ11 / K12 ) ⎤ ≈ K1 [1 − var(n1)/(En1)2], ⎦

where the last approximation is valid for small squared coefficients of variance var(n1)/(En1)2 of n1. Here, the parameters K1 and τ11 are functions of the harvesting rates h1 and h2 and En2 = En1/(δ + h2) . The optimal strategy maximizing mean annual yield when harvested individuals from the first class have value P compared to those in the second class, is then the values of h1 and h2 maximizing Ph1 En1 + h2 En2 = [Ph1 + h2/(δ + h1)]En1. The optimal harvesting rates have to be computed numerically. Some numerical illustrations will show that still h1 = 0 for small values of P and correspondingly h2 = 0 for large values. However, there is now also an intermediate interval of P for which both age classes should be harvested. The width of this interval increases and the optimal harvesting

5.1. Deterministic solution Consider first the deterministic model which under proportional harvesting has no variation of density in space. Our goal is then to find the (h1, h2 ) that maximizes the annual yield Y = PK1 (h1, h2) h1 + K2 (h1, h2) h2, where the expressions for K1 and K2 are the same as those without harvesting given by Eq. (3), but with μ 38

Mathematical Biosciences 296 (2018) 36–44

S. Engen et al.

1.0

rates decrease with increasing environmental noise.

l11 (strong)

0.9

6. Joint simulation of the density fields n1(x) and n2(x)

Autocorrelations

tion ula reg sity den ng stro

0.8

The fields n1(x) and n2(x) can be simulated using two independent spatial white noise processes with no spatial correlations, dA1(x) and dA2(x) with EdA1 (x ) = EdA2 (x ) = 0 and EdA1 (x )2 = EdA2 (x )2 = dx . Then we may express the juvenile field as z1 (x ) = ∫ ξ11 (u) dA1 (x − u) ξ11 (u) = ξ11 (−u), with giving c11 (y ) = cov[z1 (x ), z1 (x + y )] = ∫ ξ11 (u) ξ11 (y − u) du . Similarly, we express the adult density as z2 (x ) = ∫ ξ22 (u) dA1 (x − u) + ∫ ξ12 (u) dA2 (x − u), giving c22 (y ) = c12 (y ) = ∫ and ∫ ξ22 (u) ξ22 (y − u) du + ∫ ξ12 (u) ξ12 (y − u) du, ξ11 (u) ξ22 (y − u) du . We show in Appendix E how to compute the functions ξij(u). In practice, we will have to replace the dx by small squares Δx and simulate the corresponding ΔAi(x) as normal variates with zero mean and variances Δx. The above integrals are finally approximated by sums over midpoints u of all squares such as ∑ ξ11 (u)ΔA1 (x − u) giving the value of the fields in the middle of each square. Notice that the normal variates ΔA(x) do not have to be stored. For each square we can simply simulate ΔA1(x) and add its contribution ∑ ΔA1 (x ) ξ11 (x − u) to the entire field and repeat this over all squares contributing to the region we want to plot.

0.7 0.6 0.5 0.4 0.3

l11 (weak)

wea k de nsit y re gula tion

0.2 0.1 0.0

0

100

200

300

400

500

Distance Fig. 2. Autocorrelation functions c11(y)/c11(0) (solid lines), c22(y)/c22(0) (dashed lines) and c12(y)/c12(0) (dotted lines) against distance y. The lower lines marked ‘strong density regulation’ are the functions computed using the parameters in Fig. 1 where the strength of density regulation is β / δ − (μ + 1) = 0.2 . The upper lines are computed using β = 0.37 giving the much weaker density regulation of 0.05. The vertical lines show the theoretically computed spatial scales for c11(y)/c11(0) for the strong and weak density regulation. All autocorrelations describing the noise have Gaussian form.

7. Illustrative examples

for juveniles and adults for the fields shown in Fig. 1 labeled as ’strong density regulation’, together with a graph where the density regulation is weaker by the choice of a smaller birth rate β giving a smaller growth rate r when densities are small. These are computed by first computing the Fourier transforms and finally using Eq. (C1) in Appendix C to perform the backward transformation. In the weaker case, inspection of the graphs indicates considerably larger spatial scale compared to the population with stronger density regulation. This is also expressed by the vertical lines at the scales l11 computed by Eq. (6) for both cases. The fact that the scales increase with decreasing strength of density regulation is further illustrated by Fig. 3 showing the factors Gij1 and Gij2 in Eq. (6) multiplied by the strength r of density regulation, plotted against r. Over this region of r (0 < r < 0.2) these factors are nearly constant, so that the second term in Eq. (6) expressing the effect of

Fig. 1 shows a typical illustration of the spatial fields of density for both age classes simulated by using the technique described in Appendix E. The density of adults is about 5 times as high as that of juveniles in this example where the expected lifetime of adults after having reached the adult state in the corresponding average environment is on average 1/ δ = 5. It appears from the graphs that the two fields tend to follow each other due to a positive spatial autocorrelation between the two age classes. Fig. 2 provides autocorrelation functions

Fig. 3. Examples of factors Gij1 and Gij2 determining the effect of migration capacities m1 lg21 and m2 lg22 of age-class 1 and 2 on the spatial scales lij as expressed by Eq. (6). In the simple model of Lande et al. [26] without any age structure this factor is inversely proportional to the strength of density regulation. In the present model the strength of density regulation can be expressed by r = β / δ − 1 − μ, so the graphs show rGij1 and rGij2 plotted as functions of r. These would be constants if the factors were inversely proportional to the strength of density regulation. Parameter values are as in Fig. 1 except for β which is chosen as a function of r, that is, β = (1 + μ) δ + rδ .

Fig. 1. (A) Simulation of the density fields n2(x1, x2) for adults and (B) n1(x1, x2) for juveniles using the technique described in Appendix E. Parameter values are β = 0.4, σβ = 0.1, δ = 0.2, σδ = 0.04, μ = 0.8, σμ = 0.1, ν = 1, m1 = 2, lg1 = 50, m2 = 1, lg2 = 15,

l e = 50, ρβ, δ (0) = ρβ, μ (0) = −ρδ , μ (0) = −0.5 . 39

Mathematical Biosciences 296 (2018) 36–44

S. Engen et al.

Fig. 4. Optimal harvesting rates h1 for juveniles (upper panel) and h2 for adults (lower panel) as a function of the value P of juveniles relative to adults. In the example with σβ = 0.1 the other standard deviations are σδ = 0.06 and σmu = 0.04 . These standard deviations are changed proportionally to produce the other graphs. The other parameters are β = 0.8, δ = 0.2, μ = 0.8, ν = 1, m1 = 2, m2 = 1, lg1 = lg 2 = 15, l e = 50, ρβμ (0) = −0.5,

ρβδ (0) = −0.5 and ρδμ = 0.5 . All autocorrelations describing the noise have Gaussian form.

juvenile dispersal on the scale of juvenile and adult density is approximately inversely proportional to the strength of density regulation. The effect of adult dispersal given by the third term of Eq. (6) is also a considerable increase by decreasing r, although the rGij2 change up to about 30% in this interval. Altogether, this indicates that the component of the scale of densities generated by juvenile as well as adult dispersal is roughly proportional to 1/r, and all scales become very large for very small values of r (weak density regulation). Fig. 4 exemplifies computations of the optimal harvesting rates h1 (upper panel) and h2 (lower panel) as functions of the relative value P of juveniles compared to adults. For small values of P only adults should be harvested, while only juveniles should be harvested for large values of P. As we have shown, in the deterministic case the switching from one strategy to the other is discrete happening at P = δ / β = 1/ R 0 , where R0 is the lifetime reproductive success for individuals reaching the adult state. The optimal harvesting rates also have practically discrete jumps when the stochasticity is small, but for increasing stochasticity there is an interval of P of increasing length where both age classes should be harvested. This effect of increasing stochasticity would have to be investigated numerically in each case. In Fig. 4 there is assumed to be dispersal of both age classes, but as the magnitude of fluctuations in densities increases with decreasing dispersal we should expect that both age classes should be harvested in some interval for P also for smaller

Fig. 5. The carrying capacities (A) and mean annual yield (B) for the two age classes as functions of harvesting rates when harvesting the first age class only (dashed lines) as well as only adults (solid lines). Panel (C) shows the spatial scales l11 and l22 for the stochastic fluctuation in densities of the two age class under the same harvesting regimes. The other parameters are as in Fig. 4 and all autocorrelations describing the noise have Gaussian form.

(and zero) dispersal rates. Finally, increasing harvesting rates reduced the carrying capacities all the way to zero for aggressive harvesting (Fig. 5A). Fig. 5B shows the mean annual yield for harvesting the first or the second age class only with the same parameters as in Fig. 5A. As the harvesting rates increase 40

Mathematical Biosciences 296 (2018) 36–44

S. Engen et al.

Simulations of the spatial fields for both age classes using the technique derived in Appendix E (Fig. 1) illustrate the general effect that the age classes tend to end up with rather similar density patterns. This is a result of two effects. First, when juveniles pass over to the adult state they tend to bring the spatial pattern with them. Second, it is realistic to assume that there is a positive environmental correlation between the death rates of the two classes, and that these are negatively correlated with the birth rates. This is because a good environment will induce small death rates and large birth rates. In our simulations these correlations at zero distance in space are chosen as 0.5 and −0.5, respectively. The similarities in spatial patterns also appear from the spatial autocorrelations in densities (Fig. 2). Our result on spatial scale of autocovariances of densities given by Eq. (6) for our simplified noise model is a combination of transparent and less transparent results. The squared spatial scales bear striking similarities to the simple scaling formula of Lande et al. [26]. First, the squared scales all have an additive component le2 which is the squared scale of the environmental noise. With no migration this is the only term and expresses the so-called Moran effect in linear models [12,32]. Furthermore, dispersal of each age class generates additive terms proportional to mi lgi2 as in [26]. The additional constant factors in these terms are, however, rather complicated expressions. In the model of Lande et al. [26] this constant is exactly the inverse of the strength of density regulation, which is our background for looking for similar effects in our model. A few numerical examples indicate the general pattern that the factor increases with decreasing strength of density regulation which we have defined as r = R 0 − 1 − μ (Fig. 3). In most examples the factors are nearly proportional to 1/r. Thus, proportional harvesting generally reduces r and therefore induces an increase in all three spatial scales. This is demonstrated in Fig. 5, showing that the scales actually tend to infinity as the harvesting brings the carrying capacities and r toward zero. If there is no stochasticity, the density fields are constant in space and the model for proportional harvesting of both age classes is somewhat trivial. The optimal harvesting strategy in this case with our simplified age structure is always to harvest only one age class, in agreement with previous results for deterministic models with more complex structure [4,19,28]. The same solution was also reported by Sæther et al. [36] using a stochastic model with no spatial structure. When juveniles have high value relative to adults, they should be harvested at rate h1 = r/2, exactly as in a logistic model with no age structure. Interestingly, when only adults should be harvested, then h2 should also be chosen so that the growth rate of the population at small densities under harvesting is reduced to r/2. In the stochastic model the optimal solution changes qualitatively with increasing stochasticity, generating an interval of P in which both age classes should be harvested. The length of this interval increases with increasing environmental stochasticity (Fig. 4), illustrating the importance of including the effects of fluctuating environments in harvest models [3,25]. Comparative studies have revealed that harvesting strongly alters the dynamical characteristics of exploited populations [2,17,21,35]. For example, in marine systems fished populations fluctuate more in size than unharvested stocks, most likely caused by increased fishinginduced variability of key demographic parameters such as the population growth rate at small population sizes [2]. Our results show that choice of age-specific harvest strategies may also affect the spatial distribution of abundance in exploited populations (Fig. 4B). Accordingly, empirical evidence strongly indicates that harvest may affect the spatial dynamics of populations [7,16,22].

and carrying capacities decrease, the spatial scales of densities increase and approach infinity as the carrying capacities approach zero (Fig. 5C). The effects of harvesting on the spatial scaling was closely related to choice of age-specific harvest tactic (Fig. 5C). 8. Discussion The complexity of ecological systems makes them challenging to model. Their behavior is therefore generally studied through simplified models highlighting different aspects of the system. For example, agestructured models are usually analyzed ignoring density regulation. Because each class may be affected by the density of their own or other age classes in a variety of ways, models with density regulation become complicated and require many parameters [27]. Similarly, spatial agestructured models are likely to be complicated because dispersal patterns, vital rates and density regulation may all differ among age classes. For models to be realistic, they must also include stochastic effects. This has mainly been done only in models lacking both density dependence and spatial effects, because complex stochasticity in all vital rates can then be reduced to a single environmental [38] and a single demographic variance [12,14]. Although there is a large literature on harvesting of age structured populations [33], analytical results are mainly restricted to deterministic models with no density regulation. A common alternative to deriving analytical results is to perform individual-based simulations which can be done with virtually any degree of complexity. Such models may give considerable insight by studying variation in a few parameters, as exemplified by the spatial age-structured model of McCauley et al. [31], including predation and foraging. However, when including many topics in the same individualbased model, the number of parameters becomes very large, in particular when stochasticity is included. This makes it impracticable to simulate all parameter combinations and gain the overview needed to finally understand the dynamics of the system. In the present model we include age, density regulation, space, harvesting and stochasticity. With such complexity we are forced to make a number of simplifications to find some transparent analytical results. Although our model in principle could be worked out in the same way using many age classes, we have chosen only to include a juvenile and an adult stage, which is a realistic description of many organisms and considerably limits the number of parameters required. We present a rather general description of the stochasticity in this system, but assume that densities are large enough to ignore effects of demographic noise [14]. In the examples, however, we simplify the environmental stochasticity by assuming that all spatial autocorrelations in noise terms are proportional to a common environmental autocorrelation function. This reduces the number of parameters required to describe the stochastic effects. Density dependence is simplified by only including density dependence in the death rate of juveniles. Dispersal is considered as a continuous process with rates proportional to density, which is realistic at least for some species [30]. We only consider proportional harvesting, that is, harvesting rates for both age classes at any time and location are proportional to density, which may be close to using a constant harvesting effort. We also make some mathematical simplifications. First, we assume that fluctuations in population size are so small that the linearized models around the carrying capacities for both age classes provide a reasonable approximation. However, the results of Engen [11] indicate that results derived from a linearized model are remarkably accurate for moderate fluctuations and give reasonable results even for fairly large fluctuations. Second, we consider a population spread out in the entire plane to avoid dealing with complicated boundary effects, and in all examples we assume that the model is isotropic, meaning that a rotation of the space yields exactly the same model so that all directions appear equivalent.

Acknowledgments We are grateful for financial support from the Research Council of Norway (project no. SFF-III 223257/F50 and SUSTAIN).

41

Mathematical Biosciences 296 (2018) 36–44

S. Engen et al.

Appendix A. Construction of noise terms The noise terms in Eqs. (3) and (4) are given by σ1 dB1 (x ) = (σβ / δ ) dBβ (x ) − σμ dBμ (x ), so that σ12 = (σβ / δ )2 + σμ2 − 2ρβμ (0) σβ σμ/ δ , where ρβμ (y ) dt = E[dBβ (x ) dBμ (x + y )] expresses the spatial correlation between the noise in birth rates β and death rates μ of juveniles at small densities. The spatial correlations of the noise terms in Eqs. (4) and (5) are defined by the functions ρij (y ) = E[dBi (x ) dBj (x + y )]/ dt , which under the above assumptions are

ρ11 (y ) = [(σβ / δ )2ρβ (y ) + σμ2 ρμ (y ) − 2(σβ σμ/ δ ) ρβμ (y )]/ σ12,

(A1)

ρ22 (y ) = ρδ (y ),

(A2)

ρ12 (y ) = [σμ ρδμ (y ) − (σβ / δ ) ρβ, δ (y )]/(σ1).

(A3)

The noise terms dBβ(x), dBμ(x) and dBδ(x) are all temporal increments of standard Brownian motions with a spatial correlation within as well as between rates. The terms in Eqs. (4) and (5) constructed by them are environmental noise terms describing fluctuations in birth- and death-rates. A simple construction is to start with three independent environmentally determined noise increments dAi(x), i=1,2,3, with specified legal (positive definite) spatial autocorrelations ρi (y ) = E[dAi (x ) dAi (x + y )]/ dt , for example chosen to have a Gaussian form. We then define the noise in the vital rates as linear combinations such as 3

dBβ (x ) =

3

∑ ui dAi (x )/ ∑i =1 ui2 , i=1

and similarly for dBμ and dBδ with the ui replaced by other weights, say vi and wi. Then, all increments have expectation zero and variance dt as required, and all spatial covariances can be expressed easily by the ρi(y). Appendix B. Balance equations for autocovariance functions Because Eqs. (4) and (5) are linear the solution is two correlated Gaussian fields z1(x) and z2(x). Our main goal is then to derive the spatial autocovariance functions cij (y ) = cov[z i (x ), z j (x + y )], i, j = 1, 2, which then represent a complete characterization of the solution. This can be achieved by writing up three balance equations that must be obeyed under stationarity, expressing that spatial covariances do not change through time. These are cov[z i (x ), z j (x + y )] = cov[z i (x ) + dz i (x ), z j (x + y ) + dz j (x + y )], where the increments dzi(x) are given by Eqs. (4) and (5). This gives for example

cov[z1 (x ), dz1 (x + y )] +

1 cov[dz1 (x ), dz2 (x + y )] = 0. 2

Inserting the dz1 from Eq. (4) and dividing through by dt then gives the balance equation

− α11 c11 (y ) + α12 c12 (y ) + m1

∫ c11 (y − u) g1 (u) du + 12 σ12 ρ11 (y) = 0.

(B1)

Similarly, we find the two other balance equations,

α21 c12 (y ) − α22 c22 (y ) + m2

∫ c22 (y − u) g2 (u) du + 12 σ22 ρ22 (y) = 0,

0 = − (α11 + α22) c12 (y ) + α12 c22 (y ) + m1 + α21 c11 (y ) + m2

(B2)

∫ c12 (y − u) g1 (u) du

∫ c12 (y) g2 (u) du + σ1 σ2 ρ12 (y).

(B3)

Appendix C. Solution by Fourier transforms The Fourier transform of a spatial function c(y) is f (ω) =

∞ ∞ ∫−∞ ∫−∞ ei (ω1y1 + ω2 y2) c (y ) dy1 dy2 where i = −1 . For an isotropic model where

c (y ) = c͠ (r ), where r = y12 + y22 , the Fourier transform can be expressed as f (ω) = f͠ (u), where u = given by an integration in one dimension c͠ (r ) =

1 2π

∫0

∞

ω12 + ω22 . Then the backward transformation is

f͠ (u) J0 (ru) udu,

(C1)

where J0 is the Bessel function of the first kind of order zero [1]. If c(y) has an isotropic Gaussian form, 2

c (y ) = c0 e (−y1 / σ

2 + y 2 / σ 2)/2 2

= c0 e−r

2 /(2σ 2)

2 2 the Fourier transform is f͠ (u) = 2πσ 2c0 e−u σ /2 . Writing fij(ω) for the Fourier transform of cij(y) then gives, from Eqs. (B1)–(B3), three linear equations in f11, f22, and f12,

− f11 (α11 − m1 fg1 ) + α12 f12 +

α21 f12 − f22 (α22 − m2 fg 2 ) +

1 2 σ1 f = 0, 2 e11

(C2)

1 2 σ2 f = 0, 2 e22

(C3)

− (α11 + α22 − m1 fg1 − m2 fg 2 ) f12 + α12 f22 + α21 f11 + σ1 σ2 fe12 = 0.

(C4)

42

Mathematical Biosciences 296 (2018) 36–44

S. Engen et al.

(

)

1

(

1

)

Eqs. (C2) and (C3) give f11 = α12 f12 + 2 σ12 fe11 / s1 and f22 = α21 f12 + 2 σ22 fe22 / s2, where si = αii − mi fgi , which inserted into Eq. (C4) gives 1

f12 =

1

σ1 σ2 s1 s2 fe12 + 2 σ12 α21 s2 fe11 + 2 σ22 α12 s1 fe22 s12 s2 + s1 s22 − α12 α21 (s1 + s2)

.

(C5)

The covariance functions cij(y) can now be computed numerically for a given set of parameters. Appendix D. Spatial scales of autocovariance functions The spatial scales of these covariance functions are often defined as the standard deviation in a given direction defined by the covariance functions when these are scaled by an appropriate factor to become distributions [26]. Writing fij(ω1, ω2) for the Fourier transforms, these spatial scales lij along the first axis (the x1-axis) are given by [10]

lij2 = −∂2 ln fij / ∂ω12

ω1= ω2 = 0 .

(D1)

Unfortunately, even if these scales can be computed analytically in this model, the resulting formulas for the scales are too complicated for general insight and transparency. However, for the simplified model described in the main text, where all environmental spatial correlations are proportional to a common environmental correlation ρe(y) with Fourier transform fe, we can find some relevant scaling results. Under this model the Fourier transform of c12(y) then takes the simplified form 1

f12 =

1

σ1 σ2 s1 s2 ρ12 (0) + 2 σ12 α21 s2 + 2 σ22 α12 s1 s12 s2 + s1 s22 − α12 α21 (s1 + s2)

fe .

(D2)

le2

2 l12 .

To From Eq. (D2) we then immediately see, because fe is a factor, that the environmental noise gives an additive term to the squared scale find the complete solution we may consider the numerator T(ω) and denominator N(ω) of the fraction separately as each of them also give additive components. Inserting the definitions of αii we get s1 = 1 + μ + r + m1 (1 − fg1 ) and s2 = δ + m2 (1 − fg 2 ), which we write as si = ai + mi (1 − fgi ) for

i = 1, 2 with a1 = 1 + μ + 2νK1 and a2 = δ . Then si (0) = ai , si′ (0) = 0 and si′ ′ (0) = mi lgi2 . Then at ω = 0, we find T = σ1 σ2 a1 a2 ρ12 (0) + σ12 α21 a2 /2 + σ22 α12 a1/2, T ′ = 0 and T ′ ′ = σ1 σ2 a1 ρ12 (0) m2 lg22 + σ1 σ2 a2 ρ12 (0) m1 lg21 + σ12 α21 m2 lg22/2 + σ22 α12 m1 lg21/2. The contribution from the numerator is consequently

− T ′ ′/ T = −A1 m1 lg21 − A2 m2 lg22 where

A1 =

A2 =

2σ1 σ2 a2 ρ12 (0) + σ22 α12 2σ1 σ2 a1 a2 ρ12 (0) + σ12 a2 α21 + σ22 a1 α12 2σ1 σ2 a1 ρ12 (0) + σ12 α21 2σ1 σ2 a1 a2 ρ12 (0) + σ12 a2 α21 + σ22 a1 α12

,

.

Similarly, we find N (0) = a12 a2 + a1 a22 − α12 α21 (a1 + a2), N ′ = 0 and

N ′ ′ = 2a1 a2 m1 lg21 + a12 m2 lg22 + 2a1 a2 m2 lg22 + a22 m1 lg21 − α12 α21 (lg11 + lg22), giving the additive contribution from N as

N ′ ′/ N = B1 m1 lg21 + B2 m2 lg22, where

B1 =

B2 =

2a1 a2 + a22 − α12 α21 , a12 a2 + a1 a22 − α12 α21 (a1 + a2)

a12 a2

2a1 a2 + a12 − α12 α21 . + a1 a22 − α12 α21 (a1 + a2)

This finally yields the squared spatial scale of c12 as 2 l12

= le2 + (B1 − A1 ) m1 lg21 + (B2 − A2 ) m2 lg22.

(D3)

Because a1, a2, α12 and α21 are parameters not depending on the dispersal, the constants Ai and Bi are also independent of dispersal. It is straightforward (although some computations are needed) to check that Bi − Ai > 0 . 1 To find the spatial scale of c11 we consider the Fourier transform f11 = α12 f12 + 2 σ12 fe / s1 and evaluate minus its second derivative with respect to

(

)

ω1 at ω = 0 . Then the term − ln si gives the additive contribution m1 lg21/ a1 to the spatial scale. Writing Ie = fe (0) = taken over the entire two-dimensional space, and similarly I12 = f12 (0) = 2 p11 l12 + q11 le2, where q11 = 1 − p11 , and

p11 =

α12 I12 1

α12 I12 + 2 σ12 Ie

∫ ρ (y ) du, where the integral is ∫ c12 (y ) dy, we find that the contribution from the ln(α12 f12 + 12 σ12 fe ) is

, (D4)

2 which depends on the migration rates as well as the other parameters in the model. Because le2 is an additive component of l12 we see that le2 also ends 2 2 , and l 22 up being an additive component of l11

43

Mathematical Biosciences 296 (2018) 36–44

S. Engen et al. 2 l11 = le2 + p11 [(B1 − A1 ) m1 lg21 + (B2 − A2 ) m2 lg22]

(D5)

In the same way the scale of c22(y) can be expressed by

p22 =

α21 I12 1

α21 I12 + 2 σ22 Ie

, (D6)

as 2 l 22 = le2 + p22 [(B1 − A1 ) m1 lg21 + (B2 − A2 ) m2 lg22]

(D7)

Notice that if there is no migration, m1 = m2 = 0, then l12 = l11 = l22 = le , which we recognize as the Moran effect. Appendix E. Fourier transforms of weighting functions used in simulations of fields 2 From the description of the simulation procedure in the main text it follows that the Fourier transform of c11(y) is f11 = g11 , where now gij are the 2 2 Fourier transforms of the weighting functions ξij. Similarly we find f22 = g12 because the fields dA1(x) and dA2(x) are independent. Finally + g22

2 / f11 . This can then be transformed backf12 = g11 g22 . Hence the Fourier transform of the weights are g11 = f11 , g22 = f12 / f11 , and g12 = f22 − f12 wards by the integration given in Appendix C, finally giving the spatial fields n1(x) and n2(x) in a specific area we want to illustrate by integrations, in practice performed by summation over a grid of squares as explained in the main text.

[20] V. Grøtan, B.-E. Sæther, S. Engen, E.J. Solberg, J.D.C. Linnell, R. Andersen, H.B. seth, et al., Climate causes large-scale spatial synchrony in population fluctuations of a temperate herbivore, Ecology 86 (2005) 1472–1482. [21] C.H. Hsieh, C.S. Reiss, J.R. Hunter, J.R. Beddington, R.M. May, G. Sugihara, Fishing elevates variability in the abundance of exploited species, Nature 443 (2006) 859–862. [22] C.-h. Hsieh, C.S. Reiss, R.P. Hewitt, G. Sugihara, Spatial analysis shows that fishing enhances the climatic sensitivity of marine fishes, Can. J. Fish. Aquat. Sci. 65 (2008) 947–961. [23] N. Jonzén, E. Ranta, P. Lundberg, V. Kaitala, H. Linden, Harvesting-induced population fluctuations? Wildlife Biol. 9 (2003) 59–65. [24] R. Lande, S. Engen, B.-E. Sæther, Optimal harvesting of fluctuating populations with a risk of extinction, Am. Nat. 145 (1995) 728–745. [25] R. Lande, B.-E. Sæther, S. Engen, Threshold harvesting for sustainability of fluctuating resources, Ecology 78 (1997) 1341–1350. [26] R. Lande, S. Engen, B.-E. Sæther, Spatial scale of population synchrony: correlation versus dispersal and density regulation, Am. Nat. 154 (1999) 271–281. [27] R. Lande, S. Engen, B.-E. Sæther, T. Coulson, Estimating density dependence from time series of population age structure, Am. Nat. 168 (2006) 76–87. [28] H. Matsuda, A. Yamauchi, Y. Matsumiya, T. Yamakawa, Reproductive value, harvest value, impact multiplier as indicators for maximum sustainable fisheries, Environ. Econom. Policy Stud. 2 (1999) 129–146. [29] H. Matsuda, P.-A. Abrams, Maximal yields from multi-species fisheries systems: rules for systems with multiple trophic levels, Ecol. Appl. 16 (2006) 225–237. [30] E. Matthysen, Density-dependent dispersal in birds and mammals, Ecography 28 (2005) 403–416. [31] E. McCauley, W.G. Wilson, A.M. Deroos, Dynamics of age-structured and spatially structured predator-prey interactions: individual-based models and populationlevel formulations, Am. Nat. 142 (1993) 412–442. [32] P.A.P. Moran, The statistical analysis of the canadian lynx cycle. II. synchronization and meteorology, Aust. J. Zool. 1 (1953) 291–298. [33] T.J.H. Quinn, R.B. Deriso, Quantitative Fish Dynaminc, Oxford University Press, New York, 1999. [34] W.E. Ricker, Maximum sustained yields from fluctuating environments and mixed stocks, J. Fish. Res. Board Can. 15 (1958) 991–1006. [35] A.O. Shelton, M. Mangel, Fluctuations of fish populations and the magnifying effects of fishing, Proceedings of the National Academy of Sciences of the United States of America 108 (2011) 7075–7080. [36] B.-E. Sæther, S. Engen, E.J. Solberg, Optimal harvest of age structured populations of moose alces alces in a fluctuating environment, Wildlife Biol. 7 (2001) 171–179. [37] B.-E. Sæther, S. Engen, Pattern of variation in avian population growth rates, Philos. Trans. Royal Soc. B London 357 (2002) 1185–1195. [38] S.D. Tuljapurkar, Population dynamics in variable environments, Springer-Verlag, New York, 1990. [39] P. Whittle, J. Horwood, Population extinction and optimal recource management, Philos. Trans. Royal Soc. B 350 (1995) 179–188.

References [1] M. Abramowitz, I. Stegun (Eds.), Handbook of Mathematical Functions, Dover, New York, 1972. [2] C.N.K. Anderson, C.H. Hsieh, S.A. Sandin, R. Hewitt, A. Hollowed, J. Beddington, R.M. May, et al., Why fishing magnifies fluctuations in fish abundance, Nature 452 (2008) 835–839. [3] J.R. Beddington, R.M. May, Harvesting natural populations in a randomly fluctuating environment, Science 197 (1977) 463–465. [4] J.R. Beddington, D.B. Taylor, Optimum age specific harvesting of a population, Biometrics 29 (1973) 801–809. [5] R.J.H. Beverton, S.J. Holt, On the Dynamics of Exploited Fish Populations, 19 Fishery Investigations, London, 1957, pp. 1–533. Ser-2 [6] B. Bolker, S.W. Pacala, Using moment equations to understand stochastically driven spatial pattern formation in ecological systems, Theor. Popul. Biol. 52 (1997) 179–197. [7] L. Ciannelli, J.A.D. Fisher, M. Skern-Mauritzen, M.E. Hunsicker, M. Hidalgo, K.T. Frank, K.M. Bailey, Theory, consequences and evidence of eroding population spatial structure in harvested marine fishes: a review, Mar. Ecol. Prog. Ser. 480 (2013) 227–243. [8] C.W. Clark, Bioeconomic Modelling and Fisheries Management, Wiley-Interscience, New York, 1985. [9] C.W. Clark, Mathematical Bioeconomics. The Optimal Management of Renewable Resources, Wiley, New York, 1990. [10] S. Engen, A dynamic and spatial model with migration generating the log-gaussian field of population densities, Math. Biosci. 173 (2001) 85–102. [11] S. Engen, Spatial synchrony and harvesting in fluctuating populations: relaxing the small noise assumption, Theor. Popul. Biol. 116 (2017) 18–26. [12] S. Engen, B.-E. Sæther, Generalizations of the moran effect explaining spatial synchrony in population fluctuations, Am. Nat. 166 (2005) 603–612. [13] S. Engen, R. Lande, B.-E. Sæther, T. Bregnballe, Estimating the pattern of synchrony in fluctuating populations, J. Anim. Ecol. 74 (2005) 601–611. [14] S. Engen, B.-E. Sæther, Spatial synchrony in population dynamics: the effects of demographic stochasticity and density regulation with a spatial scale, Math. Biosci. 274 (2016) 17–24. [15] M.J. Fogarty, L.W. Botsford, Population connectivity and spatial management of marine fisheries, Oceanography 20 (2007) 112–123. [16] K.T. Frank, B. Petrie, W.C. Leggett, D.G. Boyce, Large scale, synchronous variability of marine fish populations driven by commercial exploitation, Proc. Natl. Acad. Sci. USA 113 (2016) 8248–8253. [17] J.M. Fryxell, C. Packer, K. McCann, E.J. Solberg, B.E. Sæther, Resource management cycles and the sustainability of harvested wildlife populations, Science 328 (2010) 903–906. [18] S.D. Gaines, B. Gaylord, L.R. Gerber, A. Hastings, B.P. Kinlan, Connecting places the ecological consequences of dispersal in the sea, Oceanography 20 (2007) 90–99. [19] W.M. Getz, R.G. Haight, Population Harvesting. Demographic Models of Fish, Forest, and Animal Resources, Princeton, Princeton University Press, 1989.

44