Does averaging overestimate or underestimate population growth? It depends

Ecological Modelling 411 (2019) 108744

Contents lists available at ScienceDirect

Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel

Does averaging overestimate or underestimate population growth? It depends Dmitrii O. Logofeta,b, a b

T

⁎

Laboratory of Mathematical Ecology, A.M. Obukhov Institute of Atmospheric Physics, Russian Academy of Sciences, Moscow 119 017, Russia Institute of Forest Science, Russian Academy of Sciences, Uspenskoe 143030, Russia

ARTICLE INFO

ABSTRACT

Keywords: Discrete-structured population Life cycle graph Matrix calibration Pattern-geometric average Dominant eigenvalue Reproductive uncertainty Stochastic growth rate Monte Carlo simulation

When a matrix population model is nonautonomous, i.e., when it represents a set of single-time-step ("annual") PPMs, L(t), t = 0, 1, …, T – 1, each corresponding to a fixed life cycle graph, then each of the annual matrices generates its own set of model results to characterize the population. In particular, λ1(t), the asymptotic growth rate, varies with t and may result in alternating predictions of population viability. A logical way to characterize the population over the total period of observations is to average the given set of T PPMs, and I have proved the correct mode of averaging to be the pattern-geometric average. It means finding a matrix, G, such that its Tth power equals the product of T annual matrices (in the chronological order), while its pattern does correspond to the given life cycle graph. In practical cases however, the mathematical problem of pattern-geometric average has no exact solution for a fundamental mathematical reason. Nevertheless, the approximate solutions have revealed a fair precision of approximation in recent case studies of alpine short-lived perennials (Eritrichium caucasicum and Androsace albana), resulting in quite certain predictions of population viability by means of λ1(G), the dominant eigenvalue of the average matrix. An alternative approach to estimate the viability leades to the concept of stochastic environment (represented by the given PPMs to be chosen at random) and the ensuing stochastic growth rate, λS. In the case studies, the λSs have been estimated by a direct method of Monte Carlo simulations, all the E. caucasicum estimates having unambiguously been less than λ1(G), whereas those for A. albana being certainly greater than λ1(G). There should be a general mathematical reason for this, too, thus a challenge to the theory of nonnegative matrices.

1. Introduction Theoretical foundations of matrix population models came from the classical Perron–Frobenius theorem for nonnegative matrices (Horn and Johnson, 1990), which introduced the concept of the dominant eigenvalue, λ1(L), and established the existence and uniqueness of λ1(L) > 0 in the spectrum of any (nontrivial) irreducible nonnegative matrix L. Since being launched into the literature as essentially linear systems with λ1s serving as the asymptotic (when t tends to infinity) growth rates (Leslie, 1945, 1948; Lefkovitch, 1965), matrix population models have been motivating nontrivial mathematical problems to be posed and often solved (surveys in Cohen, 1976; Hansen, 1986; Cushing, 1998, and Logofet and Belova, 2008; see also Cushing and Yicang, 1994; Li and Schneider, 2002; Logofet, 2002), including those motivated by the need to calibrate reliably matrix L (Logofet, 2008, 2013a,c,d, 2016,

2017, Logofet et al., 2018; Protasov and Logofet, 2014), which is nowadays called the population projection matrix (PPM, Caswell, 2001). In real populations, however, the asymptotic growth can hardly be observed, therefore apparently dooming λ1 to forever remain a pure theoretical abstraction. Nevertheless, given a matrix L, its λ1 is uniquely determined by the matrix elements (vital rates), i.e., eventually by the population data that have been served to calibrate L. Moreover, λ1(L) is mathematically known to increase/decrease together with any element of (an irreducible) matrix L. It is therefore quite logical that λ1(L) is considered to be the measure of how the local population is adapted to the environment where and when the data have been mined to calibrate the PPM. When a matrix population model deals with the ‘identified individuals’ type of data, where the ‘individuals are marked and followed over time’ (Caswell, 2001, p. 134), the PPM can already be calibrated from the data of only two consecutive time moments, namely, for the

Abbreviations: PPM, Population Projection Matrix; LCG, Life Cycle Graph ⁎ Corresponding author at: Laboratory of Mathematical Ecology, A.M. Obukhov Institute of Atmospheric Physics, Russian Academy of Sciences, Moscow 119 017, Russia. E-mail address: [email protected]. https://doi.org/10.1016/j.ecolmodel.2019.108744 Received 11 September 2018; Received in revised form 5 July 2019; Accepted 8 July 2019 0304-3800/ © 2019 Elsevier B.V. All rights reserved.

Ecological Modelling 411 (2019) 108744

D.O. Logofet

beginning (t) and the end (t + 1) of a single time step (ibidem), in accordance with the basic model equation x(t +1) = L(t) x(t), t = 0, 1, …,

(recent survey in Logofet, 2017). All of them suggest additional assumptions to be made or hypotheses be accepted. Without any additional assumptions/hypotheses, matrix calibration results only in a certain finite range for the unknown reproduction rates, hence in a certain range, [λ1min, 1max ], for the adaptation measure λ1(L) rather than a single number. When the range lies entirely to the left or to the right of λ1 = 1, it does give a certain answer about whether the population declines or increases, but it does not when λ1 = 1 falls inside the range. Both cases did occur in modelling practice (Logofet et al., 2016b, Table 4; Logofet et al., 2017b, Table 3), leaving the issue open to further general studies. Up to date, however, eight annual PPMs have been calibrated and their pattern-geometric average found luckily to have the range of λ1 located entirely to the left of 1 (Logofet et al., 2018, Table 4). The 'identified individuals’ data mined in the A. albana case study do not bear any reproductive uncertainty, thus providing for the PPM calibration in a unique way (Logofet et al., 2017c). Up to date, seven annual PPMs have been calibrated and their pattern-geometric average found (ibidem). In the present article, I report on the 8th calibrated PPM and how it modifies the average. Two outcomes from pattern-geometric averaging in particular case studies give rise to comparing this original mode of population viability estimation with the well-known concepts of ‘population dynamics in variable environments’ (Tuljapurkar, 1990) and stochastic growth rate, λS (ibidem; Caswell, 2001). In the variety of theoretical approximations to λS (ibidem; Sanz, 2019), there is none to cope with reproductive uncertainty, and, moreover, the annual PPMs for E. caucasicum and A. albana have appeared to vary so markedly with time that the validity of theoretical assumptions is dubious for them. I have therefore used direct Monte Carlo simulations to estimate1 λS of the E. caucasicum population (Logofet et al., 2018, Table 5), and I report here on the outcome of the same method applied to the eight PPMs of the A. albana population. Now, after getting the outputs from two methods of viability estimation, we compare them for the (local population of) two alpine species and see the opposite relations: whereas, in case of reproductive uncertainty, the λ1 lower bound of the average PPM for E. caucasicum is certainly greater than the upper one for its λS, the λ1 of the average A. albana matrix is markedly less than the corresponding λS. I discuss possible reasons for this controversy and formulate potential directions to study relevant issues in the theory of nonnegative matrices.

(1)

where x(t) and x(t + 1) represent the observed vectors of population structure. When data for several time steps are available, we get several calibrated PPMs, L(0), L(1), …, L(t), each obeying Eq. (1), with a variety of values for each vital rate, i.e., for each nonzero element of the matrix (Logofet et al., 2017b, Table 3). Different values of matrix elements result in respectively different characteristics (even conflicting ones) of the population, such as the asymptotic growth rate (the dominant eigenvalue, λ1), the steady-state structure (the corresponding eigenvector, x*), and so on (ibidem). The estimates of “age-specific traits from stage-specific models” (Caswell, 2001, p. 116) become different, too (Logofet et al., 2017b, Table 4). If, however, the task is to get the population characteristics ensuing from the total observation period, then we face the need to average all the given PPMs and deduce the characteristics from the average matrix. But unlike positive numbers, how to average matrices is a far not trivial issue, and this issue has not been given much attention in the ecological literature. School education prompts using the arithmetic average in practice (Klimas et al., 2012; Logofet, 2013b), while theoretical constructions just assuming the average matrix to be known (Tuljapurkar, 1986, 1990; Caswell, 2001). Averaging arithmetically poses no problems due to the linearity of matrices as operators in a vector space, but the arithmetic mean has no meaning beside of the arithmetic one. Meanwhile, averaging the PPMs should have certain population meaning, which is explained below. If t = T is the final moment of the observation period, than the observed population structure x(T) is projected from the initial one, x (0), by T consecutive projections obeying Eq. (1): x(T) = L(T – 1) L(T – 2) … L(1) L(0) x(0).

(2)

It is logical that the average matrix, G, does the same by T projections, i.e., GT = L(T – 1) L(T – 2) … L(1) L(0).

(3)

Eq. (3) bears the very idea of how to average the PPMs, and this idea leads evidently to the notion of geometric mean as was noted by Logofet et al. (2017a, 2017b). At first glance, to find G from Eq. (3) is just a matter of routine calculating the T-th power matrix root of the matrix product on the right-hand side. But matrix theory can only guarantee the root to exist in the space of complex matrices, rather than in the set of nonnegative matrices, while the complex values of vital rates make no sense (see further discussion in Logofet et al., 2018). Moreover, the average matrix G is supposed to have the same pattern (of zero/nonzero allocation) as those to be averaged, the pattern being fixed by the known life cycle graph (LCG) of the organisms under study. The above arguments have led me to the concept of pattern-geometric average for nonnegative matrices (Logofet et al., 2017a, b; Logofet et al., 2018). In this paper, I present an efficient method to solve the patterngeometric average problem computationally. The method applies in particular case studies of two alpine short-lived perennials, Eritrichium caucasicum and Androsace albana, among rare, Red-Book spp., the objects of population studies by my former (and hopefully future) botanic coauthors on the permanent sample plots in a high-mountain reserve located in the north-western Caucasus (Logofet et al., 2016a, 2017c). The data mined in the E. caucasicum case study are of a quite specific ‘identified individuals’ type, namely, ‘identified individuals with uncertain parents’ (Logofet, 2010, p. 33; Logofet et al., 2012, p. 98). The ‘uncertain parents’ cause reproductive uncertainty (ibidem), the inability to calibrate stage-specific reproduction rates in a unique way, which complicates estimating the measure of how the local population is adapted to its environment as just the dominant eigenvalue of the calibrated PPM. A number of tricks were proposed in the literature to cope with reproductive uncertainty

2. Objects and methods 2.1. Eritrichium caucasicum case study Eritrichium caucasicum is a herbaceous short-lived perennial plant inhabiting alpine lichen heaths, endemic to the Caucasus Mountains. A local population of E. caucasicum has been monitored annually on the permanent sample plots laid down in the Teberda State Biosphere Reserve (North-Eastern Caucasus), at the altitude of 2800 m asl, in the year of 2009 (Logofet et al., 2017b). Fig. 1 presents the so-called scale of ontogenesis for this species, the sequence of discrete, well-defined stages, a traditional tool in the Russian school of geobotany. Seeds germinate each spring to seedlings and then mature successively to juvenile, virginal, and generative stages, where they produce seeds. Severe wintering conditions cause some of generative plants to not flower in the upcoming season, and this phenomenon causes the backward transition va ⟵ g to appear in the LCG (Fig. 2) corresponding to the one-year time step. Botanists can distinguish between the plants at the generative stage and those at the terminal generative stage, after which the plants die out. Thus, two generative stages provide for the population recruitment at the initial stage, but no 1 Estimation means finding the attainable upper and lower bounds in case of reproductive uncertainty.

2

Ecological Modelling 411 (2019) 108744

D.O. Logofet

Fig. 1. Ontogenetic stages of E. caucasicum: pl, seedlings; j, juvenile plants; v, vegetative plants; g, generative, and gt, terminal generative plants (Logofet et al., 2016a). Fig. 2. LCG for a local population of E. caucasicum: j denotes young plants (seedlings and juveniles); va virginal plants and the adult non-flowering ones, g generative plants; gt terminally generative plants. Solid arrows indicate ontogenetic transitions occurring for one year (the lack of transition, in particular); dashed arrows correspond to annual recruitment (Logofet et al., 2016a).

method is available to distinguish between the progenies from those two stages. Reproductive uncertainty is therefore inevitable in this case, even for the ‘identified individuals’ type of data, which were mined by means of annual censuses from permanent sample plots in the period of 2009 to 2017 (Logofet et al., 2018). Note that eliminating the reproductive uncertainty by means of aggregating two reproductive stages of the life cycle into a single one may distort the true figures up to the opposite forecast of population viability (Logofet, 2017). With the population vector denoted as x(t) = [j(t), v(t), g(t), gt(t)]T and in accordance with the LCG (Fig. 2), the PPM takes on the following form:

0 c L=T+F= 0 0

0 d f k

0 e h l

0 0 0 0 + 0 0 0 0 0 0 0 0

a 0 0 0

b 0 0 0

j(t +1) = a g(t) + b gt(t), t = 0, 1, …, 7,

(5)

where j(t +1), g(t), and gt(t) are all known from the data (Table 1 in Logofet et al., 2018). Known also are all the transitions made by the ‘identified individuals” for one time step (Table 2 in Logofet et al., 2018), so that the PPMs, L(t), can be calibrated in rational numbers, with the observed frequencies of transitions (Table 1). Extremely low values of λ1 in the annual matrices L6 and L7 indicate extremely severe conditions for survival during the years of 2015–2017. The corresponding life cycles have even lost some vital transitions, resulting in zero values of certain vital rates (f, l in L6 and e, h, l in L7). This, in turn, deprives the LCG of being strongly connected (Harary et al., 1965), hence the PPM of being irreducible (Horn and Johnson, 1990), or equivalently indecomposable (Gantmacher, 1959)2 . It means that

(4) 2 A directed graph is strongly connected if any pair of its nodes is connected by a directed path along its arcs. Removing an arc may disturb this property and result in the decomposition of the associated matrix.

with two uncertain reproduction rates, a and b, constrained by the ‘recruitment balance equation’ 3

Ecological Modelling 411 (2019) 108744

D.O. Logofet

Table 1 Eight PPMs L(a, b; t) calibrated from the E. caucasicum data at the moments t, t +1 and the corresponding bounds for the measure of adaptation (expanding Table 3 from Logofet et al., 2016). Census years

2009 j=0

Matrix L(a, b; t)

0 68

L0 =

0 63

149

6

0 0

2010 j=1

2011 j=2

2012 j=3

L1 =

9

0

2

1

136

150

0

4

2014 j=5

6

2015 j=6

2016 j=7

9

0

4

14

0

3

2 3

49

0

0

0

1

73

a b 6 0

2

166

166

1

103

103 103

6

Range of λ1(Lj)

Eigenvector x *j , %

λ1min

λ1max

10a + 4b = 31; [0, 31/10], [0, 31/4]

0.9035

0.9949

9.39 79.62 9.89 1.10

23.11 68.72 7.42 0.75

9a + b = 150; [0, 150/9], [0, 150]

1.2460

1.5201

42.04 53.90 3.14 0.92

54.96 42.49 1.99 0.56

10a + 3b = 211; [0, 211/10], [0, 211/3]

1.2476

1.4948

44.58 50.88 2.64 1.90

56.16 40.93 1.72 1.19

9a + 7b = 119; [0, 119/9], [0, 17]

0.9213

1.1004

7.15 89.25 3.21 0.39

25.94 71.68 2.16 0.22

6a + b = 99; [0, 99/6], [0, 99]

0.7864

0.8588

56.49 40.13 2.69 0.69

62.82 34.63 2.00 0.55

11a + 4b = 49; [0, 49/11], [0, 49/4]

0.8376

0.9119

38.81 48.88 7.30 5.01

48.57 42.25 5.57 3.61

17a + 8b = 73; [0, 73/17], [0, 73/8]

0.6719

0.6987

0 98.85 0 1.15

12.62 86.41 0 0.97

a + b = 13; [0, 13], [0, 13]

0.4541

0.5970

35.17 57.46 6.14 1.23

61.74 34.86 2.83 0.57

x *jmin

x *max j

0

9

4

0

0 0

296

0 5

0

a b 9 0

296

128

0

0

10

0 0

128

0

10

1

128

39

b 0

0

0 86

10

6

181

166

0 L7 =

a

296

103

0

9

0

4

0

99

0 L6 =

0

9

181

139

0

9

9

0

119

0 L5 =

a b 9 0

6

0 153

0

22

0

129

0 L4 =

10

129

0

23

0

129

0 2013 j=4

10

0 101 7

211

b 0

1

136

0

137

10

1

136

0

L3 =

3

80

106

0

L2 =

5

0

31

76

0

80

0

0 17

Equation (5); range of a, range of b

a 3 3 5

11

b 0

11

0

11

0

a 8 1 0

17

b 0

17

0

17

0

a b 1 0

0 0 0

1

0

1

0

regular routines to calculate eigenvalues and eigenvectors may fail in finding the dominant one (Logofet, 2013b) and we need a closer look at the strongly connected components into which the graph has disintegrated and to check which of the corresponding matrix blocks bears the dominant eigenvalue. It should preferably be a diagonal block that contains reproduction rates if such a block does exist. It does not, for example, in matrix L6 when a = 0 or b = 0 and the matrix decomposes just into its diagonal elements, so that λ1 becomes the maximal of them. Zero components in the dominant eigenvectors do not contradict the Perron–Frobenius theorem in this case as the theorem applies only to irreducible matrices (see Caswell, 2001, §§ 4.5.1, 4.5.4, 7.5.2 for greater detail). Among the 8 calibrated PPMs, there are 2 having the range of λ1 located entirely to the right of λ1 = 1, 5 to the left of 1, and one matrix (L3) having 1 inside of the range. Intuition prompts that both the average matrix and the stochastic growth rate should indicate population decline in this case.

2.2. Averaging under reproductive uncertainty In general, the averaging matrix equations like Eq. (3) have no exact solutions as overdetermined systems of scalar equations. For instance, Eq. (3) takes on the following form: G8 – L7 L6 L5 L4 L3 L2 L1 L0 = 0

(6)

for E. caucasicum, with matrices Lj being determined in Table 1. While corresponding to the LCG, the average matrix G has only 9 unknown elements (see Eq. (4)), whereas the system of element-wise scalar equations for them (the one equivalent to the matrix Eq. (6)) consists of 16 equations, which can therefore not be consistent. So, we have to resort to an approximate solution, the one that minimizes the difference between the (chronological) product of 8 PPMs and the 8th power of the average matrix. The task to find the pattern-geometric average of the given matrices reduces therefore to a constrained minimization problem for a norm of 4

Ecological Modelling 411 (2019) 108744

D.O. Logofet

Fig. 3. Ontogenetic stages of Androsace albana: pl, plantules; j, juvenile plants; im, immature plants; v, adult vegetative plants, and g, generative plants (Logofet et al., 2017c).

the left-hand side of Eq. (6), the constraints being determined by matrix G being a PPM of the fixed pattern (4) and within the bounds of the PPMs to be averaged. Greater detail and routine descriptions can be found in Logofet et al. (2018). Reproductive uncertainty complicates the task as we have to deal with the whole one-parametric families3 instead of single annual matrices, with certain ranges of their dominant eigenvalues (Table 1), hence the average matrix has to also represent a respective family with a certain range of its λ1s. To avoid a variety of mathematical settings for the pattern-geometric average problem in this case, I have proposed a heuristic approach using the T + F decomposition as in Eq. (4) (Logofet et al., 2018). We first solve the constraint maximization problem for the transition parts Tj (with a high approximation accuracy guaranteed by their substochastic property), then we add the average of Fjs found by means of interval computations (e.g., Neumaier, 1990). This TF-method has been applied to solve the pattern-geometric average problem for 7 PPMs of E. caucasicum (Logofet et al., 2018, Table 2), with the range [λ1min(GTF), 1max (GTF)] lying luckily on the right of λ1 = 1 (Logofet et al., 2018), and I report here the result of adding the 8th, recently calibrated PPM, to the former set of annual PPMs to be averaged.

corresponding set of vital rates (matrix elements). Stochastic variations in time are therefore reduced to random choice from the given set of PPMs in a long (theoretically infinite) sequence of transformations, Lj, that act on the population vector. Such a sequence was long ago proved to converge (Furstenberg and Kesten, 1960; Oseledec, 1968; Cohen, 1979a, 1979b) to the limit that defined what is now called the stochastic growth rate, λS:

log s= lim

1

log N ( ) = lim

1

log L

1…L 0 x (0) 1

(7)

Here N(τ) denotes the total population size and || … ||1 the norm as the sum of vector components in modulus. Evidently, the value of λS relative to 1 indicates exponential growth or decline in the total population size, thus serving as an index of population viability. Theoretically, the limit of a converging sequence can be accurately approximated by any distant enough member, and this is the method I apply to estimate λS in a number of Monte Carlo tests. The random choice is i.i.d. (independent, identically distributed) at each step τ, while the reproductive uncertainty forced me to choose uncertain reproduction rates, a and b, at random, too. Therefore, I implement two series of tests: the first one at the uniform probability distribution over the Lj-specific range of a (hence b being uniquely determined after a has been chosen), the second series at the normal distribution truncated on the same finite range. Repetitions of each series reveal the minimal and maximal estimates of λS for each version of simulated reproductive uncertainty.

2.3. Stochastic growth rate as an index of population viability In the theory of population growth in random environments (Tuljapurkar, 1990; Caswell, 2001), the set of calibrated annual PPMs is considered to reflect the spectrum of potential variations in the environment, each particular matrix being responsible for the whole complex of environmental conditions that determined the

2.4. Androsace albana case study Androsace albana is also a herbaceous short-lived perennial plant inhabiting alpine lichen heaths and other habitats (Logofet et al., 2017c). It is listed in the "Red Book of the Krasnodar Territory" (2007) and the "Red Book of the Adygea Republic" (2012) as a rare species. Its

3 Two parameters, a and b, are bound by the linear Eq. (5), whereby the second can be expressed via the first or vice versa.

5

Ecological Modelling 411 (2019) 108744

D.O. Logofet

Fig. 4. LCG for a local population of A. albana: the same notations as in Fig. 3; dashed arrows indicate the recruitment found in the corresponding stages (Logofet et al., 2017c).

scale of ontogenesis is shown in Fig. 3 and the LCG in Fig. 4 bearing the same features as the LCG of E. caucasicum. In particular, delays ↺ in the im and v stages explainable by the fact that the harsh conditions of high mountains force plants to resort to the "space-holder strategy" (Körner, 2003) meaning to stay or grow in one place for as long as possible; accelerated transitions pl →im as a manifestation of polyvariant ontogeny in A. albana under conditions of the alpine zone in the northwestern Caucasus. However, the species is monocarpic unlike E. caucasicum, which excludes any chance of reproductive uncertainty. The ‘identified individuals’ data have been mined from permanent sample plots each August in the 2009–2017 period (Logofet et al., 2019). The ‘individuals’ exhibit different rates of ontogeny, so that population recruitments are found in either of the stage-specific groups pl, j, or im at the moment of census. In accordance with the LCG (Fig. 4), the PPM takes on the following form:

0 d L= e 0 0

0 0 f 0 0

0 0 h k 0

0 0 0 l m

a b c ; a, b , …,l, m 0 0

Table 2 Eight PPMs L(t) calibrated from the A. albana data at the moments t, t +1, the corresponding dominant eigenvalues and their normalized eigenvectors (expanding Table 3 from Logofet et al., 2017c). Census years 2009 t=0

Matrix L(t)

0

0

0

0

30

37

0

0

0

40

0

3

8 2

2010 t=1

0

Again, eight PPMs, L(t), have been calibrated in rational numbers, with the observed frequencies of transitions (Table 2). Now 3 calibrated PPMs have their λ1s located to the right of λ1 = 1 and 5 to the left of 1, while the favourable/unfavourable years seem to synchronise for both species (see Table 1). The same intuition prompts again that both the average matrix and the stochastic growth rate should indicate population decline.

2012 t=3

0

0

0

0

19

0

0

0

31

0

0

Table 3 presents the outcome of solving the pattern-geometric average problem by the TF-method explained in Section 2.2. The range [λ1min(GTF), 1max (GTF)] is still (after adding the 8th annual PPM to those 7 averaged before, Table 2 in Logofet et al., 2018) on the left of λ1 = 1. The steady-state population structures (dominant eigenvectors) are apparently within the bounds determined by the annual PPMs.

30

0

0

0

55

1

0

0

19

0

0

85

0

25

35

21

45

0

0

0

0

0

43

10

48

43

4

0

57

0

57

0

0

0

0

0

0

0

0

0

0

0 16

19

0

137

0

0

0

0

0

14 6

0

0

0

16

2

98

0

0

0

0

6 4

95

0 0 16

34

4

50

0

50

0

0

0

10

0

0

0

0

29

0

7

3

3

19

2

0 0

10

14

10

2

0

4

3

0

0

0

0

8

0

0

0

0

0

0

4

13 13

0

14.45

0.3988

11.71 56.64 11.69 17.46 3.50

12.64

0.8382

11.13 32.27 21.60 31.27 3.73

8.17

0.8843

7.78 20.75 1.96 64.92 4.59

43.79

6

6

0

20

0

2

7.76 47.54 14.34 28.51 1.85

4

0

29

0.6345

0

20

0 2

13.06

4

0 0

6

3 19 3 0 3

0

6.01 43.15 29.67 19.56 1.60

0 4

0

1.2641

0

73

73

34

22.15

0

2

0

17.20 30.40 39.39 12.45 0.55

0

0

0

1.5779

1

0

0

18.97

1

98

3

0

1

16

44

95

15.90 31.99 18.25 32.83 1.03

1

0

0

1 1

19 4 136 4 1 4

49 10 45 39 0 49 86 87 0 0 28 87 45 58 6 0 0 0 58

0

1.2283

100x (t ) x (t ) 1

0

26

0

19

39.68

0

26

0

0

6

23

0

0 2016 t=7

55

34

10.61 18.20 16.99 51.59 2.60

Δ(

0

35

0

0

Random environments are represented by the 8 annual PPMs calibrated from the ‘identified individuals with uncertain parents’ data (Table 1). Chosen at random (i.i.d.), the matrices form fairly long sequences to estimate limit (7). Being random, different realisations of the finite sequence give different, yet close to the limit, estimates of λS, the greater the number of realisations, the wider should be the range of the estimates. Table 4 represents the results illustrating this point.

1

0.5661

13

0

35

49

2

3.2. E. caucasicum growth in random environments

17

48

0

2014 t=5

2015 t=6

22

30

99

Vector x*, %

13

0

2

3.1. Averaging the annual E. caucasicum PPMs

19

0

1

2013 t=4

3. Results

99

0

6

(8)

7

0

14

1

28

110

0

4

2011 t=2

22

37

13

λ1(L(t))

14 1

2 2

16

0

16

0

– x*)

Ecological Modelling 411 (2019) 108744

D.O. Logofet

As was expected, all the estimated ranges of variations in the estimates locate to the left of 1, indicating population decline. It is natural that different realizations of the random sequence N(τ) give different estimates of λS and their range narrows when the “length” of the finite sequence increases (Table 4). Conversely, the increase in the number of realizations within 100, while preserving the “length”, somehow widens the range of variations (Table 4). The ranges obtained in Series 2 (with normal distributions over the ranges of uncertainty in the reproduction rates) are consistently narrower than those in Series 1 (with uniform distributions). Comparing these ranges with that of the patterngeometric average (Table 3), we see the latter overestimating all the ranges of λS thus obtained.

Table 3 Averaging the annual E. caucasicum PPMs (Table 1) in the form of G(a, b) = Tav + Faa(a, b). Faa(a, b)

0 0 0 0

a, b Average: Tav

Approxim. errora λ1(G)min λ1(G)max * x min ,%

x

max

*

a

2a/23+ b/47 = 1 pattern-geometric, Tpg 0 0 0 0.366717 0.632832 0.422436 0 0.060145 0.217835 0 0.008751 0.069178 5.192…⋅10–7 0. 945004 0.989230 [40.79 54.04 4.35 0.82]T [44.69 50.75 3.84 0.72]T

,%

0 0 0 0

0 0 0 0

a 0 0 0

b 0 0 0

arithmetic, Taa 0 0 0 0 0.3604 0.6670 0.4554 0 0 0.0490 0.1595 0 0 0.0119 0.1346 0 0.000…0 0. 959403 1. 014798 Not calculated

3.3. Averaging the annual A. albana PPMs

Not calculated

Table 5 presents the outcome of solving the pattern-geometric average problem by minimizing directly the norm of the left-hand side in Eq. (6) for the 8 PPMs presented in Table 2. The minimal value of the norm is considered as the approximation error, and it has turned out fairly small, equal to 2.941…⋅10–3. Again, the arithmetic mean of

The minimal value found in the constraint minimization problem.

Table 4 Estimating the stochastic growth rate, λS, of the local E. caucasicum population by the direct Monte Carlo technique in two series of experiments (adapted from Logofet et al., 2018). Product “length”a)

Number of random realisations

Range of variations in the estimates of λS; range length Series 1

1 × 10

5

13 33 100 13 33 100 13 33 100 13 33 100

2 × 105 3 × 105 5 × 105

[0.933260, [0.933260, [0.932331, [0.933893, [0.933893, [0.933284, [0.933736, [0.933736, [0.933718, [0.934529, [0.934194, [0.934156,

Series 2 0.936427]; 0.937678]; 0.937678]; 0.937214]; 0.937214]; 0.937214]; 0.936813]; 0.936813]; 0.936813]; 0.935926]; 0.935949]; 0.936192];

0.003167 0.004418 0.005347 0.003321 0.003321 0.003930 0.003077 0.003077 0.003095 0.001397 0.001755 0.002036

[0.935261, 0.938907]; 0.003646 [0.935108, 0.938907]; 0.003799 [0.935108, 0.939106]; 0.003998 [0.935965, 0.937391]; 0.001426 [0.935965, 0.937845]; 0.001881 [0.935593, 0.937899]; 0.002306 [0.935766, 0.936946]; 0.001181 [0.935766, 0.936946]; 0.001181 [0.934975, 0.937149]; 0.002174 [0. 936473, 0.937200]; 0.000727 [0.936261, 0. 937200; 0.000938 [0.935831, 0. 937225]; 0.001394

a The number of a finite member the sequence {N(τ)}that approximates its limit (7) as τ → ∞; it coincides with the number of cofactors in the product of randomly chosen matrices Lτ–1 Lτ–2 … L1 L0 that still retains vector x(τ) from becoming the computer zero due to normalizing, at each step, by coef = 0.93497…, a specially selected "scaling factor" (for Series 2 coef = 0.93599…); greater detail in Logofet et al. (2018).

Table 5 Pattern-geometric average of the annual A. albana PPMs (Table 2) in comparison with the arithmetic average. Average:

pattern-geometric, Lav = G

arithmetic, Lav = Laa

Lav

0 0 0 0 2.7087 0.3678 0 0 0 6.1163 0.0433 0.3321 0.1474 0 0.0000 0 0 0.1802 0.7940 0 0 0 0 0.0439 0 –3 2.94122⋅10 0.83580 [8.89 23.98 12.13 52.26 2.74]T

0 0 0 0 10.3822 0.0945 0 0 0 23.37&42 0.1172 0.2904 0.2883 0 3.4455 0 0 0.2415 0.6929 0 0 0 0 0.0655 0 0.00…0 0.97524 [15.44 36.27 25.25 21.59 1.45]T

Approxim. error λ1(G) x*, %

7

Ecological Modelling 411 (2019) 108744

D.O. Logofet

references therein), the direct method reducing to the calculation of a distant enough member of the sequence should be regarded as the most reliable and consistent one: its consistency is the same as the fundamentals of real analysis in mathematics. The only point to be criticized is the random i.i.d. choice of the current PPM, which is, of course, a caricature of reality and should be modified to mimic the real course of “good and bad” years (Sanz, 2019). But even two caricatures should look alike when drawn by the same technique. As regards the pattern-geometric average, it is a novel concept wellgrounded at least in the population projection terms (Introduction, Section 2.2) and motivating the proper mathematics (Eqs. (3) and (6)). Thinking just mathematically, one might (and some reviewers did!) perceive the task as solving the mere equation

Table 6 Estimating the stochastic growth rate, λS, of the local A. albana population by the direct Monte Carlo technique. Product “length”a)

Number of random realisations

Range of variations in the estimates of λS; range length

1 × 105

13 33 100 13 33 100 13 33 100 13 33 100

[0.872940, [0.872293, [0.872293, [0.873760, [0.873759, [0.873759, [0.874439, [0.874238, [0.874110, [0.874859, [0.874859, [0.874138,

2 × 105 3 × 105 5 × 105

a

0.879010]; 0.879010]; 0.879010]; 0.876120]; 0.877356]; 0.877719]; 0.876584]; 0.876949]; 0.877002]; 0.875923]; 0.876113]; 0.876138];

0. 006070 0. 006717 0.006717 0.002360 0.003597 0.003960 0.002145 0.002711 0.002892 0. 001064 0.001254 0. 002000

G = (L7 L6 L5 L4 L3 L2 L1 L0)1/8,

(9)

which is formally equivalent to Eq. (6) and represents the idea of geometric mean for the 8 given matrices. However, the matrices are generally not commutative in the product4, unlike real numbers, and it is the chronological order only (right to left in the multiplication) that makes sense in the population projection context (Eq. (2)). Even theoretically, to certify that a proper matrix root does exist in the variety of formal solutions to Eq. (9) is a hard problem, and matrix theory has not yet achieved a great progress in the problem of geometric mean. Certain progress has only been achieved for those positive definite matrices that commute in the product (Ando et al., 2004) and a few results found for those that do not (Bhatia and Holbrook, 2006). By no means can the PPM be positive definite already because of its not being symmetric. The problem to find the matrix root of a given natural power has been permanently attracting certain attention in the literature on matrix theory (Jain et al., 1979; Denman, 1981; Reams, 1997; Higham and Lin, 2009; McDonald et al., 2014; Politi and Popolizio, 2015; Tam and Huang, 2016), treated mostly as a particular case of finding the matrix function of a given matrix via its Jordan canonical form (Horn and Johnson, 1990, Ch. 9; Higham, 2008, Ch. 7). In particular, it was noted that many questions about nonnegative matrices with a positive eigenvector could be reduced to the questions about stochastic matrices (Horn and Johnson, 1990, § 8.7.P1), linking to applications with discrete-time Markov chain models. However, in computation practice, even if you succeeded in getting the arithmetic value of the T-th power root from the diagonal Jordan canonical form of the PPM product and thus obtain the real-valued (or even positive-valued) matrix root, there would be nothing to guarantee the result having the same matrix pattern as the PPMs to be averaged. Moreover, such a matrix does generally not exist as is noted in Section 2.2, with regard to Eq. (6). In these circumstances, the approximate solution to Eq. (6) as a matrix whose pattern does correspond to the given LCG, while the elements being within the bounds predetermined by the annual PPMs to be averaged, represents a reasonable decision. The constraint minimization problem that has implemented the decision finds fairly low values of the approximation error (Tables 3 and 5), much lower than those that might explain the change in relationships between the viability indices under concern by the inaccuracy of approximate solutions. Thus, the methods themselves can hardly be the source of controversy, and we have to seek out along another direction, from the Methods to Objects, i.e., the local populations of two species. There are obvious distinctions in the species biology, A. albana being monocarpic, but E. caucasicum being polycarpic, with two reproductive stages in the life cycle (Fig. 2) and the lack of known traits to distinguish the recruitment from the parents at different stages (Logofet et al., 2017b). Ensuing reproductive uncertainty in the field data modifies slightly the experiment design in Monte Carlo simulations, but can hardly cause any

See the footnote to Table 4; coef = 0.874966….

annual PPMs gives markedly greater value of λ1, with a markedly differing dominant eigenvector. 3.4. A. albana growth in random environments Random environments are presented by the 8 annual PPMs calibrated from the ‘identified individuals’ data (Table 2). The same technique as before, yet without random choice of the reproduction rates, results in the estimates of λS shown in Table 6. Again, all the estimates locate to the left of 1, indicating population decline. It is natural that different realizations of the random sequence N(τ) give different estimates of λS, and the range narrows (from 6 to 1 thousandth parts) with increasing the sequence “length”, illustrating a fundamental lemma of real analysis (Table 6). Conversely, the increase in the number of realizations within 100, while preserving the “length”, somehow widens the range of variations, illustrating just common sense (Table 6). However, the A. albana stochastic growth rate turns out markedly greater that λ1(G) (Table 5), the dominant eigenvalue of the patterngeometric average of the same annual PPMs, in contrast to the opposite relation in E. caucasicum. This observation is what the common sense cannot explain. 4. Discussion and conclusion Two species populations of herbaceous plants are monitored annually in the same alpine habitat, and it is not surprising that their PPMs calibrated from the data of (almost) the same type react on temporal variations of environment conditions in a concerted manner: the years when the E. caucasicum PPMs have (at least the maximum of) their λ1s greater/less than 1 are the same as those for A. albana (cf. Tables 1 and 2). As the most of PPMs have their λ1s essentially less that 1, it is also not surprising that both the λ1 of their pattern-geometric averages and the stochastic growth rates, these indices of local population viability, are markedly less than 1 (cf. Tables 3 and 5), predicting the populations to decline. Surprising is the fact the species-specific investigations of those indices reveal the opposite order relationships between them in the populations under concern (cf. Tables 3,4 and 5, 6), which deprives the pattern-geometric averaging (as a less resourceconsuming method) of an ability to serve as an a priori estimate of the stochastic growth rate. The first direction to seek out the reason for such an opposition concerns the methods that both indices are calculated by. As regards the stochastic growth rate, the method that I have applied ensues directly from the mathematical definition of λS as the limit of the infinite sequence in the left-hand side of Eq. (7). While there exist theoretical approximations of (or bounds for) the λS (see Sanz, 2019, and

4

8

It means that a permutation of cofactor matrices does change the product.

Ecological Modelling 411 (2019) 108744

D.O. Logofet

principal distinctions in the outcome of the direct method to estimate the λS. Hence, the reason of controversy should be sought in averaging the annual PPMs, among the effects caused by reproductive uncertainty. Unfortunately, matrix theory cannot, at the moment, propose any theoretical results that might elucidate how, given a set of annual PPMs, the pattern of the LCG would effect the λ1 of G, their pattern-geometric average matrix as a solution to the constraint minimization problem (Section 2.2). The mathematical problem itself is not trivial (Logofet, 2008), which complicates theoretical research. It is still not clear whether the solution ought to belong the boundary or the interior of B, the set of constraints, and this uncertainty diminishes reliability in the computational procedures that should find the global minimum of the approximation error on B. In view of high demands the direct method of lambda estimation places on computing resources, a “dream” result would be a theorem that could enable us, given a particular LCG to judge whether the pattern-geometric average would give the upper of the lower bound for the corresponding stochastic growth rate. This is a great challenge to the theory of nonnegative matrices and an issue worthy of research efforts in a number of aspects discussed.

33, 183–212. Leslie, P.H., 1948. Some further notes on the use of matrices in population mathematics. Biometrika 35, 213–245. Li, C.-K., Schneider, H., 2002. Application of Perron–Frobenius theory to population dynamics. J. Math. Biol. 44, 450–462. Logofet, D.O., 2002. Three sources and three constituents of the formalism for a population with discrete age and stage structures. Matematicheskoe Modelirovanie (Mathematical Modelling) 14 (12), 11–22 (in Russian, with English summary). Logofet, D.O., 2008. Convexity in projection matrices: projection to a calibration problem. Ecol. Modell. 216, 217–228. Logofet, D.O., 2010. Svirezhev’s principle of substitution and matrix models of the dynamics of populations with a complex structure. Zhurnal Obshchei Biologii (J. General Biology) 71 (1), 30–40 (in Russian, with English summary). Logofet, D.O., 2013a. Complexity in matrix population models: polyvariant ontogeny and reproductive uncertainty. Ecol. Complex. 15, 43–51. Logofet, D.O., 2013b. Projection matrices in variable environments: λ1 in theory and practice. Ecol. Modell. 251, 307–311. Logofet, D.O., 2013c. Calamagrostis model revisited: matrix calibration as a constraint maximization problem. Ecol. Modell. 254, 71–79. Logofet, D.O., 2013d. Projection matrices revisited: a potential growth indicator and the merit of indication. J. Math. Sci. 193 (5), 671–686. Logofet, D.O., 2016. Estimating the fitness of a local discrete-structured population: from uncertainty to an exact number. Ecol. Modell. 329, 112–120. Logofet, D.O., 2017. Aggregation may or may not eliminate reproductive uncertainty. Ecol. Modell. 363, 187–191. Logofet, D.O., Belova, I.N., 2008. Nonnegative matrices as a tool to model population dynamics: classical models and contemporary expansions. J. Math. Sciences 155 (6), 894–907. Logofet, D.O., Ulanova, N.G., Belova, I.N., 2012. Two paradigms in mathematical population biology: an attempt at synthesis. Biol. Bull. Rev. 2 (1), 89–104. Logofet, D.O., Belova, I.N., Kazantseva, E.S., Onipchenko, V.G., 2016a. Local population of Eritrichium caucasicum as an object of mathematical modelling. I. Life cycle graph and a nonautonomous matrix model. Zhurnal Obshchei Biologii (J. General Biology) 77 (2), 106–121 (in Russian, with English summary). Logofet, D.O., Ulanova, N.G., Belova, I.N., 2016b. Polyvariant ontogeny in woodreeds: novel models and new discoveries. Biol. Bull. Rev. 6 (5), 365–385. Logofet, D.O., Ulanova, N.G., Belova, I.N., 2017a. From uncertainty to an exact number: developing a method to estimate the fitness of a clonal species with polyvariant ontogeny. Biol. Bull. Rev. 7 (5), 379–396. Logofet, D.O., Belova, I.N., Kazantseva, E.S., Onipchenko, V.G., 2017b. Local population of Eritrichium caucasicum as an object of mathematical modelling. I. Life cycle graph and a nonautonomous matrix model. Biol. Bull. Rev. 7 (5), 415–427. Logofet, D.O., Kazantseva, E.S., Belova, I.N., Onipchenko, V.G., 2017c. How long does a short-lived perennial live? A modeling approach. Zhurnal Obshchei Biologii (J. General Biology) 78 (5), 63–80 (in Russian, with English summary). Logofet, D.O., Kazantseva, E.S., Belova, I.N., Onipchenko, V.G., 2018. Local population of Eritrichium caucasicum as an object of mathematical modelling. III. Population growth in random environments. Zhurnal Obshchei Biologii (J. General Biology) 79 (4), 249–261 (in Russian, with English summary). Logofet, D.O., Kazantseva, E.S., Belova, I.N., Onipchenko, V.G., 2019. Disappointing survival forecast for a local population of Androsace albanain in a random environment. Zhurnal Obshchei Biologii (J. General Biology) 80 (in press) (in Russian, with English summary). McDonald, J.J., Paparella, P., Michael, J., Tsatsomeros, M.J., 2014. Matrix roots of eventually positive matrices. Linear Algebra Appl. 456, 122–137. Neumaier, A., 1990. Interval Methods for Systems of Equations (Encyclopedia of Mathematics and Its Applications, vol. 37 Cambridge Univ. Press, Cambridge, UK. Oseledec, V.I., 1968. A multiplicative ergodic theorem: ljapunov characteristic numbers for dynamical systems. Trans. Moscow Math. Soc. 19, 197–231. Politi, T., Popolizio, M., 2015. On stochasticity preserving methods for the computation of the matrix pth root. Math. Comput. Simul. 110, 53–68. Protasov, V..Yu, Logofet, D.O., 2014. Rank-one corrections of nonnegative matrices, with an application to matrix population models. SIAM J. Matrix Anal. Appl. 35 (2), 749–764. Reams, R., 1997. A Galois approach to mth roots of matrices with rational entries. Linear Algebra Appl. 258, 1287–1294. Red Book of the Krasnodar Territory, 2007. In: Litvinskaya, S.A. (Ed.), Red Book of the Krasnodar Territory (Plants and Mushrooms)), 2nd ed. Design Bureau, Krasnodar No.1 (in Russian). Red Book of the Adygea Republic, 2012. In: Zamotaylov, A.S. (Ed.), Red Book of the Adygea Republic: Rare and Endangered Objects of Fauna and Flora: in 2 Parts, Part 1: Plants and Fungi, 2nd ed. Kachestvo, Maykop (in Russian). Sanz, L., 2019. Conditions for growth and extinction in matrix models with environmental stochasticity. Ecol. Modell (in press). Tam, B.-S., Huang, P.-R., 2016. Nonnegative square roots of matrices. Linear Algebra Appl. 498, 404–440. Tuljapurkar, S.D., 1986. Demography in stochastic environments. II. Growth and convergence rates. J. Math. Biol. 24, 569–581. Tuljapurkar, S.D., 1990. Population Dynamics in Variable Environments. Springer, New York.

Acknowledgements The study was supported by the Russian Foundation for Basic Research, grants #16-04-00832, #19-04-01227. Roberto SalguroGomez and two anonymous reviewers made certain corrections and suggestions that have improved the quality of the original manuscript. References Ando, T., Li, C.-K., Mathias, R., 2004. Geometric means. Linear Algebra Appl. 385, 305–334. Bhatia, R., Holbrook, J., 2006. Noncommutative geometric means. Math. Intell. 28, 32–39. Caswell, H., 2001. Matrix Population Models: Construction, Analysis and Interpretation, 2nd ed. Sinauer, Sunderland, MA, 722 pp. Cohen, J.E., 1976. Ergodicity of age structure in populations with Markovian vital rates, I: countable states. J. Am. Stat. Assoc. 71 (354), 335–339. Cohen, J.E., 1979a. Ergodicity of age structure in populations with Markovian vital rates. II. General states. Adv. Appl. Prob. 9, 18–37. Cohen, J.E., 1979b. Ergodicity of age structure in populations with Markovian vital rates. III. Finite-state moments and growth rate; an illustration. Adv. Appl. Prob. 9, 462–475. Cushing, J.M., 1998. An Introduction to Structured Population Dynamics (CBMS-NSF Regional Conference Series in Applied Mathematics, 71). SIAM. Cushing, J.M., Yicang, Z., 1994. The net reproductive value and stability in matrix population models. Nat. Resource Model. 8, 297–333. Denman, E.D., 1981. Roots of real matrices. Linear Algebra Appl. 36, 133–139. Furstenberg, H., Kesten, H., 1960. Products of random matrices. Ann. Math. Statist. 31, 457–469. Gantmacher, F.R., 1959. Matrix Theory. Chelsea Publ., New York, NY. Hansen, P.E., 1986. Leslie matrix models: a mathematical survey. In: Csetenyi, A.I. (Ed.), Papers on Mathematical Ecology. Karl Marx Univ. Economics, Budapest, pp. 54–106. Harary, F., Norman, R.Z., Cartwright, D., 1965. Structural Models: an Introduction to the Theory of Directed Graphs. John Wiley, New York 415 pp. Higham, N.J., 2008. Functions of Matrices: Theory and Computation. SIAM, Philadelphia, PA. Higham, N.J., Lin, L., 2009. On Pth Roots of Stochastic Matrices. MIMS EPrint, pp. 21. http://www.manchester.ac.uk/mims/eprints. Horn, R.A., Johnson, C.R., 1990. Matrix Analysis. Cambridge University Press, Cambridge 561 pp. Jain, S.K., Goel, V.K., Kwak, E.K., 1979. Nonnegative mth roots of nonnegative 0-symmetric idempotent matrices. Linear Algebra Appl. 37–51. Klimas, C.A., Cropper Jr, W.P., Kainer, K.A., de Oliveira Wadt, L.H., 2012. Viability of combined timber and non-timber harvests for one species: a Carapa guianensis case study. Ecol. Modell. 246, 147–156. Körner, C., 2003. Alpine Plant Life: Functional Plant Ecology of High Mountain Ecosystems, 2nd ed. Springer, Berlin. Lefkovitch, L.P., 1965. The study of population growth in organisms grouped by stages. Biometrics 21, 1–18. Leslie, P.H., 1945. On the use of matrices in certain population mathematics. Biometrika

9