Matrix Models for Size-Structured Populations: Unrealistic Fast Growth or Simply Diffusion?

Matrix population models are widely used to study population dynamics but have been criticized because their outputs are sensitive to the dimension of the matrix (or, equivalently, to the class width). This sensitivity is concerning for the population growth rate () because this is an intrinsic characteristic of the population that should not depend on the model specification. It has been suggested that the sensitivity of to matrix dimension was linked to the existence of fast pathways (i.e. the fraction of individuals that systematically move up a class), whose proportion increases when class width increases. We showed that for matrix population models with growth transition only from class to class , was independent of the class width when the mortality and the recruitment rates were constant, irrespective of the growth rate. We also showed that if there were indeed fast pathways, there were also in about the same proportion slow pathways (i.e. the fraction of individuals that systematically remained in the same class), and that they jointly act as a diffusion process (where diffusion here is the movement in size of an individual whose size increments are random according to a normal distribution with mean zero). For 53 tree species from a tropical rain forest in the Central African Republic, the diffusion resulting from common matrix dimensions was much stronger than would be realistic. Yet, the sensitivity of to matrix dimension for a class width in the range 1–10 cm was small, much smaller than the sampling uncertainty on the value of . Moreover, could either increase or decrease when class width increased depending on the species. Overall, even if the class width should be kept small enough to limit diffusion, it had little impact on the estimate of for tree species.


Introduction
To model the dynamics of a population, there are two main options depending on the level of the population: individual-based models, where the trajectory of every individual is monitored; and distribution-based population models, where individual attributes are summarized by their population-level distribution [1,2]. Among the latter, four types of models can be distinguished depending on whether the distribution and time are modeled as continuous or discrete: matrix population models (discrete distribution and discrete time [3]), integral projection model (continuous distribution and discrete time [4]), continuous-time Markov chain (discrete distribution and continuous time, e.g. [5]), and partial differential equations (continuous distribution and continuous time [6]). There are two diverging opinions on the choice of a modeling approach. On one hand, some have highlighted the advantages and limitations of each approach, thus suggesting that some approach may intrinsically be superior to the others. For instance, matrix population models have been criticized for the arbitrariness of the class division [7][8][9][10] and integral projection model put forward as a solution to this issue [11][12][13]. On the other hand, others have put the emphasis on the theoretical connections that exist between all these modeling approaches [14][15][16][17][18], thus suggesting that the choice of a modeling approach should be a pragmatic choice that marginally affects the predictions [3,18.5].
Matrix models have been criticized for their inability ''to incorporate variation among individuals within a size category'' [7, p.346]. Because matrix models operate on size-(or age-, or stage-) structured populations, differences of growth among individuals due to size are accounted by the model, but size may indeed be an incomplete predictor of growth. Nevertheless, in this case, additional predictors of growth can be added as structuring variables of the population, with a subdivision of the categories of the transition matrix [3, 18.4] [19]. To address autocorrelation in growth, second-(or higher-) order Markov chains (and corresponding transition matrices) can also be considered [20]. Finally, residual error in growth which results from random variability in individual growth can also be addressed in matrix modeling using random shocks [21] or as a diffusion process (i.e. by adding transition rates off the main diagonal of the transition matrix [3, p.199]).
Second, matrix models have been criticized because their outputs (population growth rate l, i.e. the temporal rate of change of the population number of individuals on the long term; elasticities, i.e. the relative rate of change of an output with respect to a parameter; age estimates) are sensitive to the dimension of the matrix (or, equivalently, to the width of the size classes for a given range of size) [8,9,13,22]. This sensitivity to matrix dimensionality is concerning when the outputs are intrinsic population characteristics that should be defined irrespective of the mathematical model of population dynamics and, thus, irrespective of the particular division of size (or age) into classes. The dependence of elasticities on matrix dimensionality is consistent with the fact that the relative importance of growth compared with stasis (i.e. remaining in a size-class for more than one time step) changes with the class width [23]. Moreover, the dependence of the population growth rate l on matrix dimensionality has been questioned by other studies [24,25]. In some cases, the influence of matrix dimensionality on model outputs has been investigated using different populations with different models fitted to them [10,13,26]. In this case, the effect of matrix dimension on model outputs can be confounded with the effect of population differences.
As a possible explanation of the dependence of l on matrix dimensionality, Zuidema et al. wrote [7, p.346]: ''As the transition probabilities in a matrix model depend only on the current situation, there is no obstruction for unrealistically fast pathways through the life cycle. For instance, in matrix models with 10-cmwide diameter categories and small progression probabilities, a small fraction may reach 50 cm diameter in five time steps, something that is clearly impossible biologically (and physically). This fraction contributes strongly to population growth and probably causes the high estimates of l for small matrix models.'' Because the population growth rate l is an intrinsic characteristic of the population that is often used in population viability analysis and that should not depend on the model specifications, understanding why l depends on matrix dimensionality in matrix population models would indicate to which extent different modeling approaches are reconcilable. In case of significant dependence, it should provide guidance to the modelers about which modeling approach to use.
In this study, we will assess the dependence of the population growth rate l on class width, and how fast pathways (i.e. the fraction of individuals that systematically move up a class at each time step) possibly contribute to this dependence. We will show that if there are indeed fast pathways, there are also in about the same proportion slow pathways, i.e. the fraction of individuals that remain in the same class for a long time. As a consequence, fast pathways do not directly bias growth rates towards higher than expected values, but rather, in combination with slow pathways, act as a diffusion process with limited impact on the estimate of l. Diffusion here is defined as the movement in size of an individual whose size increments are random following a normal distribution with mean zero. At the population level, these individual random walks flatten the size distribution and make it uniform (in the statistical distribution sense). This study is based on a data set from a tropical rain forest in central Africa, including different tree species whose dynamics are modeled by a Usher transition matrix (i.e. a transition matrix with non null elements on its main diagonal, on its lower subdiagonal and on its first row only).

The M'Baïki Forest
The M'Baïki experimental site is located in the south of the Central African Republic (3 0 549N, 17 0 569E), at the northern limit of the rain forest of the Congo basin. It is dedicated to studying the effects of logging damage on stock recovery [27] and lies in a terra firme rain forest. The experimental design of the site consists of two blocks of three and one block of four 300|300-m permanent sample plots with a 50-m inner buffer zone. In each central 200|200-m square, all trees over 10-cm diameter at breast height (dbh) were identified and georeferenced. Since 1982, girth at breast height, standing deaths, treefalls, and newly recruited trees over 10-cm dbh have been monitored annually except in 1997, 1999and 2001. Between 1984and 1985, two silvicultural treatments were applied: three plots (including the buffer zone) were logged, and four plots were logged and thinned. The three remaining plots were left as controls. For this study, we used the data of the control plots only, between 1982 and 2006 with a time step of 2 years. In total, 66,749 tree records in 53 species were used in this study.

Changing the Dimension of a Transition Matrix
Given a transition matrix U for a size-structured population with K size classes, different techniques have been proposed to derive a transition matrix U 0 for K'vK size classes. When the data used to fit U are available, a data-driven approach consists in refitting the transition matrix using these data and K' classes. This approach has the advantage that the possible resulting change in the estimate of the population growth rate l readily corresponds to what is observed when fitting the matrix model. The limitation is that it is not possible to disentangle what is specifically due to the data set used and what is due to the properties of matrix modeling in general.
A second technique consists in computing the transition rates of the U 0 matrix from those of U. For instance, when K'~K=2 and when combining every two successive classes i and iz1 into a single one j, it has been proposed to compute the upgrowth transition rate p' j , the mortality rate m' j and the recruitment rate r' j of U K' as [8,9]: where p i , m i and r i are the transition rates of U K for growth, mortality and recruitment respectively, f i is the number of individuals in size class iƒK, and f ' j~fi zf iz1 is the number of individuals in size class jƒK'. This approach has the advantage that no additional information beyond the population-level characteristics are needed to change the size of the transition matrix. A first limitation is that this technique can be used only when U 0 is nested into U (i.e. when the K' classes are obtained by merging together some of the K classes). A second limitation, more theoretical, is that the relationships between the transition rates of U 0 and those of U depend on the number f i of trees in the classes. Because f i changes with time, this implies that the U 0 matrix derived from U will not be the same depending on the time step considered (even if U is stationary), which is not consistent. Moreover, other algebraic relationships than (1) could be used to collapse U into U 0 . In particular, it is mathematically feasible to collapse a transition matrix into a smaller matrix while maintaining the same dominant eigenvalue and eigenvectors [10,25,[28][29][30]. If such algebraic relationships were to be used rather than (1) to collapse matrices (and there is no theoretical reason for not doing it) then, by construction, l would not depend on the class width.
A third technique to change the dimension of the transition matrix dates back to [14,15] and is based on the connection between matrix population model and continuous partial differential equations [3,18.1.4]. In particular, it has been the basis for optimizing the width of size classes in matrix models for sizestructured populations [31,32], considering that the matrix model is a discrete approximation of a continuous partial differential equation with a bias/variance trade-off to optimize. It has the advantage that it can deal with any change of the class limits in a theoretically consistent manner [33]. It has the limitation that size must be continuous and that a model of continuous-size dynamics must be assumed as a prerequisite. The matrix population model can be seen as a discretization of a McKendrick partial differential equation [16,[34][35][36]: with the boundary condition:   where f (x,t) is the continuous size distribution at time t, such that the number of individuals with a size between x and xzd for any infinitesimally small d is f (x,t)d, a is the size growth rate, m is the mortality rate, r is the recruitment rate, and x 0 is the minimum size for inventory. Let us consider a numerical scheme to solve (2) [37, chapter 20]. Size is discretized using a size spacing of d: x i~x0 z(i{1)d where x 0 is the minimum size and i~1,…, I. Time is discretized using a time spacing of t: which can be written as: where F(n)~½f n 1 ,…, f n I is the vector of length I that contains the number of individuals at discrete time n in each size class with width d and lower bound x 0 z(i{1)d, and U is a Usher transition matrix: where the stasis rate q i (i.e. the probability for an individual to stay alive in class i between two consecutive time steps), the upgrowth transition rate p i (i.e. the probability for an individual to grow up from class i to class iz1), and the recruitment rate g i are given by: Equation (5) corresponds to the Usher matrix model for sizestructured populations [38,39], which has developed in forestry independently from the McKendrick equation (e.g. [40][41][42]). Given an individual-based model of size growth a(x), an individual-based model of death probability m(x), and an individual-based model of recruitment r(x), equation (7) defines the transition rates of the Usher matrix U for any partition of size into classes of width d. Notice that p i ƒ1 for all i implies: which corresponds to the Usher [38] hypothesis that individuals cannot grow by more than a single class in one time step.

Population Growth Rate
The asymptotic population growth rate, l, is the dominant eigenvalue of the transition matrix U as given by (6) [3]. In the general case, there is no explicit expression for l. However, in the particular case of a Usher transition matrix with constant mortality rate m and constant recruitment rate r (i.e. m(x) and r(x) do not depend on size x), and irrespective of the variations of the growth rate a(x) with size, there is an analytical expression for l [24,43,44]: l~1{mzr ð9Þ which mathematically proves that l does not depend on the class width in this case. This result is also valid when transition matrices are collapsed using (1), because m i :m and r i :r for all i imply m' j~m and r' j~r for all j. A corollary of this result is that variations of l with class width can occur only if the mortality rate or the recruitment rate varies with size. Therefore, in this study, we will only consider matrix models with either size-dependent mortality or size-dependent recruitment.
When dealing with forest ecosystems, due to the complexity of sexual and asexual reproductions and variability of elapsed time for germinated seeds to become recruited trees, it is not possible to assign a newly recruited tree as originating from a given size class [13]. Therefore, when dealing with forest dynamics, an average recruitment rate (the same for all size classes) is generally estimated as the ratio of the number of newly recruited trees over the number of trees at the previous time step (see [22,45] for exceptions). If such is the case, variations of l with class width can occur only if the mortality rate varies with size.
Practically, to assess how l varied with class width for the tree species at M'Baïki, the following analyses were performed. Trees were classified into diameter at breast height (dbh) classes with equal width d, ranging from a minimum dbh for inventory of 10 cm to a maximum dbh of 150 cm. The number K of dbh classes correspondingly varied proportionally to 1=d. The time step of the matrix model was t~1 year. For each species with at least 300 observations, a constant dbh growth rate a was estimated from the M'Baïki data base as the empirical mean of the dbh increments (including negative increments) over 2 years divided by this period of 2 years. The variations of the growth rate with dbh were not considered here because they are not a condition for l to vary with class width d. A constant recruitment rate was estimated for each species as the ratio of the number of newly recruited at year t over the number of living trees at year t{2, divided by this period of 2 years.
The dependence of mortality on dbh was necessarily accounted because it is a condition for l to vary with class width d. The tree mortality rate was modeled for each species as a function of tree dbh using one of the three following models [46,47]: where logit {1 (z)~(1z exp ({x)) {1 is the inverse logit function and a, b, c are parameters to estimate. Models (10) and (11) were fitted using the generalized linear model (command glm in R software) whereas model (12) was fitted using the generalized nonlinear model (package gnm in R software). The three models were compared using the Akaike Information Criterion (AIC) and the one with the lowest AIC was retained. Given the growth rate a, the mortality rate m and the recruitment rate r for each species, the class width d was changed from d min~m axf1, tag cm to d max~1 0 cm; for each value of d, the Usher transition matrix was computed using (7), and the population growth rate l was computed as the dominant eigenvalue of this matrix. The variations of l with d were compared to the sampling variability of l for the smallest class width d min (i.e. for the matrix model that is the closest to the McKendrick equation (2)). The sampling variability of l represents the uncertainty on l due to the finiteness of the data set used to estimate l. For each species, a 95% confidence interval of the estimate of l for the smallest class width d min was computed using 500 bootstrap replicates [48,49].

Fast and Slow Pathways
The McKendrick equation corresponds to a propagation of the diameter distribution at a speed defined by a(x). Hence, implicitly with the McKendrick equation, all trees with the same diameter grow at the same rate and there is no fast pathway. If the Usher matrix model was an exact scheme to solve the McKendrick equation, there would be no fast pathway either. Therefore, the fast pathways depicted by Zuidema et al. [7] correspond to the approximation brought by the discretization of (2) into (5).
When focusing on the upgrowth part of the McKendrick equation (i.e. setting the mortality and the recruitment rates to zero), a von Neumann stability analysis [37,120.1.1] shows that the Usher matrix model is a stable numerical scheme to solve the McKendrick equation provided that condition (8) is met. Therefore, the Usher condition that no individual can grow by more than one class in a single time step identifies with the Courant-Friedrichs-Lewy condition in numerical analysis. The von Neumann stability analysis also shows that the Usher matrix model is numerically dissipative unless a(x i )t=d~1 for all i. This means that in the Fourier transform of the distribution f , all terms with a wave number k such that kd&1 will be inaccurately calculated. This dissipative effect of the Usher scheme can be intuitively understood by considering that the numerical scheme (4) can be rewritten as: where b 2 i~a (x i )d for all i. The numeric scheme (13) corresponds to the forward-time centered-size differencing scheme for the Fokker-Planck equation: with b(x) 2~a (x)|d. When compared to the McKendrick equation (2), the Fokker-Planck equation has an additional term (the second-order partial derivative in (14)) that corresponds to diffusion. As a first result, the fast pathways pointed out by [7] can be interpreted as a diffusion process and go along with slow pathways, i.e. individuals that remain in the same size class longer than expected. The diffusion process has a simple biological interpretation and may be a desirable feature [40,42]. It relates to the individual variability in growth. More precisely, b(x) 2 t can be interpreted as the variance of the individual size increments during an infinitesimally small time interval t for trees with size x. Hence, the diffusion process will be realistic provided that b(x) is a realistic model for the standard deviation of annual tree size increments. The Usher scheme implies that b(x) 2~a (x)|d. Hence, the diffusion generated by the Usher scheme will remain biologically realistic as long as where s(x) is the standard deviation of tree size increments during the time interval of trees t years.
To visualize fast and slow pathways, we considered the transient dynamics of a single even-aged cohort of trees uniformly distributed between x 0 and x 0 zD where D is the dbh amplitude of the cohort, i.e.: f (x,0)~(N=D) 1(x 0 ƒxvx 0 zD), where N is the initial number of trees and 1(p) is the indicator function of proposition p (~1 is p is true and 0 if p is false). To focus on the transient dynamics of this cohort, recruitment was set to zero (r:0). In the particular case where the growth rate a is constant, the analytical solution of the McKendrick equation (2) is known [16, p.45] [36] and corresponds at time t to a displacement of the cohort by a dbh at with an attenuation of the number of trees with dbh x by exp ({ Ð t 0 m(xzas)ds), i.e.: In the particular case where m is given by model (10), By comparing the exact solution (16) of the McKendrick equation to the prediction of the Usher matrix model, fast and slow pathways due to discretization can be identified. We calculated the proportion of fast pathways as the proportion of F(t=t) that was above x 0 zatzD, and the proportion of slow pathways as the proportion of F(t=t) that was below x 0 zat.

Variations for Celtis zenkeri
To illustrate the dependence of the population growth rate l on the class width, we first describe the variations of l for just one species taken from the Mbaki data set. We chose Celtis zenkeri Engl. (Ulmaceae) because it was the most abundant species of the data set (with 7295 observations between 1982 and 2006) and representative of the most commonly observed pattern of variation of l. Its average growth rate was 0.263 cm yr {1 (standard deviation: 0.30 cm yr {1 ). Its recruitment rate was 1.028% yr {1 . Its mortality rate was best modeled by model (10) with a~{5:348 (std. dev.: 0.272) and b~0:026 (std. dev.: 0.009). Therefore, the mortality rate for C. zenkeri was an increasing function of tree dbh, ranging from 0.6% yr {1 for a dbh of 10 cm to 2.3% yr {1 for a dbh of 62 cm (that is the 99.5% percentile of dbh for C. zenkeri at M'Bïaki).
To illustrate fast and slow pathways, we considered the dynamics of an even-aged cohort of 100 trees uniformly distributed between 10 and 15.0 cm at t~0. When d~0:263 cm and t~1 yr, the condition at=d~1 was met and the Usher scheme for growth was not dissipative ( Figure 0A): there were neither fast pathways nor slow ones in this case. With the exception of the last class, non-null transition rates defined one-toone connections between classes. The only difference between the exact solution of the McKendrick equation and the Usher model followed from the difference between Ð t 0 exp ({m(xzas))ds (for the exact solution) and P t=t i~1 (1{m(x i )t) (for the Usher model) for the attenuation of the number of trees. For C. zenkeri, this difference was actually so small that it is not visually perceptible in Figure 1A.
When d~1:0 cm and t~1 yr, at=dv1 and the Usher scheme became dissipative ( Figure 1B): there were some fast and slow pathways in this case. Several classes at initial time contributed to the number of trees in any class at final time. For C. zenkeri, the proportion of slow pathways (22.8%) was greater than the proportion of fast pathways (15.9%). If the mortality rate m was constant, then the proportion of slow pathways would have been exactly equal to that of fast pathways. The class width d~1:0 cm brought the same diffusion as a random growth with standard deviation (ad=t) 0:5~0 :51 cm yr {1 , which is greater than the observed standard deviation of the growth rate (0.30 cm yr {1 ). The dissipative effect of the Usher scheme increased as the class width increased from 1.0 to 2.5 cm ( Figure 1C). For this latter class width, the proportions of slow and fast pathways were 30.9% and 20.1%, respectively, and the dissipation was equivalent with that produced by a random growth with standard deviation 0.81 cm yr {1 .
The population growth rate of C. zenkeri increased from l~1:00451 for a class width of 1 cm to 1.00457 for a class width of 10 cm. In comparison, the 95% confidence interval of the estimate of l for d~1 cm was 1.00212-1.00709. Therefore, the amplitude of the 95% confidence interval of the estimate of l for d~1 cm was 96 times greater than that of the variations of l for d varying from 1 to 10 cm (Figure 2), even though C. zenkeri was the species with the largest number of observations and the narrowest 95% confidence interval.

Variations Across Species
After excluding those species for which the mortality rate did not significantly vary with dbh (and thus with a population growth rate l that did vary with d), there were 53 species left. Models (10), (11) and (12) for mortality were selected for 85%, 9% and 6% of the species, respectively (Table S1). For 74% of the species, the population growth rate l increased with class width d; for 24% of the species, l decreased with d; and for 2% of the species, the change of l when d varied from 1 to 10 cm was less than 10 {6 . There was no one-to-one relationship between the direction of change of l and the shape of the mortality model, with all combinations of mortality model and direction of change of l being observed. Nevertheless, when the mortality rate m(x) was an increasing function of dbh x, l most often (but not always) increased with d.
On average across species, the amplitude of the 95% confidence interval of the estimate of l for d~1 cm was 31 times greater than that of the variations of l for d varying from 1 to 10 cm (Figure 3). No species had a population growth rate l for 1ƒdƒ10 cm that went outside the 95% confidence interval of the estimate of l for d~1 cm.
The relationship across species between the variance s 2 of the growth rate and the one-year dbh increment at could be modeled by a power relationship: s 2~0 :775(at) 1:434 (Figure 4). Combining the Usher/Courant-Friedrichs-Lewy condition (8), condition (15) and this power relationship gives the following approximate interval for the class width: atƒd * v 0:775(at) 0:434 Only very fast growing species, with a dbh growth rate greater than 0.64 cm yr {1 , cannot meet this condition. At M'Baïki, only four species (Musanga cecropioides R. Br., Trilepisium madagascariense DC., Macaranga paxii, and Ricinodendron heudelotii) had a mean growth rate greater than 0.64 cm yr {1 . For all other species (92% of the species), it would be possible to set a class width that is consistent with the variability in growth of the species. However, the resulting class width (less than 0.64 cm for a time step of 1 year) would be much less than the class widths commonly used in matrix modeling for forest dynamics.

Is l Sensitive to Class Width?
The Usher matrix model can be seen as a discrete approximation of a continuous-size distribution model. This discretization induces a diffusion, with fast pathways (i.e. fractions of trees that grow up classes faster than expected), but also slow pathways (i.e. fractions of trees that remain in the same class longer than expected). Diffusion in itself is appropriate since it corresponds to the individual variability in growth. Hence, slow pathways represent the fraction of individuals with the lowest growth (possibly including those with negative growth). There are also instances when fast pathways are appropriate, in particular for stage-structured populations when individuals are able to skip intermediate stages in their ontogenic development [50]. However, in size-structured populations, the strength of diffusion is directly related to the class width, and the class widths often used in matrix modeling in forestry (often in the range 3-10 cm for dbh; [51]) induces a diffusion that is much stronger than that solely due to the individual variability in growth.
Although the diffusion due to the discretization in classes was much greater than what would be realistic, the matrix model predictions of l at M'Baïki were particularly robust to changes of the class width d, for variations of d as large as from 1 to 10 cm. The median of the difference Dl between the maximum value of the population growth rate and its minimum for a class width ranging from 1 to 10 cm was 1:6|10 {4 at M'Baïki. In comparison, Enright et al. [8] found a difference Dl of zero for three tree and a grass species (including an imaginary tree species), and a difference of 10 {3 for the tropical conifer Araucaria cunninghamii. Ramula and Lehtilä [9] found differences Dlƒ0:01 for 19 tree species and significantly larger differences for 18 herbaceous species. Zuidema [22] found a maximum difference of 0.009 for four tree species and 2:5ƒdƒ10 cm.
Higher variations of l with class width may be obtained when the recruitment rates are not constant with size [8,9]. Ramula and Lehtilä [9] reported that changes in l were significantly larger for herbaceous than for woody species, maybe because the former have much less classes that are based on development stage rather than on size. Transition matrices do not have the same structure for herbs and trees, with higher recruitment for herbs and more frequent retrogression [52]. Therefore, the robustness of l to class width that was observed at M'Baïki for tree species is also linked to the specific structure of the transition matrix and may not be extrapolated to other types of transition matrices.
The direction of variation of l with the class width depended in a complex and non-systematic way on the demographic rates (growth rate, mortality rate, recruitment rate). The population growth rate at M'Baïki could increase or decrease as the class width increased, whether or not the mortality was an increasing function of size. Ramula and Lehtilä [9] also observed that l could increase or decrease when matrix dimensionality was reduced, with the former situation being more frequent than the latter [22].
At M'Baïki, the variations of l due to the choice of the class width were negligible with respect to the sampling variability of the estimate of l (i.e. the uncertainty on the estimate of l due to the finiteness of the data available). Therefore, at M'Baïki, it can be concluded that the choice of the class width (in the range 1-10 cm) had little importance to estimate the species-specific population growth rates l. The sampling uncertainty on the estimate of l may have been overlooked in matrix modeling (but see [49,[53][54][55][56]). For instance, Ramula and Lehtilä [9] gave the example of a change of l with matrix dimension by 0.006 for a Primula veris population and insisted on the difference in population size that this change induces in 50 years, while the standard error on the estimate of l was 0.026, which induces an uncertainty on the population size after 50 years that is much larger. When comparing the Dl values between tree and herbaceous species, these authors also considered only the between-species variability in Dl and disregarded all the within-species uncertainty on the estimates of Dl. Zuidema [22] also presumably considered only the between-species variability and disregarded the within-species uncertainty on l when comparing the values of l for different class widths (but the exact test was not specified). Enright et al. [8] did not consider the sampling uncertainty on l, as if transition rates were known exactly.

Are Matrix Models Reconcilable with other Models?
The way transition matrices are collapsed seems to influence the dependence of the population growth rate l on matrix dimensionality. The greatest dependence of l on matrix dimensions was observed when matrices were collapsed using (1) [8,9]. Equation (1) implicitly means that transition rates are estimated as proportions using data from the class of interest only, i.e. using the maximum likelihood estimator of the underlying Markov chain [57]. However, this proportion estimator of transition rates is known to have a large variance [58]. When compared to other modeling approaches like integral projection models, matrix models may then be expected to under-perform in terms of the precision of predictions [9,13]. Moreover, equation (1) to collapse transition matrices has the limitation that it depends on the distribution of individuals among classes. This may lead to undesirable results, e.g. the magnitude of changes in l with reduced matrix dimensionality is affected by the distance from the stable class distribution [9]. Theoretically, l depends on the transition matrix alone and not on the current distribution of individuals among classes [3].
Using (7) to collapse transition matrices implicitly means that transition rates are estimated on the basis of individual-based regressions for growth, mortality and recruitment over the entire size range [7]. In particular, with this approach, the number of parameters to estimate does not depend on the number of classes, which means that the variability of predictions does not depend on matrix dimension (contrary to [13]). As pointed out by Zuidema et al. [7], the use of regressions over the entire size range to estimate transition rates is also the basis of integral projection modeling, thus establishing a close connection between matrix models and integral projection models (IPM). In the same way as the Usher matrix model (5) is a discrete scheme to solve the continuous McKendrick equation (2), the numerical calculation of continuoussize IPM requires some discretization whose expression is a big transition matrix model [7,12]. In the same way as the discretization of the McKendrick equation into the Usher matrix model brings an error, the discretization of the IPM into a transition matrix model brings an error [7]. In fact, matrix models, IPM and the Fokker-Planck equation (depending on whether size and time are discrete or continuous) are equivalent in some limit [14,15], with the implication that all estimates that are dependent on class width (like age, see [22]) should tend to the same value when class width is small enough.
Another lesson to learn from this approach is that, although much attention has been devoted to the influence of the class width on predictions (i.e. the discretization of size), the time step (i.e. the discretization of time) may also influence the predictions. At M'Baïki, we collapsed bisannual transition data into annual transition rates so that the matrix model with the smallest class width be close to the McKendrick equation. Although the choice of the time step may bias predictions in the same way as the choice of the class width, and although changing the time step of a matrix model raises issues that are similar to collapsing its dimension [59], the influence of the time step in IPM and matrix models does not seem to have been studied.

Conclusion
We concur to conclude that matrix models should be used with narrow size classes, to be nearly equivalent with a continuous-size McKendrick equation [7]. The use of regressions over the entire size class and of equations (7) to estimate transition rates allows the modeler to decrease the class width d, with the only constraint on the lower bound of d that condition (8) must be met. At M'Baïki like in other studies [8,9,22], the choice of the matrix dimensionality had little influence on the population growth rate l. We showed that this influence was similar to that of a diffusion process, and did not act as a systematic bias towards fast pathways. Moreover, the change of l due to the class width was much less than the sampling uncertainty on the estimate of l. Therefore, the bias of l due to matrix dimensionality is not the only statistic to consider; the variance of the estimator of l should be considered as well. Making a parallel with the histogram that is more often used than continuous kernel estimators to estimate the density of distribution from a sample of data, searching for a trade-off between bias and sampling variance might lead to matrix models with size classes that are not so narrow.

Supporting Information
Table S1 Characteristics of population dynamics for 53 tree species at M'Baïki, Central African Republic. (PDF)