An Evolutionary and Ecological Community Model for Distribution of Phenotypes and Abundances among Competing Species

Here, we propose a theory for the structure of communities of competing species. We include ecologically realistic assumptions, such as density dependence and stochastic fluctuations in the environment, and analyze how evolution caused by r- and K-selection will affect the packing of species in the phenotypic space as well as the species abundance distribution. Species-specific traits have the same matrix G of additive genetic variances and covariances, and evolution of mean traits is affected by fluctuations in population size of all species. In general, the model produces a shape of the distributions of log abundances that is skewed to the left, which is typical of most natural communities. Mean phenotypes of the species in the community are distributed approximately uniformly on the surface of a multidimensional sphere. However, environmental stochasticity generates selection that deviates species slightly from this surface; nonetheless, phenotypic distribution will be different from a random packing of species. This model of community evolution provides a theoretical framework that predicts a relationship between the structure of the phenotypic space and the form of species abundance distributions that can be compared against time series of variation in community structure.


Introduction
Factors affecting the coexistence of competing species have received a great deal of attention by ecologists for nearly a century (Lotka 1924;Volterra 1926;Gause 1934). Still, how different processes such as interspecific competition affect local coexistence of ecologically similar species and hence influence the structure of communities is poorly understood (Chase 2010;Ellner et al. 2019). A major reason for this is that extrapolation of results from models of pairwise interactions between species to networks of interacting species has proved challenging (Levine et al. 2017). Furthermore, most models used in community ecology do not include realistic assumptions about factors affecting intraspecific variation, for example, with respect to abundance or distribution of phenotypes (Vellend 2010;Bolnick et al. 2011). Finally, many community models are purely ecological, which prevents analyses of the consequences of phenotypic evolution for the structure of the communities. Here, our aim is to propose a model for organization of communities that includes ecologically realistic assumptions about inter-and intraspecific dynamics and evolution based on general quantitative genetics theory for generating packing of species in the phenotypic space.
Theories for the structure of communities of competing species are based on the assumption that competing species reduce the population growth rate of one another through the exploitation of some critical common resource (McPeek 2017a(McPeek , 2017b. These consumer-resource interactions distribute the species along resource axes determined by the value of the species-specific carrying capacities and the competition coefficients (MacArthur 1970). The degree of limiting similarity between the competing species specifies the difference in the utilization curves that promote coexistence and is generally a function of variation in the resource supply (MacArthur and Levins 1967;Roughgarden 1974Roughgarden , 1976Loreau 2010, pp. 26-39). These analyses led MacArthur (1972, p. 237) to the general conclusion that, in competitive communities, species will reach equilibrium population sizes such that their resource utilization best matches the production of limiting resources. He further argued that not only should the number of individuals of the different species be maximized but they should also be uniformly spaced within the phenotypic space, allowing the presence of a maximum number of coexisting species. Thus, the general outcome would either be coexistence of the competing populations or spatial exclusion resulting in a checkerboard distribution of species in phenotypic space, dependent on their respective carrying capacities and the competition coefficients (Diamond 1975). However, empirical examination of this bottom-up theory of community structure has been difficult, mainly due to difficulties in measuring the required parameters, such as competition coefficients and carrying capacity in the absence of competitors under natural conditions, and accounting for alternative processes to competition, such as predation (Freckleton et al. 2009;Terborgh 2015;Aschehoug et al. 2016;Chesson 2018;Hart et al. 2018). MacArthur's (1970MacArthur's ( , 1972) models of packing of competing species were purely ecological and did not include evolutionary processes. They were based on Gause's (1934) presumption that no two species can occupy the same niche. These models were expanded by Roughgarden (1976) to include coevolution among competing species, based on a generalization of MacArthur's (1962) models for densitydependent selection acting within a single population to include multiple species. Another important theoretical development was provided by Abrams (1986), who included specific assumptions about evolutionary trade-offs between parameters describing the ecological dynamics. These models proved especially useful in analyses of character displacement, which was defined as mutual phenotypic divergence of species in response to interspecific competition (Brown and Wilson 1956). Slatkin (1980) then developed a quantitative genetic framework to analyze phenotypic divergence of two competing species. This and subsequent work have shown that character displacement can occur regardless of interspecific competition for limiting resources (Abrams and Cortez 2015). Recently, models of phenotypic evolution of competing species based on basic concepts from quantitative genetics have been expanded (Ellner et al. 2011;Govaert et al. 2016) and now also include consumerresource dynamics (see summaries in Vellend 2016, pp. 63-66;McPeek 2017bMcPeek , 2019. Still, the theory for phenotypic evolution in communities of competing species is based on simplified assumptions about important processes affecting the ecological dynamics. For example, fluctuations in population size may generate density-dependent selection ). This causes a trade-off between genotypes or phenotypes favoring rapid growth at low population densities, resulting in r-selection, and those that have higher fitness at population sizes close to the carrying capacity, causing K-selection (MacArthur and Wilson 1967;Pianka 1970Pianka , 1972Roughgarden 1971Roughgarden , 1979MacArthur 1972;Gilpin et al. 1976;Fel-senstein 1979;Boyce 1984;Mueller 1997Mueller , 2009). Such r-and K-selection strongly affects phenotypic evolution (Lande et al. 2009(Lande et al. , 2017Engen et al. 2013Engen et al. , 2020 and coevolution among competing species (Roughgarden 1976).
Another process affecting eco-evolutionary dynamics is random fluctuations in the environment (Lande 2007;Govaert et al. 2019). The effects of environmental variation on the outcome of competitive interactions between species was in particular highlighted by Chesson (1994Chesson ( , 2000aChesson ( , 2018Chesson ( , 2019, who proposed that temporal fluctuations in the environment can promote coexistence of competing species through the storage effect (i.e., that fitness benefits acquired during favorable periods are stored to resist unfavorable periods of population decline). Empirical analyses of the species composition as function of different functional traits along ecological gradients have also revealed that environmental filtering strongly affects the coexistence of potentially competing species by favoring those best adapted to the abiotic conditions (Freckleton and Watkinson 2001;Adler et al. 2013Adler et al. , 2018Kraft et al. 2015).
An important characteristic of communities, in addition to the packing of the species in the phenotypic space, is the distribution of abundances. Species abundance distributions, such as the classical Fisher's logarithmic series distribution (Fisher et al. 1943) and Preston's (1960) lognormal model, are now commonly described by a variety of different dynamic population models, such as Hubbell's (2001) neutral model driven by demographic stochasticity and models of single-species dynamics dependent on environmental fluctuations Lande 1996a, 1996b;Saether et al. 2013;Engen et al. 2017). Models including environmental stochasticity can be classified by their degree of heterogeneity, that is, the relative amount of variation in abundance among species determined by differences among species in dynamic parameters, often found to be large (Engen et al. 2002;Lande et al. 2003;Grøtan et al. 2012Grøtan et al. , 2014Solbu et al. 2018). Heterogeneity implies that differences among species are sufficiently large to generate significant differences in the parameters describing their population dynamics. Such differences may ecologically be a consequence of niche differentiation (Hutchinson 1959;May 1974), whereas at an evolutionary timescale it is the result of character displacement (Schluter 2000;Dayan and Simberloff 2005;Abrams and Cortez 2015).
Here, we develop a model based on a generalization of the concept of r-and K-selection (MacArthur and Wilson 1967;MacArthur 1972;Boyce 1984;Mueller 1997Mueller , 2009Reznick et al. 2002;Lande et al. 2009;Engen et al. 2013;Bassar et al. 2013;Saether et al. 2016) applied to a system of an arbitrary number of competing species with the degree of interspecific competition expressed by the generalization of the Lotka-Volterra model used by May (1974) and Roughgarden (1976). The dynamics of each species is determined by parameters given by their vital rates, which in turn are functions of the phenotype of each individual expressed as a vector of traits of an arbitrary number of dimensions. In addition to including intraspecific r-and K-selection, the model generates character displacement (Brown and Wilson 1956) by embedding interspecific competition coefficients that decrease as the phenotypes of species become more different, as in MacArthur's (1970MacArthur's ( , 1972 and May's (1974) classical models. The model includes competition coefficients that are not symmetric for pairs of species, which is required for a realistic description of the dynamics of interspecific competition (Roughgarden 1976(Roughgarden , 1979Brown and Vincent 1987;Case and Taper 2000;Goldberg and Lande 2006;McPeek 2017aMcPeek , 2017b. Finally, the model is stochastic and can be used to study the effects of fluctuations in the environment on the degree of character displacement as well as how variation in the population dynamics affects the species abundance distribution. Since the species end up with a distribution of traits among them, the abundance models are in general heterogeneous, with the degree of heterogeneity depending on consumerresource interactions (McPeek 2017a(McPeek , 2017b(McPeek , 2019 and environmental stochasticity as well as the evolutionary response to phenotypic selection.
Our approach contains the following steps. We start by presenting the model very generally through a direct generalization of the single-species model of Engen et al. (2013). Then, as is commonly done in quantitative genetic models (e.g., Charlesworth et al. 1982;Chevin et al. 2017), we express the fitness-related parameters by general seconddegree polynomials in the phenotype components, corresponding to MacArthur's (1970MacArthur's ( , 1972 classical competition models. For this model, we derive the mean Malthusian fitnesses (Saether and Engen 2015) for the species in the system. This is required to describe the species' dynamics as well as the phenotypic evolution using the general gradient formula for the multivariate response to selection (Lande 1976(Lande , 1979. However, because the model is complex and may require a huge number of parameters, we show some transparent illustrations based on simplifying approximations while still keeping some degree of ecological realism. These are based on transformations of the phenotypic space so that a number of matrices occurring in the description are proportional to the identity matrix. These assumptions lead to a transformed model with properties that can be investigated through simulations by varying only a small number of parameters, providing results that in the end can be transformed backward to express the findings in terms of the original phenotypes. This approach produces predictions about the distributions of abundances, as well as phenotypes, among the different species in the communities dependent on specific assumptions about processes affecting the joint dynamics of all species.

General Model
Consider a system of n species with dynamics given by their d-dimensional phenotypic vectors z i , i p 1, 2, ::: , n, where z i p (z i1 , z i2 , ::: , z id ) T ; in general a superscript T denotes matrix transposition (definitions of variables for the general model are provided in table 1). We assume that the number of phenotypic traits, d, is large enough and the traits appropriately chosen so that z i accurately determines individual fitnesses. The phenotypes have multinormal distributions f i (z i ) among the individuals of each species with mean values z i and the same covariance matrix P, both omitted to simplify the notation. The breeding values are also multinormally distributed with the same mean values as the phenotypes with identical additive genetic covariance matrices G across species. Using the Ito formulation of the stochastic differential equation (Karlin and Taylor 1981), the joint dynamics of a hypothetical community, where species i is of population size N i and with no variation in phenotypes (P p 0) within species, is solely determined by the set of phenotypic vectors, as is the Malthusian fitness of species i, r(z i ) is the deterministic growth rate in the absence of density regulation and competition, g ij (z i , z j ) for i ( j are the competition coefficients, and g ij (z i , z i ) is the effect of density when the species is alone (we denote z p (z 1 , z 2 , ::: , z n ) T and N p (N 1 , N 2 , ::: , N n ) T ). The quantities dB(t) and dB i (t) are independent increments of standard Brownian motions with zero mean and variance dt (Karlin and Taylor 1981); the first term is a noise term common for all species, while subscripts i indicate independent noise components specific for each species. The model can also be formulated in discrete time, as must be done for simulations, and then dt is replaced by small time interval Dt, while the corresponding stochastic components DB and DB i are normal variates with zero means and variances Dt. The populations are assumed to be large so that demographic stochasticity and random genetic drift can be ignored (Lande et al. 2003). The variance of the total noise acting on species i during time dt is j 2 e r dt 1 j 2 e (1 2 r) dt p j 2 e dt, so that j 2 e is the total environmental variance acting on each species. The correlation between the stochasticity for two species is E[j 2 e r dB(t) 2 ]=(j 2 e dt) p r. We assume that the noise parameters are not genetically determined and hence not subject to selection. Including such variation, as done by Lande (2007), Saether and Engen (2015), and Engen et al. (2020), would make the model slightly more complex and less transparent.
To describe the dynamics of z and N for covariance matrices P ( 0, we first define the mean parameters where bars in general denote mean values over all individuals of a species. Then it is shown in appendix A that the mean fitness for species i is so that the dynamics of population densities are given by because environmental stochasticity, by definition, acts equally on all individuals in a population (Engen et al. 1998).
Assuming constant phenotypic and additive genetic covariance matrices P and G, the response to selection of the phenotypic vectors for all species i p 1, 2, ::: , n given the population vector N and all mean phenotypes z is where ∇ i denotes the gradient with respect to z i . Because of the stochastic terms in the dynamic equations for the log population sizes (eq. [2]), the population vector N will fluctuate and generate stochastic components in the d z i , defining z as a stochastic process.

Competition Structure
We write the coefficients g ij in the general form The w are symmetric components of competition coefficients, w(z i , z j ) p w(z j , z i ), likely to decrease as the vectors z i and z j evolve to become more different, so that the two species diverge phenotypically and thus have less overlap in resource utilization (MacArthur 1970).
The assumption that w(z i , z j ) ! 1 for i ( j ensures that g ii (z i , z i ) 1 g ij (z i , z j ), so that the same equality holds for the mean coefficients. Hence, intraspecific competition for species i is always larger than its competition with species j ( i. We are not including description of resources in the present model but will still refer to w(z i , z j ) as a function depending on the diversity of resources. If there is a large diversity of available resources, then individuals as well as species can more easily show niche differentiation and avoid, or at least decrease, the strength of competition (Loreau 2010). Large diversity is then expressed mathematically by a w(z i , z j ) decreasing rapidly with increasing distance between z i and z j .
Phenotypic vector of species i (i p 1, 2, :::, n); z i,j is the phenotype distribution of species i for trait j z p (z 1 , z 2 , ::: , z n ) T d # n Phenotypic (two-dimensional) vector of community N i 1 Population size of species i N p (N 1 , N 2 , ::: , N n ) T n Population vector of community Environmental variance r(z i ) 1 Maximum growth rate of species i (i.e., in the absence of density regulation and competition) Symmetric components of competition coefficients, related to resource diversity g(z i ) 1 Effect of density and competition on individuals of species i r 1 Correlation between noise of two species Since the w are symmetric-and to get realistic modelsthe factor g(z i ) should be related to the intrinsic growth rate r(z i ), implying a trade-off so that individuals with large intrinsic growth rate are also strongly affected by density. At an evolutionary timescale, this will generate a balance between r-and K-selection (Lande et al. 2009;Engen et al. 2013Engen et al. , 2020. In other words, individuals that have genotypes favoring high growth rates will be most negatively affected by other individuals, either of the same or of competing of species. This trade-off may not be strict with exactly the same ordering of r(z i ) and g(z i ) among individuals, but realism requires some pattern of positive covariation of these parameters within species. If no such trade-offs were present, evolution would simultaneously increase r(z i ) and decrease g(z i ), implying unrealistic unlimited increase of individual fitness.
With these general assumptions, the competition between any pair of individuals or species is not symmetric; that The trade-off between r-and K-selection will make individuals and species with large values of r(z i ) and/or r( z i ) more negatively affected by other individuals and species than those that grow slowly under no influence of variation in population density. In general, individual z i affects the individual fitness of individual z j by a different strength than individual z j affects that of an individual with phenotype z i (asymmetric community matrix).

Specific Competition Model
Suppose, first, that all of the competition coefficients w(z i , z j ) are zero (for i ( j). The evolution of z i is then determined by the model analyzed by Engen et al. (2013). With appropriate trade-offs between the r(z i ) and g(z i ) preventing unlimited growth of fitness, evolution will then end up in stationary fluctuations of z i around a central value z * i maximizing [ r( z i ) 2 j 2 e =2]= g( z i ), which is actually the expectation in the stationary distribution of N i for a given constant value of z i (for G p 0).
We shall assume that the growth rates at small densities r(z i ) are second-degree polynomials in z i ; more precisely, where m p (m 1 , m 2 , ::: , m d ) T may be any vector of dimension d defining the phenotype z i p m with the maximum growth rate and A is some positive definite symmetric (d#d) matrix defining a p A 21 (we will use the same notation, with small and capital letters, for other matrices; see table 2). We may alternatively consider this as a second-order approximation using the Taylor expansion to the second order at the maximum m of r(z i ). The multiplicative deterministic growth rate at small densities, l p e r , is then proportional to the multivariate normal density with mean m and covariance matrix A (see app. C). This expresses how the growth rates r(z i ) decrease as phenotypes diverge from the phenotype m that gives the maximum growth rate r 0 . If there is no genetic variation in the density regulation (g is a constant ), then in the absence of stochasticity the species will be purely r selected, the growth rate will be maximized, and the mean phenotype stabilizes at z i p m.
More generally, we assume that density regulation acting through the factor g is genetically determined and ; maximum (i.e., at low densities) growth rate of species i as a quadratic function of its phenotype z i m d Phenotype of maximum growth rate a p A 21 d # d Symmetric positive definite matrix, sensitivity of maximum growth rate to phenotype g(z i ) 1 g 0 exp(2 1 2 (z i 2 n) T b(z i 2 n)); effect of community (i.e., intra-and interspecific) density on species i as a quadratic function of its phenotype z i n d Phenotype of maximum density effect b p B 21 d # d Symmetric positive definite matrix; b and a define a trade-off between maximum growth rate and effect of density w(z i , z j ) 1 exp(2 1 2 (z i 2 z j ) T c(z i 2 z j )); effect of competition between species i and species j as a quadratic function of phenotype difference z i 2 z j c p C 21 d # d Symmetric positive definite matrix, variability of resource strength, inversely proportional to strength of competition Additive genetic covariance matrix also has a Gaussian form (which can also be considered as a second-order approximation); that is, so that the effect of population density is at its maximum for z i p n. The matrices a and b and vectors m and n define some trade-off between r(z i ) and g(z i ), since phenotypes with very small growth rates r will also have very small g (z i deviating much from both m and n). This means that genotypes leading to slow growth under small densities (small r) are also those leading to a large intraspecific ability to compete (small g; e.g., Saether et al. 2016Saether et al. , 2021Kentie et al. 2020). However, note that this trade-off may in general not be strict in the sense that a change in phenotype corresponding to an increase in r does not necessarily lead to an increase in g, and the phenotype with the largest r in general is different from that with the largest g. However, the trade-off is strict if m p n and b p ba, where b 1 0; this is the model we later will use for illustrative purposes.
The factors w of the competition coefficients are expected to be smaller as the phenotypes z i and z j become more different. We also express this by a Gaussian function (or second-order approximation): where C p c 21 is a positive definite symmetric matrix. For a given z j , this competition coefficient, affecting the dynamics and evolution of z i , then has a Gaussian form in z i , with the maximum at z i p z j . As mentioned above, the form of this factor is dependent on the diversity of resources. Because we are working only with large populations, ignoring demographic stochasticity, the actual amount of resources acting through the parameter g 0 is not important in structuring the model. It only acts as a scaling of the abundances because the rescaled densities g 0 N j gives the same dynamics for abundances and mean phenotypes for any value of g 0 . However, if the diversity of resources decreases, it will become more difficult for individuals or species to find resources not heavily utilized by other individuals or species. Hence, if the elements of the matrix c increases, the species will experience less competition, which is consistent with more variability in the resource spectrum or higher diversity of resources.
The solution for this, requiring a number of complex integrations over functions proportional to the multivariate normal density, is given in appendix B. The result can be summarized by the gradients of the form For i p j, the same formula applies by plugging in c p 0. This determines the response to selection (from eq. [3]): N). The changes in population sizes are given by equation (2). Here, d i p z i 2 n and d ij p z i 2 z j , while J, K, and L are (d#d) matrices expressed by P, A, and B, as derived in appendix D (which also provides the derivation of the mean competition coefficients).
Under the assumption of weak selection, we show in appendix B that the responses to selection simplify to so that the complex matrices J, K, and L are not needed.

An Illustration of the Simplified Model
Although we can simulate this model for any set of parameters, we need to make simplifications to obtain some transparent information about the properties of the community. For a given phenotypic covariance matrix P, it is commonly observed that along axes where the phenotypic and additive genetic variances are large there will be small variation in fitness by a given change in phenotype; this phenomenon is often referred to as demographic buffering (Hilde et al. 2020). This observation is related to Fisher's (1930) fundamental theorem of natural selection, which ensures that the additive genetic variance in fitness at stasis is zero, while under environmental fluctuations around an equilibrium it will undergo small fluctuations likely to have rather similar amplitudes in any direction of the phenotypic space. If fitness changes much in a given direction, selection will also be strong in that direction, and individuals deviating much from the mean will have smaller fitnesses and thus be selected against (Schluter 1996) and generate small phenotypic variation in that direction. In directions with small variation in fitness, it is more likely that large phenotypic and additive genetic variances will be maintained. As we assume that the phenotypes within a species are multinormally distributed with covariance matrix P, the distribution is proportional to exp(2(1=2)(z i 2 z i ) T p(z i 2 z i )), determining the variance in any direction. Ignoring competition there are two fitness components, l(z i ) p exp(r(z i )) and g(z i ), assumed to have Gaussian forms with mean values m and n and covariance matrices A and B, respectively (for descriptions of variables, see table 2). In accordance with the above remarks, we obtain a class of models suitable for transparent illustrations of the theory by assuming that A and B equal P multiplied by some constants, A p a 21 P and B p b 21 P. Furthermore, assuming m p n then yields a model with a strict trade-off between l and g in the sense that phenotypes order their values equally. By this we mean that if one phenotype has a larger value of l than another, then it also has the larger value of g. More accurately, by the above assumptions, g is proportional to l b=a . We can also, without loss of generality, choose a new measurement scale by subtracting m from all phenotypes so that the maxima of l and g are located at z i p 0. Finally, assuming that the competition coefficients, which also have Gaussian forms, have the same structure, the model is given by a p ap, b p bp and c p cp, where a, b, and c are positive scalars, as well as m p n p 0. We emphasize that we have no strong empirical support for these assumptions, but we still believe that they represent a nice simplified model with rather few parameters that is useful for transparent illustrations likely to be relevant for the more general model. To analyze this model, it is appropriate first to introduce transformed variables by a Cholesky decomposition z i p Du i , as explained in appendix F, so that the phenotypic covariance matrix for the new variables u i is the identity matrix I. To avoid technical complications causing distraction from the main message, we also make the simplifying assumption that G p h 2 P, meaning that the variance of the environmental component of the standardized phenotype in any direction equals a proportion 1 2 h 2 of the phenotypic variance. In other words, the heritability is independent of direction in the transformed phenotype space. Then, the additive genetic covariance matrix for the transformed phenotypes u i is simply h 2 I. It is shown in appendix F that the dynamics of ln N i and u i then takes the simplified form where g( u i , u j ) is given by equation (F1), while the response to selection is where for i ( j while for i p j the same formula applies by plugging in c p 0. Here, J, K, and L are now scalars, defined in appendix F (they are multiplicative factors of the identity matrix, allowing us to get the associated matrices J, K, and L). Although this model is based on a number of simplifying assumptions, it still captures the basic structure of the original model but maintains a transparent structure appropriate for further investigations. This occurs because the only parameters included are r 0 , a, b, and c, in addition to j 2 e and r defining the stochasticity in the dynamics of N i and the predetermined number of phenotypic dimensions and initial number of species, d and n, respectively.

Scaling of Densities
All density-regulating terms are proportional to the parameter g 0 . Writing W for the total amount of available resources, we can define the population density of individuals with respect to the amount of resources available as the number of individuals per resource unit, N j =W. This is a general approach in models of resource-consumer interactions (McPeek 2017a(McPeek , 2017bEngen et al. 2020). Hence, by choosing g 0 p 1=W, the model for the dynamics of the densities is equivalent to that for N with g 0 p 1. This shows that g 0 is a pure scaling parameter. We are then free to consider N as densities or to choose N j to be the number of individuals for some specific value of g 0 . In the latter case, g 0 can be changed to regulate the size of the community, which on average is inversely proportional to g 0 . Hence, this parameter is crucial for the extinction rate but otherwise does not affect the properties of the competition system.

Distribution of Mean Phenotypes in the Absence of Competition
Next consider the parameters in the diagonalized model. The parameters r 0 , a, b, and j 2 e define the single-species dynamics in the case of no competition, leaving the competition parameters c and the noise parameter r as specific for the competition system. For very large values of c, it appears that g( u i , u j ) for j ( i is small and can be neglected, meaning that dynamics and evolution for all species are practically the same as for single species considered alone. The biological interpretation of this is that there is quite large resource diversity so that even species with quite similar phenotypes can find practically nonoverlapping resource niches, corresponding to the results from classical competition without any evolutionary dynamics (MacArthur and Levins 1967;MacArthur 1970MacArthur , 1972. If there is no competition, the dynamics of each species is described by the logistic model, and evolution tends to make the phenotypic vector of each species fluctuate around a value u * i maximizing the expected population E(N i ) size , determined by E(d ln N i ) p 0. This gives, from equation (F1), It appears that any vector u * i satisfying is a solution by maximizing E(N i ), showing that evolution will end up fluctuating around some phenotype close to the d-dimensional sphere with radius R, corresponding to the initial mean phenotypes z i being located at some arbitrary point close to the surface of an ellipsoid. In other words, u * is not uniquely determined in this model but may take any value on the sphere surface. The stochastic fluctuations of R will increase with increasing environmental stochasticity and decreasing strength of selection, occurring, for example, when r 0 gets smaller Engen and Saether 2019). When there is practically no competition, which occurs when there is large resource diversity (large values of c), the species will follow approximately the same logistic dynamics in ecological time, which is a model known to give species abundance distributions of the gamma type (Engen and Lande 1996b). For small shape parameters in the gamma, the corresponding distribution of log abundances is known to be skewed to the left, as observed in many empirical studies (Williams 1964;Ulrich et al. 2010;Saether et al. 2013).

Competition in Relation to Resource Diversity
This preliminary investigation gives some intuitive understanding of the more general model with smaller resource diversity c. The model construction is performed in such a way that any competitor with mean phenotype z j generates a component of d z i that increases the distance between the two. Generally, the mean of one species will be affected by the others in a way that on average makes its mean phenotype diverge from the others, although the effect of each specific competitor differs as a result of variation in mean phenotypes as well as population sizes among the species. As a consequence, for small values of c, the presence of competitors will increase the density-regulating term in the dynamic equations, which on average will make abundances smaller. Species that become most rare are then likely either to stay quite rare or to become extinct if some extinction barrier is introduced. When a sufficient number of species have become rare or have disappeared, the remaining phenotypic vectors z i for the more abundant species in the community are likely to end up being approximately uniformly distributed on the sphere because of the abovedescribed divergence of phenotypes. However, even in the deterministic case a fully uniform distribution of points on a sphere does not actually exist except for some specific number of species. Moreover, our stochastic model is complex, so the radius R is also likely to be somewhat affected by competition, and the stochastic fluctuation in N generates stochasticity in the evolution so that the points are likely to undergo stochastic fluctuations away from a fixed sphere, generating variation in mean population size among species.

Stochastic Simulations
Because of the complexity in the model, the dynamics of the community composition can be investigated only by extensive stochastic simulations. Figure 1a shows an example from such simulations of a distribution of two components of u (two out of five transformed mean phenotypic components for all species) after simulating the system for a long period of time. From this graph we cannot see how the radius R is distributed among species because of the stochasticity in the model, since it is only a projection onto a lower-dimensional (!d) space (namely, a plane here), but the actual distribution of R is shown in figure 1b. Since the variance in R is rather small, the phenotypic projections take values approximately located within a circular disk, as expected, meaning that the original mean phenotypes are roughly distributed within some ellipse. Accordingly, species may look approximately evenly separated if we observe only two out of many dimensions, while in reality they are expected to be approximately evenly spaced on the surface of a multidimensional sphere. The log population sizes for four of the species for a period of 1,000 years, after first simulating the community for 10,000 years, are shown in figure 2, while the evolution of one phenotypic component for eight randomly chosen species for the following 10,000 years are shown in figure 3.
If distributions of mean phenotypes among species transformed to give the diagonalized model are estimated from data, it may be of interest to compare it with a multivariate normal distribution with independent components with the same mean and variance. If the components have variance j 2 and are translated to have zero means, the distribution of R 2 =j 2 is known to be x 2 distributed with d degrees of freedom. Accordingly, R=j has the so called x distribution (Johnson et al. 1994). The mean and variance of R is then E(R) p j ffiffi ffi 2 p Γ[(d 1 1)=2]=Γ(d=2) and var(R) pj 2 [d 2 (E(R=j)) 2 ]. In the simulations shown in figure 1, it appears that the mean of R is approximately 2.81, which, with d p 5, means that j ≈ 1:79 and thus SD(R) ≈ 1:23, which is considerably larger than what ap-pears from figure 1, which is about 0.02. This means that with the parameters used in this example, competition forces species to be distributed much more closely to a sphere surface than under a multivariate normal distribution of mean phenotypes among species. Weaker selection (smaller a and/ or b) and larger environmental variance will increase the variance of R. Two species may coexist even if they have very similar values of one phenotype component (fig. 3). The reason for coexistence is then that they have different values of other components affecting fitness so that their phenotype vectors have different directions in space.

Species Abundance Distributions and Heterogeneity of Communities
These coevolutionary processes generate specific distributions of the log abundances of the species in the community ( fig. 4a). To find the average abundance of each species we could in principle calculate the expectation of the right side of the dynamic equations and solve for the expected values of the N j . However, we would then have to invert a very large matrix likely not to give stable solution, perhaps also even negative values for some mean abundances. However, we are not primarily interested in mean abundances but rather mean log abundances. Hence, we have used the more direct approach of simulating the system over an additional 1,000 years, including neither evolution (h 2 p 0) nor environmental stochasticity (j 2 e p 0), and then recording the mean log abundance of each species. This distribution determines the genetic differences of species with respect to the mean log abundance they code for, reflecting the heterogeneity of the community. The variance of the distribution of log abundances ( fig. 4a) can then be decomposed into a component, var(ln E(N)), due to genetic differences among species (fig. 4b) and one due to the environmental fluctuations in log abundances for given mean phenotypes. These components add up to the total variance var(ln N) (see fig. 2). The ratio between var(ln E(N)) ( fig. 4b) and var(ln N) ( fig. 4a), which in this example is 0.56, is the fraction of the variance in the species log abundance distribution explained by genetic differences among species and is called the heterogeneity of the community.
The major difference between this model and the evolutionary model of r-and K-selection by including only a single species is the modeling of competition through the competition coefficients w(z i , z j ). These terms generate character displacement (Brown and Wilson 1956;Dayan and Simberloff 2005), reducing the effects of interspecific competition (e.g., Roughgarden 1974Roughgarden , 1976 as well as lowering the risk of extinction of species that, through stochastic effects, evolve phenotypes incapable of keeping up large abundances. This enables us to examine how the  Figure 1: a, Distribution among species of transformed mean phenotype components u i1 and u i2 in the diagonalized model for a community of 119 species. Simulations were first run for 10,000 years without extinction barriers, starting with 200 species, and then were run for another 1,000 years with extinctions at density 0.5. The first period is a burn-in simulation to stabilize phenotypes and population sizes, while the next step with an extinction barrier is arbitrary but has been chosen large enough to make the community approximately stable with a very small extinction rate. During this time, 81 species reached a density below 0.5 and were considered extinct, and they did not later affect the remaining species. The extinction barrier is arbitrarily chosen because the N i are densities scaled by the parameter g 0 , so that the actual number of individuals is not really in the model. Parameter values are d p 5, h 2 p 0:5, r 0 p 0:2, g 0 p0:0001, a p 0:02, b p 0:5, c p 0:15, r p 0:5, and j 2 e p 0:01. b, Corresponding distribution of R i p À P d rp1 u 2 ir Á 1=2 . variance of log species abundances as well as the number of species are affected by resource diversity expressed by the parameter c. In the diagonalized model, the variance in log abundances decreases with increasing c, hence reducing the strength of competition ( fig. 5a). In other words, increasing resource diversity generates larger species diversity. This occurs not only via a reduction of variance in log abundances but also through an increase in the number of species ( fig. 5b) that can be kept up by these resources. Stochastic variation in the environment strongly affects the dynamical consequences of interspecific competition (Chesson 2000a(Chesson , 2000b. Accordingly, the variance in log abundances tends to increase with increasing environmental variance j 2 e ( fig. 6a) and reduces the number of species ( fig. 6b).

Discussion
Here we have developed a theoretical framework that integrates coevolutionary processes of interacting species using an evolutionary model based on ecological realistic assumptions. This general stochastic model includes the dynamics of n competing species (eq. [1]) as well as the response to selection with respect to a d-dimensional phenotypic vector for each species (eq. [2]). Accordingly, the model embeds basic processes affecting intraspecific var-iation in abundance as well as interactions among species. Two general conclusions appear from our analyses. First, the species tend to distribute themselves surprisingly stably within the phenotypic space, often approximately uniformly on the surface of a d-dimensional sphere-an alternative to the pattern of species packing suggested by MacArthur (1972). Second, species interactions as well as the potential for niche differentiation also affect another important characteristic of community structure, the species abundance distribution. A general pattern that appears is that increased resource diversity increases the number of species in the community and also reduces the variance in log abundance (fig. 5). A similar effect occurs if environmental stochasticity is reduced ( fig. 6).
Although the model given by equations (1) and (2) is very general, our illustrations are based on several simplifying assumptions to increase transparency. Phenotypic evolution is determined by matrices-matrices occurring in the Malthusian fitness computation as well as the additive genetic covariance matrix-that are proportional to the phenotypic covariance matrix. Furthermore, we describe the evolutionary dynamics by the simple Breeder's equation (Lande 1976) G p h 2 P, so that the heritability in any direction in the transformed phenotypic space is the same. In ecology, Hubbell (2001) proposed a different approach, assuming that metacommunities of identical species had dynamics affected only by demographic stochasticity. This model produced several simple relationships, which were difficult to refute by empirical analyses and provided important general insights into processes affecting variation in species diversity by highlighting the importance of ecological drift (Rosindell et al. 2011(Rosindell et al. , 2012. Recently, processes affecting intraspecific variation have received increased focus in analyses of community structure (Bolnick et al. 2011;Levine et al. 2017). In our approach, we include a specific model for density dependence acting within the species, describing how individuals respond to variation in the number of conspecifics. The densityregulating terms could have been modeled in many different ways (Holt 2020), for example, by choosing competition coefficients g rather than ln g to be a second-degree polynomial. This would actually be simpler to handle mathematically, but this form was chosen to give the model realistic structure for any legal set of parameter values. This ensures that the population growth rate l is always positive (since l b=a is proportional to g). Furthermore, we assume that phenotypes with the highest growth rates are also those most strongly affected by variation in population size (i.e., m p n, which leads to a strict trade-off between the deterministic multiplicative growth rate l and the common factor g in the density regulating terms; see Travis et al. 2013). Hence, comparing two individuals, the one with vital rates at small densities giving the largest l is also the one for which the growth rate is most strongly reduced by increasing population density (for empirical examples, see Saether et al. 2016Saether et al. , 2021Kentie et al. 2020). These assumptions enable linear transformation of the phenotypic vectors giving a diagonalized model with only seven parameters regardless of the value of phenotypic dimensions, a, b, c, r 0 , h 2 , j 2 e , and r, in addition to a pure scaling parameter, g 0 , with little effect on the dynamics of the system. Thus, these assumptions increase the transparency of the model dramatically and have made it possible for us to reveal, by extensive stochastic simulations, some general relationships between phenotypic evolution and community structure. For example, the underlying eco-evolutionary dynamics of species strongly affect the species abundances distribution (fig. 4). The shape of the distribution of log abundances ( fig. 4) is slightly skewed to the left, which is typical for many natural communities (e.g., Williams 1964;Saether et al. 2013). It also includes a component due to permanent heterogeneity among species of similar magnitude, as estimated in natural communities (Engen et al. 2002;Lande et al. 2003;Grøetan et al. 2012Grøetan et al. , 2014Solbu et al.  The first simulations were performed for 10,000 years, leaving 113 species remaining above the extinction barrier, and the model was then simulated for another period of 10,000 years. The dashed lines show two of these species that became extinct during this period. Note that figure 1a, which is a different simulation of the same model with the same parameters, shows a snapshot of the distribution among species of two phenotypic components after 10,000 simulations, each varying over the same range as the components shown in this figure. 2018). Furthermore, one general prediction from our model is that species should be approximately evenly spaced out, close to the surface of a d-dimensional sphere in the phenotypic space ( fig. 1a). Environmental stochasticity tends, however, to move species away from the surface, but still the variance in the distribution of the phenotypes will be far less than expected from a random species packing ( fig. 1b).
Another general outcome of the model is that character displacement (Brown and Williams 1956;Dayan and Simberloff 2005;Abrams and Cortez 2015) generates a surprisingly stable distribution of phenotypic traits ( fig. 3), with transformed phenotype components distributed among species approximately on a d-dimensional sphere ( fig. 1), which after a backward transformation in general forms an ellipsoid. Since the radius R shows small interspecific variation ( fig. 1b), the five-dimensional phenotypes are approximately located uniformly on the surface of a sphere. This is the expected strategy for optimal reduction of com-petition by phenotypic divergence so that the species on average deviate as much as possible from each other while still having growth rates large enough to avoid extinction (MacArthur 1972;Loreau 2010). In fact, if d increases towards infinity and R is constant (var(R) p 0), then a uniform distribution of points on the sphere approaches exactly the multivariate normal distribution of Euclidian components. Thus, our model provides an extension of previous ecological models of species packing (e.g., MacArthur and Levins 1967;MacArthur 1970;Loreau 2010, pp. 26-32) that is based on modifications of the Lotka-Volterra equations but includes no evolution. Similarly, most evolutionary models for interspecific coevolution of traits based on consumer-resource or predator-prey models are deterministic (Vellend 2016;McPeek 2017aMcPeek , 2017b, ignoring stochastic fluctuations in the environment. Still, we assume that c describes competition continuously in phenotype differences. To illustrate the power of our approach, we have here focused on two parameters describing ecologically important processes, the parameter c occurring in the competition coefficients and the environmental variance j 2 e expressing the effects of fluctuations in the environment (Chase 2010). We have called c the diversity of resources since it is likely to be associated with the strength of interspecific competition, with no competition (competition coefficients close to 0) for very large values of c and very strong competition (coefficients close to 1) for small values of c. With increasing variation in resource types and hence in the number of available niches, competition decreases, which can be expressed by an increase in the value of c. However, in many cases the niche space will, to some extent, consist of discrete entities so that clusters of species (guilds) may evolve within the same niche type and coexist even if they are quite similar compared with other species in the community (Scheffer and van Nes 2006). This will cause clumping of species within the niche space. In addition, variation in c will also affect the species abundance distributions. In general, species diversity is commonly considered to have two major components (Magurran 2004), the number of species and the equitability (the evenness in abundances among species), both contributing to large species diversity. Variance in log abundance decreases (equitability increases; fig. 5a) and the number of species remaining in the community after 1,000 years increases ( fig. 5b) with increasing c. Hutchinson (1978, pp. 155-161), Chase and Leibold (2003, pp. 14-16), and Peterson et al. (2011, pp. 10-15) have made the important distinction between components of the niche space that are dynamically unaffected by competition and the components modified by species-specific interactions. Here, we show that both species richness and equitability increase with increasing c. On the other hand, the amount of resources will mainly The variance in log abundance is 3.04, considerably larger than the variance of the mean log abundances, which is 1.70. Mean log abundances were found by simulating the community for another 1,000 years without selection and environmental stochasticity (h 2 p 0, j 2 e p 0) and recording the mean log abundances for all species. The degree of heterogeneity of the community, which is the fraction of the total variance explained by genetic differences among species, is the ratio of 0.56 between these variances. For parameter values, see figure 1. affect the parameter g 0 , which is a pure scaling parameter regulating population and community size. Consequently, this will also have an effect on species diversity by regulating the extinction processes in the sense that smaller populations are more likely to become extinct (Lande et al. 2003). According to the Red Queen hypothesis (Van Valen 1973;Dawkins and Krebs 1979;Stenseth 1979;Stenseth and Maynard Smith 1984), species must continuously adapt to an ever-changing environment because of continued evolutionary changes occurring in all of the competing species. Extinctions and speciations may thus produce an everlasting "arms race" between species with no final winners.
Fluctuations in the environment generate temporal variation in population size, which in turn affect the evo-lutionary process given by equation (2). Increasing the environmental variance reduces the species diversity by generating an increase in the variance of log abundances ( fig. 6a) and a decrease in species number ( fig. 6b). This corresponds to the environmental filtering effect, which arises from correlations between environmental gradients and distribution of traits within communities (Kraft et al. 2008;Cornwell and Ackerly 2009;Adler et al. 2013). The heterogeneity among species, however, which is expressed by the ratio of the two variances shown in figure 6a, is rather stable over a large range of environmental fluctuations.
In our model, the dynamics of the community are determined by the set of phenotypic vectors. This enables us to examine how specific assumptions about the effects of trait vectors on ecological dynamics, as illustrated by our numerical example, will affect the distribution of species in the phenotypic space ( fig. 1). Furthermore, these underlying assumptions will also influence the form of the species abundance distributions (figs. 4, 5a). We therefore expect that comparative analyses of temporal variation in community dynamics (see examples in Saether et al. 2013) will provide relationships between the pattern of distribution of species in the phenotypic space and the shape of the species abundance distribution. Relating such distributions to general ecological characteristics, such as environmental stochasticity and the potential for niche differentiation, should make it possible to make inferences about processes most strongly affecting variation in structures of communities of competing species. This means that our model can be used to include historical effects (Chase 2003), for example, making the ghost of competition past (Connell 1980) more visible than can be revealed by using null models for the simultaneous occurrence of competing species (Chase and Leibold 2003).
To summarize, although the model presented here is purely theoretical, the ideas and general findings can be important in identifying processes affecting variation in the structure of natural communities. One important task that emerges from this study is to analyze relationships between structures in phenotype distributions within communities and species abundance distribution as well as relating those patterns to general characteristics affecting variations in availability and diversity of critical resources. For example, we expect that temporal variation in the environment will affect the pattern of species packing in the phenotypic space as well as the form of species abundance distribution. A challenge will be that empirical examination of the patterns derived from these theoretical explorations will require a long-term time series of variation in community composition that includes estimates of both abundances and structures of the phenotypic space.

Acknowledgments
We are grateful to R. Salguero-Gómez for valuable input. The study was funded by a Centre of Excellence grant from the Research Council of Norway (SFF-III, project 223257).

Statement of Authorship
S.E. and B.E.-S. developed the idea. S.E. with support from C.F.D.C. developed the model, and V.G. wrote the R programs. All authors contributed to the interpretation of the results from the analyses and to the writing of the article.

APPENDIX A The General Expression for Mean Malthusian Fitnesses
Now, for a continuum of phenotypic values within species (P ( 0), consider the subpopulation N i f i (z i )dz i of species i with phenotypes in (z i , z i 1 dz i ) (N i is the total number of individuals of species i). The reduction in growth rate of this population due to competition with individuals of species j with phenotypes in (z j , z j 1 dz j ) is assumed to be g ij (z i , z j )N j f j (z j ) dz j . Summing the effects of all individuals of species j then yields the reduction equal to the integral of this quantity over all z j . The reduction in growth rate is The total mean logarithmic growth rate of the subpopulation The corresponding mean logarithmic growth rate of the total population for species i is then the mean of these growth rates with respect to the distribution of z i : which is equivalent to equation (1).

APPENDIX B Solution for the Specific Competition Model
Our goal is to find the results required to describe and perform joint simulation of the dynamics and evolution as expressed by equations (2) and (3); that is, we need to compute the species-specific mean Malthusian fitnesses m i ( z , N) and their gradients ∇ i m i ( z , N) in order to apply the general gradient formula for response to selection derived by Lande (1979). Appendix C gives some required basic general results for the multivariate normal distribution. For our specific second-degree model, the mean of r(z i ) is The expression for the mean value of g(z i , z j ) is most transparently expressed by the mean value deviations d i p z i 2 v and d ij p z i 2 z j and the matrices q p p 1 b 1 c and v p p 1 c 2 cQc (see definition of Q and q in app. D). For i ( j, this yields the competition coefficients where and J, K, and L are (d#d) matrices expressed by functions of P, A, B, and C, derived in appendix D (eqq. [D9]-[D11]). The within-species coefficients for density regulation g( z i , z i ) are given by the same formula by plugging in c p 0.
The gradients have the form and for i ( j For i p j, the same formula applies with c p 0. This determines the response to selection, d z i p G∇ i m( z , N), while the changes in population sizes are given by equation (1). Under the assumption of weak selection-that is, small elements of P compared with those of B and C-giving Q ≈ P and V ≈ P as well as the exponent in equation (B3) being small compared with the exponent in equation (B2), it appears that Γ(d i , d ij ) ≈ 1, so that g( z i , z j ) ≈ g( z i , z j ), in accordance with standard first-order theory. The responses to selection then simplifies to Under stronger selection, Γ(d i , d ij ) expresses a more complex effect on the mean Malthusian fitness of phenotypic variation within species, making the coefficients of d i and d ij smaller than those appearing in the formula for weak selection.

APPENDIX C Some General Properties of the Multivariate Normal Distribution
If a multidimensional variable z i (in this appendix, we use the notation for the phenotype of species i, as per our model) has multivariate normal distribution N ( z i , P), with mean vector z i and positive definite covariance matrix P with inverse p p P 21 , then its distribution has density: The moment-generating function E[exp(t T z i )] p exp(t T z i 1 (1=2)t T Pt) or E[exp(t T x i )] p exp((1=2)t T Pt), where x i p z i 2 z i and t p (t 1 , t 2 , ::: , t d ) T . Equivalently, this yields the integral Ð exp(t T x 2 (1=2)x T px) dx p (2p) d=2 jPj 1=2 exp((1=2)t T Pt).
From this, it appears that when x is multivariate normal with zero mean and covariance matrix P, then for any positive definite (d#d)-matrix M where m p M 21 , r p p 1 m, and R p r 21 .

APPENDIX D
Derivation of Matrices J, K, and L Taking into account the variation in phenotypes among individuals, the dynamics are in general given by the mean Malthusian fitnesses m i ( z , N) p where f is the multivariate normal density and the integral runs over all (d#n) dimensions. For the model given by equation (1), this yields m i ( z , N) p r( z i ) 2 j 2 e 2 2 X n jp1 g( z i , z j )N j : First, the mean value of r is r( z i ) p r 0 2 E 1 2 ( z i 2 m) T a( z i 2 m) ! , p r 0 2 1 2 X uv P uv a uv 2 1 2 ( z i 2 m) T a( z i 2 m): To find the mean Malthusian fitness, we also need to find the mean values of the competition coefficients for i ( j, We compute this using double expectation E[g(z i , z j )] p E j E i [g(z i , z j )jz j ], where subscripts i and j denote expectations with respect to z i and z j , respectively. Rewriting the terms in the exponent of equation (D4) yields where d i p z i 2 v and 1 2 (z i 2 z j ) T c(z i 2 z j ) p 1 2 x T i cx i 1 ( z i 2 z j ) T cx i 1 1 2 ( z i 2 z j ) T c( z i 2 z j ): Using these expressions, the exponent of equation (D4) takes the form where u j p 2bd i 2 c( z i 2 z j ) and a j p d T i bd i 1 ( z i 2 z j ) T c( z i 2 z j ). Now, conditioning on z j so that u j and a j are constants not depending on x i , using equation (C2) and writing q p p 1 b 1 c then gives E i [g(z i , z j )jz j ] p g 0 jQj 1=2 jPj 1=2 exp 2 1 2 (a j 2 u T j Qu j ) : Finally, this must be integrated over the distribution of x j . The quantities a j and u j expressed by x j p z j 2 z j are a j p x T j cx j 2 2d T ij cx j 1 d T ij cd ij 1 d T i bd i , where d ij p z i 2 z j , and m j p cx j 1 w, where w p 2bd i 2 cd ij , giving m T j Qu j p x T j cQcx j 1 2w T Qcx j 1 w T Qw: The exponent in equation (B5) is then 2 1 2 (a j 2 m T j Qu j ) p 2 where s p c2cQc, y p c(d ij 1Qw) p (c2cQc)d ij 2cQbd i , and b p d T i bd i 1 d T ij cd ij 2 w T Qw p d T i (b 2 bQb)d i 1 d T ij (c 2 cQc)d ij 2 2d T i bQcd ij . Again, using equation (C2) in combination with equation (D5) and the fact that x j is multivariate normal with zero mean and covariance matrix P then gives where v p s 1 p. Evaluating the exponent in equation (D6) then yields g( z i , z j ) p g 0 exp 2 where J p bQb 1 bQcVcQb, ðD9Þ K p cQc 1 (c 2 cQc)V(c 2 cQc), ðD10Þ L p bQc 2 bQcV(c 2 cQc): Matrices J, K, and L are of dimensions (k#k). The within-species density regulation g( z i , z i ) is computed by using the same formula and plugging in c p 0.

Computation of the Response to Selection
According to the general theory presented in the main text the response to selection of species i is From equation (D3) we get ∇ i r( z i ) p 2a( z i 2 m): From equation (D6) we see that and from equations (D6) and (D7) we find for i ( j that For i p j the same formula applies by plugging in c p 0.

APPENDIX F The Simplified Illustrative Model
To analyze the simplified model defined in the main text, it is appropriate first to introduce transformed variables by a linear transformation z i p Du i , where D is an upper triangular Cholesky matrix obeying D T D p P so that the covariance matrix for u i is the (d#d) identity matrix I. This construction is as follows. For the symmetric definite positive matrix P, we want to generate an upper triangular Cholesky matrix D, such that D T D p P. From D, it is then possible to get u i p D 21 z i , so that the covariance matrix for u i is identity matrix I.
To construct matrix D, we will actually look for upper unit triangular matrix E (diagonal is made of 1s) and diagonal matrix M, such that P p E T ME. The following recursive iteration allows the construction of such matrices J and E: From matrices M and E, one gets D p M 1=2 E With this transformation, the components of u i for any species i is accordingly a vector with independent normally distributed components among the individuals of the species, with unit variance and mean u i p D 21 z i p P 21 D T z i . With the additional simplifying assumptions that G p h 2 P, the additive genetic covariance matrix for the transformed phenotypes u i is simply h 2 I. Now, expressing our model in terms of the u i , it appears to maintain the same mathematical structure as before, but all matrices occurring in the solution given by equations (4) and (5) are multiples of I. For this completely diagonalized model we find, by plugging into the general solution given in appendixes B and C, that q p qI, v p vI,