Potential-Growth Indicators Revisited: Higher Generality and Wider Merit of Indication

: The notion of a potential-growth indicator came to being in the ﬁeld of matrix population models long ago, almost simultaneously with the pioneering Leslie model for age-structured population dynamics, although the term has been given and the theory developed only in recent years. The indicator represents an explicit function, R ( L ), of matrix L elements and indicates the position of the spectral radius of L relative to 1 on the real axis, thus signifying the population growth, or decline, or stabilization. Some indicators turned out to be useful in theoretical layouts and practical applications prior to calculating the spectral radius itself. The most senior (1994) and popular indicator, R 0 ( L ), is known as the net reproductive rate, and we consider two others, R 1 ( L ) and R RT ( A ), developed later on. All the three are different in terms of their simplicity and the level of generality, and we illustrate them with a case study of Calamagrostis epigeios , a long-rhizome perennial weed actively colonizing open spaces in the temperate zone. While the R 0 ( L ) and R 1 ( L ) fail, respectively, because of complexity and insufﬁcient generality, the R RT ( L ) does succeed, justifying the merit of indication.


Introduction
The concept of the potential-growth indicator (PGI) was developed in the theory of matrix population models (MPMs) for the dynamics of discrete-structured biological populations [1,2]. This kind of model represents the basic tool in the mathematical demography of plant and animal populations structured with regard to a certain classification trait such as the age, size, or developmental stage of individuals in a local population of a given species [1,2]. Mathematically, the MPM is a system of difference equations, x(t + 1) = L(t) x(t), t = 0, 1, 2, . . . , (1) where the vector of population structure, x(t) ∈ R n + , belongs to the positive orthant of the n-dimensional vector space and L(t) is a nonnegative n × n matrix called the population projection matrix (PPM) [1,2]. Each component of x(t) is the (absolute or relative) number of individuals in the corresponding status-specific group at time moment t, while the elements of L, called vital rates [1], carry information about the rates of demographic processes in the population.
The pattern of (i.e., the allocation of zero/nonzero elements in) matrix L corresponds to the associated directed graph [3], which is called the life cycle graph [1] (LCG) as it represents graphically the biological knowledge of life histories involved into the model and the way the population structure is observed in the field or laboratory (Figure 1 gives an example). The LCG may be strongly connected [3], signifying a certain integrity of the individuals' life history and providing for the PPM being irreducible [4,5] (or indecomposable in the other In theoretical layouts and practical applications, it turned out to be useful to consider the PPM as the sum, of its parts responsible respectively for the transitions (T) between individual statuses and population recruitment (F) [1,8]. In particular, we have for the LCG shown in Figure 1.
When calibrated reliably from data, matrix L(t) gives rise to a rich repertoire of qualitative properties and quantitative indices characterizing the population under study at the place where and the time when the data were mined [1,2]. In particular, its dominant eigenvalue λ 1 (L) > 0, a positive eigenvalue coincident with ρ(L) (the spectral radius of L) and existing by the classical Perron-Frobenius theorem [4][5][6], gives a quantitative measure of how the local population is adapted to its environment [7], thus serving as an efficient tool of comparative demography [1] and enabling a forecast of population viability [9]. This ability ensues from the dynamics of trajectory x(t) as t tends to +∞ when (a primitive [5]) L(t) = L does not change with time: where x* is a positive eigenvector corresponding to λ 1 [6]. Therefore, we have so that the location of λ 1 relative to 1 is of principal importance. Of biological importance is one more characteristic, namely, the net reproductive rate which turned out to meet all of the following conditions [1,8,9]: where ⇔ means 'if and only if'. Matrix L belongs a set of PPMs, i.e., it is representable as the sum (2) where T is substochastic (none of the column sums exceeds 1) but not pure stochastic (each of the column sums equals 1). (6) suggests the convergence of the power series in T, it requires a stronger restriction on L than just T being substochastic. The convergence requires lim t→∞ T t = 0, which mathematically is equivalent to ρ(T) < 1 and biologically means nonzero mortality in the population.

Remark 1. Since Expression
Equivalence (7) was later called indicator property [10] (as it indicates the location of λ 1 relative to 1). Any function R(L) that meets conditions (7), or equivalently was called a potential-growth indicator (PGI), and an indicator R 1 (L) was proposed that was much simpler than R 0 (L) (6) [10]: where I denotes the identity matrix.
As an immediate consequence from the theory of rank-1 corrections of nonnegative matrices [11], the indicator property of R 1 (L) was proved for the class of PPMs in which the fertility matrix F has rank 1. A particular case is cited in (3), meaning a single stage where the population recruitment appears (Figure 1), or F may have only one nonzero column meaning a single reproductive stage in the life cycle ( Figure 2 and matrix (3) in [12]).
These two particular cases are, however, not exhaustive among PPMs, and since R 1 (L) has revealed a certain "merit of indication" [2] also motivated by the current modeling practice [13], the challenge is to expand R 1 (L) to a wider class of PPMs or to propose another PGI simple enough in calculation. Aimed at responding to this challenge, we present here a construction of a more general PGI by means of a general algorithm developed recently by V.N. Razzhevaikin and E.E. Tyrtyshnikov [14,15] for nonnegative matrices. We expose this algorithm in Section 2. Section 3 reports on a case study of Calamagrostis epigeios, a perennial herbaceous graminoid, a traditional object to study in the Russian botanical school [16], with the PPMs beyond the class of rank-1 F. The "Wider Merit of Indication" by PGIs are thereafter discussed in Section 4.

Razzhevaikin-Tyrtyshnikov's Algorithm to Construct a PGI
Note first that R 1 (L) (9) cannot be expanded for an L with rank(F) = 2, which is illustrated in Section 3.3. Therefore, the task to develop a more general PGI is really urgent, and a more general framework for that [15] is presented below.
The above notations are standard, and we introduce some original ones below.

Definition 1.
We call renumbering in a square matrix any permutation of its rows and the same columns.

Corollary 1.
Renumbering in a matrix does not affect its spectrum, hence the spectral radi us, too.
Matrix A being irreducible is then equivalent to the lack of a renumbering that might transform A to a block-triangular form with two or more nonempty diagonal blocks.

Definition 2.
A submatrix B of matrix A is called a principal submatrix if it is formed of the elements that are symmetric about the principal diagonal. There is a renumbering in A that places B into the lower right position. When the renumbering is needless, the principal submatrix is called the lower principal submatrix.
Definition 3. Let f 1 and f 2 be two functions defined on the same domain Ω and taking values in the spaces of real matrices sized n 1 × n 1 and n 2 × n 2 , respectively. We call f 1 and f 2 to be R-equivalent if the equality sign(ρ(f 1 (ω)) − 1) = sign(ρ(f 2 (ω)) − 1) holds true for any ω ∈ Ω.

Remark 2.
Domain Ω determines certain sets within the spaces defined by the elements of f 1 and f 2 . In particular, if those elements depend (explicitly or implicitly) on some parameters, then Ω can be regarded as a set in the space of the parameters.

Remark 3.
Definition 3 reduces to the PGI definition (7) when Ω is the domain where R(L) does possess the PGI property and n 2 = 1.
The following definition generalizes the PGI notion to a general set Ω.
Definition 4. Let n 2 = 1 and function f 2 (ω) defined on Ω be R-equivalent to f 1 . Then, the absolute value of f 2 (ω) (or f 2 (ω) itself when f 2 (ω) ≥ 0 for all ω ∈ Ω) is called a stability indicator of matrix f 1 over the Ω set.
When Ω is a cone of nonnegative matrices, Definition 4 reduces to the PGI notion (7) by the classical Perron-Frobenius theorem for nonnegative matrices [5].
The proof is given in [15], as well as those of the other theorems providing for the construction of a stability indicator for an arbitrary nonnegative matrix.
If the theorem conditions are met at each successive step of the dimensionreducing chain: then function f 1 = a 1 11 can serve as a stability indicator for matrix A ≥ 0. In what follows, we refer to this indicator as R RT (A).
The following obstacles may occur to perform the next step in Chain (12) when applied beyond the terms of Theorem 1.
(1) The next a k kk > 1. In this case, ρ(A k ) > 1, hence we have ρ(A n ) > 1 by virtue of Theorem 1. (2) The next a k kk = 1 and a k ii = 1 for some i < k. In this case, a proper renumbering in matrix A k relocates a k ii to the lower right angle and, if a k ii < 1, we can continue the chain, but if a k ii > 1, we come to point 1) above.
2) There exists a renumbering in A k that brings it to a triangular form with units on the principal diagonal. In this case, ρ(A k ) = 1, so that ρ(A n ) = 1 by virtue of Theorem 1.
When the terms of Theorem 1 are met at every step of Chain (12), the resulting indicator has the form of a multi-story fraction of the matrix elements. This fraction can be simplified by means of the following downgrade procedure: we replace the fraction gained at the current step in the form of α k−1 β k−1 with a new indicator in the form of α k−1 − β k−1 + 1. This procedure eventually enables an indicator in the form of a polynomial in the matrix elements to be obtained.
It turns out (see [15]) that this polynomial is equal, in the above case, to The aforementioned alternatives allow a similar consideration applied to the lower principal l × l submatrices B l of matrix A used to calculate the elements a k kk with k = n + l − 1 to be compared with 1 at the l-th step in the Chain (12). In the robust case, i.e., for a k kk > 1 and a j jj < 1 for j < k, the inequality 1 − p(1; B l ) > 1 holds true, while the inequalities a k kk < 1 holding for all k ≤ n imply that R RT (A) < 1. Here, we have where B runs the chain of n lower principal submatrices B j , j ≤ n, of matrix A (beginning with the matrix itself). Thus, Expression (13) is appropriate as a stability indicator in case the strict inequalities hold. Since renumbering in a matrix does not affect its spectrum, one can also use a chain of upper principal submatrices in Expression (13).
Problematic cases arising for a fixed numbering of rows and columns in the construction of a stability indicator of matrix A due to the equality a k kk = 1 holding true at the next step force us to use the same expression instead of Expression (13), albeit with B that runs through all possible principal submatrices of matrix A (there are 2 n − 1 such Bs, see [15] for greater detail).

Material and Method
Reed grasses of the Calamagrostis genus are perennial herbaceous polycarpic graminoids growing in a wide range of ecological conditions, mainly in forests and meadows, a traditional research object in the Russian school of population botany. Their ontogenies are well studied, including that of Calamagrostis epigeios [16][17][18], where the ontogenetic stages of individual plants are detectable from the morphology of the above-ground parts. Moreover, the chronological age (in years) of the plant can also be determined from the annual increments in its tillering zone [18].
As a pioneer species of plant succession, C. epigeios invades open plant communities, such as spruce clear-cut areas, due to extensive population growth by means of vegetative propagation. However, the mosaic of its colonies (systems of clonal plants) has not yet formed a continuous cover during the initial years of colonization, so the colonies have distinguishable boundaries [16]. Annual observations (late August, when the herbs complete their development) show that different plants spend different numbers of years at the same stage of ontogeny, i.e., there are individuals of different ages among the plants at the same stage. The diversity of individual pathways of development among the plants within a single colony can be represented in the form of LCGs shown in Figure 2. These graphs are constructed on a two-dimensional integer-valued "lattice" of all possible states that an individual plant may have in terms of age-stage. Botanists call this diversity polyvariant ontogeny and consider it as a major mechanism of plant adaptation to the environment [16,[18][19][20].
The numerical values of ontogenetic transition rates assigned to the solid arrows in Figure 2 were obtained from the data of "identified individuals" type [1], p. 134, while the reproduction rates are shown on the dashed arcs of the LCG as symbolic parameters to be calibrated rather as a tribute to the tradition of reproductive uncertainty, i.e., inability to detect the status of parents for each recruited plant, hence inability to calibrate the status-specific reproduction rates in a unique way [19]. In fact, the 2017 project aimed to eliminate the reproductive uncertainty inherent in the analysis of above-ground vegetation [ibidem], and for this purpose the whole reed colony (two of them in fact) was dug up, while preserving the entire system of rhizome parent-offspring links. The fragments of the linkage system related to each parent plant (Figure 2 in [16]) were analyzed and summarized over the colony, thus resulting in a total picture of parent-offspring links in a colony. Figure 3 represents the diagrams of offspring survival for all the age-stage-specific parental groups in each C. epigeios colony [16]. The integer-valued parameters a, b, ..., n, o are the numerators of the reproduction rates that are shown on the arcs of the LCG (Figure 2), and the denominators are the (absolute) sizes of the corresponding parental groups in 2014, i.e., one year before the excavation ( Table 2 in [16]). Thus, the vector-matrix Equation (1) takes on the form of

Calibration of the C. epigeios PPMs
which is satisfied for each of 11 (for Colony I) or 14 (for Colony II) components of vector x(t) with matrix L presented in Table 1. Unlike Equation (1), which describes the "projection" of the population structure one step into the future, Equation (14) describes the transition to the current state of the colony from its last year's state, retrospectively restored based on the excavation data. However, this circumstance does not affect the properties of the model, and in particular, the dominant eigenvalue, λ 1 (L rec ), of the reproductive-core submatrix of matrix L (outlined by dashed lines in Table 1, with deletion of unnecessary rows and columns) still serves as a quantitative measure of the colony's adaptability to environmental conditions at a given time step. For comparison purposes, this measure is modified by a "scaling factor" S 0 /S 1 accounting for the expansion of colony area in one year [16].   Figure 3), the corresponding λ 1 (L)s *, and the refined adaptation measures µ(L) ** (according to Figure 2I). Calibrated without any reproductive uncertainty, the C. epigeios PPMs give rise to population characteristics of certain biological interest [16]. In particular, the values of λ 1 much greater than 1 (Table 1) witness the high speed of population growth in the young colonies.

Parent
As we can see from Table 1, the fertility matrices have ranks 2 and 3 for Colonies I and II respectively, and we therefore cannot apply the PGI R 1 (L). Expression (13) for the R RT (L) produce explicit formulae for Colonies I and II (Appendix A, (A6), (A7)), which result in the numeric values shown in Table 1 just to illustrate the PGI property.

Calibration from the Above-Ground Data Alone
The excavation of colonies is a much more laborious method of field study than the above-ground observations of marked individuals on permanent sample plots bearing reproductive uncertainty [19]. Though revealing valuable information on the clonal rhizome system, the excavation destroys the field object, thus ceasing its further study. Therefore, developing reliable methods to cope with reproductive uncertainty is still an urgent challenge [19,[21][22][23]. The question is how the PGI (if found in a relevant form) could help respond this challenge.
It follows from Figure 2 that the reproductive-core submatrices for Colonies I and II are, respectively, and Although the ontogenetic-transition matrices T 1 and T 2 were calibrated from the excavation data [16], they could also be calibrated, in a unique way, from the above-ground data of the "identified individuals" type [1,19]. These data are, however, unable to determine the status-specific reproduction rates (the elements of F), leaving them as uncertain parameters, which have only to satisfy certain linear constraints. These constraints express the observation that population recruitment at a status-specific group consists of the contributions by parental plants at several generative statuses, i.e., for Colony I and    a + b + c + d + e + f = 207, for Colony II. It is known that during the vegetative propagation of long-rhizome grasses, the young rhizomes of the first-year virginal plants develop most actively [24], and this activity decreases over the years. The following hierarchy of inequalities expresses this expert knowledge for Colony I: where C( . . . ) denotes the total number of alive daughter plants produced by the parents of this status-specific group. A similar hierarchy for Colony II looks as follows: In terms of our integer-valued parameters, these hierarchies take on the following forms: for Colony I, and for Colony II. Now, the basic Equation (14) is satisfied with the 8-xparameter family of PPMs constrained by (17), (21) for Colony I, and with the 13-parameter family constrained by (18), (22) for Colony II.
The integer-valued feature of these families lead us to Diophantine systems, which may have a lot of solutions in our case, following the number of ways to split a given sum into a given number of summands (diminished by the additional constraints) multiplied with similar numbers for the second given sum and for the third. However, the principal point is that this number is finite (no matter how great), which means that all these solutions can be found by enumerating all the feasible options. Every solution thus found, i.e., a set of nonnegative integer numbers a, b, ..., n, o, satisfying the system of equations and inequalities, forms its own reproductive-core submatrix L rec (as shown in Table 1) with the corresponding value of λ 1 (L rec ), and we find those ones among (a finite number of) such matrices that give λ 1min and λ 1max , the minimum and maximum value of λ 1 , respectively. Thus, we obtain a certain range [λ 1min , λ 1max ] for uncertain λ 1 (L rec ), or the fitness bounds (lower and upper) for the local population of a given clonal species with polyvariant ontogeny ( Table 2).   1 Multiplied with a scaling factor S 0 /S 1 equal to 0.7109 [16]. 2 R 0 (L rec ) for the minλ 1 set of parameters. 3 R 0 (L rec ) for the maxλ 1 set of parameters. 4 R 1 (L rec ) for the minλ 1 set of parameters. 5 R 1 (L rec ) for the maxλ 1 set of parameters. 6 Occasional coincidence with 4) within 4 decimal digits; further digits reveal a difference. 7 R RT (L rec ) for the minλ 1 set of parameters. 8 R RT (L rec ) for the maxλ 1 set of parameters. 9 Parameters in the next right column. 10 Parameters in the second right column.
To illustrate the relationships of λ 1 s with various potential-growth indicators, Table 2 also shows the corresponding PGI values. In particular, we see that 1 < λ 1 (L rec ) < R 0 (L rec ) in both Colonies, thus illustrating Theorem 3.3 of Li and Schneider [9]. The negative values of R 1 demonstrate the PGI inability of R 1 (L) when L has a rank-2 or rank-3 fertility matrix F. The two columns on the extreme right show that λ 1 and R RT attain their extreme values at quite different parameter points.

Discussion
When models look complicated, such as those shown in Figures 2 and 3, the additional attributes such as PGIs are of logical use, rather than attempting to simplify the model, for instance, by aggregating several variables. However, the complexities shown in Figures 2 and 3 reflect the reality of empirical data, while simplifications may distort the real picture.
For example, although the LCG shown in Figure 1 is fairly simple, it bears what we call "reproductive uncertainty" [25], p. 43, i.e., the inability to calibrate the status-specific reproduction rates, a and b, from field data because only the total contribution to the population recruitment by the two generative groups, g and gt, is observed in the field. As a result, we could only obtain a certain range of uncertain λ 1 values rather than a particular number [26].
Among various, poorly efficient, ways to eliminate reproductive uncertainty, there was an attempt to aggregate the two generative groups into one, thus eliminating the very basis of uncertainty and providing for the PPM calibration in a unique way [27]. A logical step was to obtain the λ 1 of the aggregated model within the range of the original one, and this happened indeed for several observation years. However, there were also a pair of years (2013, 2014; Table 3, ibidem) that gave the aggregated λ 1 outside the range. Moreover, it was greater than 1, whereas the range was located entirely to the left of 1 (Table 2, ibidem), i.e., the aggregation resulted in a principally wrong forecast of population viability.
Another mode of eliminating reproductive uncertainty was to accept an optimality hypothesis: the uncertain reproductive contributions by status-specific parental groups are distributed in a way that maximizes the resulting λ 1 [2,22,25]. The calibration problem was reduced, then, to a nonlinear constraint maximization problem over a polygon of variables to be optimized, while the indication by R 1 (L) enabled, using the corresponding linear-programming problem, all the pertinent vertices of the constraining polygon to be enumerated [ibidem]. Although the question remained whether or not biological evolution leads to a maximum population growth rate, or more generally, whether it obeys any optimization principle [28,29], the hypothesis did allow a single tool to be used in comparative demography of the species with either vegetative of seed mode of reproduction [2,22,25].
The laborious colony excavation method presented in Section 3.2 does eliminate the reproductive uncertainty but destroys the object of study, while the alternative method of integer-valued calibration (Section 3.3) turns out to be resource-consuming, too. Here, we expected another merit of indication for saving computer time when solving the Diophantine systems (17), (21) and (18), (22) by means of enumeration. Using a computer algorithm to find all the solutions for Colony I (eight unknowns, Appendix A) poses no problem, but for Colony II (13 unknowns, Supplementary Materials), it does. The corresponding 13 nested cycles of the algorithm contain 57 × 14 × 29 × 9 × 4 × 24 × 50 × 25 × 9 × 4 × 24 × 6 × 6 = 7.7739 × 10 14 repetitions of the inmost check, and a crude estimate suggested 14 to 24 days of continuous operation with a 2.3 GHz Intel Core i9 processor. Since the Colony was growing fast, there was an idea to conduct a preliminary check for R RT (L 2rec ) > 1 and exit otherwise in order to accelerate the inmost cycle. However, the idea failed once we found that R RT (L 2rec ) > 1, hence λ 1 (L 2rec ) > 1, already at the lower bounds of all the parameter ranges (see (A7) in Appendix A). Since λ 1 is a monotone function of matrix elements [5], it follows that λ 1 (L 2rec ) > 1, hence R RT (L 2rec ) > 1 over the total ranges of parameters, and this renders the check useless in each particular case.
However, the idea was fruitful, in this regard, when the task was to investigate the effects that uncertain parameters of the seed bank exert on the model outcome and to find the bounds of uncertainty in a case study of Androsace albana, another alpine short-lived perennial [13]. The R 1 (L) turned out to be linear with regard to the uncertain parameters, and this linearity caused a great "merit of indication" while checking whether λ 1 (L) was less or greater than 1.
Similarly, the calculation of R RT (A) turned out to be fruitful in the evolutionary optimality problems [30] where ρ (L) and R RT (L) were considered as selection criteria and an additional condition, ρ (L) = R RT (L) = 1 held true. The condition guaranteed the optimal values being reached at the same sets of parameters-unlike the extremal points in our current problem (Table 2)-thus reducing the optimal search to that for a polynomial function.
In addition, the use of simple algorithms was highly efficient to solve stability problems for highly sparse nonnegative matrices [31].

Conclusions
The history of PGIs in matrix population models, which began with the simplest expression for the Leslie matrix three quarters of a century ago [32], has now come to its logical end with the explicit formulas for R 0 (L), R 1 (L), and R RT (L), where L is a PPM representable as L = T + F (2). The theme has been developing along two directions, the first direction tending to simplicity, the second one to generality.
While R 0 (L) (6), the net reproductive rate, confines only to the convergence of a series in the definition (6), hence ρ (T) < 1, thus leading (till 2020) in the generality, it has been loosing in simplicity as the spectral radius represents a much more difficult computational task than the determinant in R 1 (L) (9) or a polynomial of matrix elements in R RT (L) (cf. (A4)-(A6)).
The dimension-reducing algorithm by Razzhevaikin and Tyrtyshnikov [14,15] has reached the peak of generality, as it can be applied to any square nonnegative matrix. Although, it fails in the monotone PGI property (when inequality R(L 1 ) > R(L 2 ) being equivalent to λ 1 (L 1 ) > λ 1 (L 2 )), so that the extremal points over a compact set are different for R(L) and λ 1 (L) (see, e.g., Table 2), the indicators can still be useful for solving constraint optimization problems [25].
To conclude, the current theory of PGIs offers three explicit functions-R 0 (L), R 1 (L), and R RT (L)-of matrix elements, which complement the mathematical toolbox of matrix population models. At certain levels of generality and simplicity, these PGIs are able to cover a vast variety of particular situations. Particular problem formulation, be it a calibration task under parameter uncertainty, or stability study in a certain class of PPMs, or anything else with nonnegative matrices, should prompt a pertinent PGI to obtain "more merits of indication" in solving the problem.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
Both the Expressions (A4) and (A5) are obviously much more complicated than (A1), at least because of square-root terms. It is also more complicated than any component of the (A6) We have eliminated a pair of original rows as some current ones evidently majorized them.
The expression for R RT (L 2rec ) as a polynomial function of 13 parameters is similar to (A6) but cumbersome enough (Supplementary Materials) to be not reproduced in a publication. However, its value at the lower bounds of the parameter ranges is R RT (L 2rec ) = 5.0870 > 1, hence λ 1 (L 2rec ) > 1 due to the PGI property, and it is so everywhere within the ranges due the monotone property of λ 1 [5].