A matrix model for density-dependent selection in stage-classified populations, with application to pesticide resistance in Tribolium

The study of eco-evolutionary dynamics is based on the idea that ecological and evolutionary processes may operate on the same, or very similar, time scales, and that interactions of ecological and evolutionary processes may have important consequences. Here we develop a model that combines Mendelian population genetics with nonlinear demography to create a truly eco-evolutionary model. We use the vec-permutation matrix approach, classifying individuals by stage and genotype. The demographic component is female dominant and density- dependent. The genetic component includes random mating by stage and genotype, and arbitrary effects of genotype on the demographic phenotype. Mutation is neglected. The result is a nonlinear matrix population model that projects the stage×genotype distribution. We show that the results can include bifurcations of population dynamics driven by the response to selection. We present analytical criteria that determine whether one allele excludes the other or if they persist in a protected polymorphism. The results are based on local stability analysis of the homozygous boundary equilibria. As an example, we use a density-dependent stage-classified model of the flour beetle Tribolium castaneum . Our model permits arbitrary life-cycle complexity and nonlinearity. Tribolium has developed resistance to the pes- ticide malathion due to a dominant allele at a single autosomal locus. Using parameters reported from laboratory experiments, we show that the model successfully describes the dynamics of both resistant and susceptible homozygotes, and the outcome of a selection experiment containing both alleles. Stability analysis of the boundary equilibria confirms that the resistant allele excludes the susceptible allele, even in the absence of malathion, agreeing with previously reported results.


Introduction
The demographic processes of birth and death drive changes in gene frequencies and changes in population density and structure. Demography is therefore central to understanding ecology and evolution, and eco-evolutionary analyses always strive to incorporate the fundamental demographic processes of birth, death, and development (e.g., Coulson et al., 2006;Metcalf and Pavard, 2007). de Vries and Caswell (2019b) recently introduced an eco-evolutionary framework that combines matrix population models with basic Mendelian genetics. Here we use this framework to explore density-dependent selection.
Density dependence occurs when the per-capita vital rates (rates of birth, mortality, and development) depend on population size or density. Density dependence may be negative or positive (Allee effects). In models that contain demographic structure, "density" is a multivariate concept. Vital rates may depend on the abundance of a particular stage or age class, or on a weighted density that gives distinct weights to different stages rather than the total population density (Caswell et al., 2004). For example, density-dependent effects due to difficulty in finding mates leads to fertilities that depend on the abundances or densities of reproducing stages. Cannabilism is usually restricted to large individuals eating smaller conspecifics, so density-dependent effects due to cannibalism are size-or stage-specific. We will assume that genotypes can differ in age-or stage-specific rates of development, survival, or fertility anywhere in the life cycle. Thus, we assume that pleiotropic effects are the rule rather than the exception.
(1971) extended the logistic equation to multiple competing genotypes, and showed that selection leads to an increase in population density in a constant environment because only alleles with heterozygote advantage in the carrying capacity can invade. This result leads to the ideas of r-and K-selection. Charlesworth (1971) extended this to any fitness function that decreases with population densities. Charlesworth (1994) used high-order difference equations to model density-dependent selection in age-structured populations. In this paper, we will incorporate age or stage structure by using matrix population models rather than scalar difference equations.
In this paper we show how to construct a density-dependent Mendelian matrix population model, based on genotype-specific demographic measurements. We show how to use that model to project the stage × genotype distribution, and derive analytical conditions that determine whether alleles will coexist in a genetic polymorphism, or if one or another allele will go to fixation. We apply the analysis to a study of pesticide resistance in Tribolium beetles, a species that is a pest of stored grain products. We investigate the effect of incomplete dominance on the speed of invasion and on the outcome of invasion using evolutionary stability analysis.

Model construction
Individuals are jointly classified by stage (1, …, ω), and by genotype (1, …, g). Each genotype is characterized by a survival and transition matrix and a matrix of fertility rates. Both the transition matrices and the fertility matrices are assumed to be density dependent. In this section, we will consider general density dependence without specifying a functional dependence. In Section 3, we apply the general results obtained to a specific model with density-dependent demographic rates.
We make the important assumption of female demographic dominance; i.e. we assume that enough males are always present to fertilize all the females and that the number of offspring produced in a mating is not affected by the stage or genotype (i.e. the i-state) of the male. This assumption can be relaxed by introducing a marriage function, but we do not explore that here.
We also assume that males and females have the same survival and transition rates, and that male and female offspring are produced in equal proportions. These two assumptions imply that the male and female genotype × stage population vectors remain equal provided they start equal (de Vries and Caswell, 2019a). Therefore, we can treat the female vector as representative of both the male and the female populations, and we can calculate allele frequencies in the breeding population from the female population vector. These assumptions make it possible to model sexual reproduction in a one-sex model.
The matrices, vectors and mathematical operations used in this paper are listed in Table 1. Matrices are denoted by bold upper case letters (e.g., U) and vectors by bold lower case letters (e.g., n). Following Caswell et al. (2018), we use a tilde (e.g., Ũ , ñ ) to denote block structured matrices and vectors whose entries describe both demographic stage and genotype.
For a single locus with 2 alleles, say A and a, we will identify genotypes 1, 2, and 3 as AA, Aa, and aa, respectively. The population state vector is where e.g., n AA contains the numbers or densities of stages 1, …, s for genotype AA.
The population vector ñ is projected from time t to time t + 1 by a denisty-dependent matrix A ñ (˜), so that where Ũ describes survival and transition rates and F reproduction. The population projection matrix Ã depends on ñ because of genetics (the genotypes of offspring depend on gene frequencies of parents) and because of ecological nonlinearities due to density dependence. The genetic component of the model depends on the normalized population distribution vector, t p ( ), defined as where ||·|| denotes the 1-norm. The normalized population distribution vector consists of three genotype-specific population vectors:

The components of the population projection matrix
The population projection matrix is constructed from two sets of matrices describing demographic transitions for each genotype, and one set of matrices describing the parent-offspring map for each stage: map the genotype of a mother in stage j to the genotypes of her offspring. The (k, ℓ) entry of H j is the probability that the offspring of a genotype ℓ mother, of stage j, has genotype k. We assume that mating is random with respect to stage and hence that the parent-offspring map is the same for all stages, i.e. = H n H n (˜) (˜) j (assortative mating by stage would lead to differences among the H j ). The The vec operator, which stacks the columns of an m × n matrix X into a mn × 1 vector.
C. de Vries, et al. Ecological Modelling 416 (2020) 108875 matrix H n (˜) is discussed in Section 2.3. The model also formally contains matrices describing the transitions of individuals among genotype classes for each age or stage (Caswell et al., 2018), but since individuals do not change genotypes these are identity matrices.

Mating: from genotypes of parents to genotypes of offspring
Non-reproductive (e.g., immature) stages play no role in mating. Hence, following de Vries and Caswell (2019b), we define the breeding population by an indicator vector c j , for j = 1, …, g, that shows which stages of genotype j take part in mating. The ith entry of c j is 1 if individuals of stage i and genotype j reproduce, and 0 otherwise. The size of the breeding population is then where e i is a vector (g × 1) with a 1 in position i and zeros elsewhere, and ⊗ indicates the Kronecker product. Breeding stages are allowed to differ among genotypes in order to study the fate of traits that change reproductive schedules. In the special case where the genotypes do not differ in their reproductive stages, c i = c for all genotypes i and where T 1 g is a vector of ones of dimensions 1 × g. The genotype frequency vector within the breeding population is with E ii a matrix of dimension g × g with a 1 in the (i, i) location and zeros elsewhere. If the breeding vectors are the same for all genotypes, c i = c, then The genotype frequency vector for genotype i is (trivially) p i = e i . The gene frequencies in an individual of genotype i, q i , and the gene frequencies in the breeding population, q b , are functions of the genotype frequencies, so that where for the 2-allele case, Combining Eqs. (8) and (12) yields the following expression for the gene frequencies in the breeding population We set mutation rates to zero; see de Vries and Caswell (2019b) for details on how mutation could be included.

The matrix H
When a female randomly picks an allele from the gamete pool, she will pick allele A with probability q A b , and allele a with probability q a b .
These probabilities therefore determine the distribution of genotypes in her offspring, which is captured in the matrix H p (˜), The first column of H p (˜) contains the genotype distribution in the offspring of an AA mother; she produces an AA offspring with probability q A b and an Aa offspring with probability q a b . The second and third columns give the genotype distributions for mothers of genotypes Aa and aa. The matrix H is a homogenous function of degree zero of the population vector ñ since = H n H n (˜) (˜). Thus we can write it equally as a function of ñ or p . For a step by step derivation of H, see de Vries and Caswell (2019b, Section 2.3 and Appendix A).

The population projection matrix
To project the eco-evolutionary dynamics, the component matrices, U n (˜) i and F n (˜) i must be incorporated into the population projection matrix, A ñ (˜) (e.g., Caswell et al., 2018). To do so, create a set of blockdiagonal matrices , , and that contain the corresponding demographic matrices on the diagonal, i.e.
The fertility matrix F ñ (˜) is constructed from the block matrix containing genotype-specific fertility rates and from the parent-to-offspring genotype map, where K is the vec-permutation matrix (Henderson and Searle, 1981), which changes the arrangement of the vector from stages-within-genotypes to genotypes-within-stages. From right to left, the block-diagonal matrix first produces offspring, possibly of different birth stages (e.g., seedlings of different sizes) as a function of the genotype of the mother. When they appear, these offspring are associated with the genotype of the mother. The vec-permutation matrix K rearranges the vector, and then the block-diagonal matrix p (˜) allocates the offspring to their genotypes, based on the genotype of their mother and the genotype distribution of the rest of the population. Finally, K T returns the vector to its original arrangement. This vec-permutation approach to constructing matrix population models was introduced by Hunter and Caswell (2005) and is described in more detail by Caswell (2012) , gives the production of AA offspring by AA females that randomly pick allele A from the gamete pool, which happens with probability q A . The next block down, q F n (˜) a AA , contains the production of Aa offspring by an AA female picking allele a from the gamete pool. Similarly, the second and third row blocks contain offspring produced by an Aa female and an aa female, respectively.
Since individuals do not change their genotype once they are born, C. de Vries, et al. Ecological Modelling 416 (2020) 108875 the survival matrix is the block diagonal matrices, or written in terms of the genotype-specific block matrices, We project the stage × genotype dynamics with U ñ (˜) in Eq. (22) and F ñ (˜) in Eq. (20), using Eq. (3).

Stage × genotype dynamics of Tribolium
Flour beetles of the genus Tribolium have been used extensively to study population dynamics and population genetics, see Costantino et al. (2005) for a review. Tribolium is an economically important pest of flour and stored grain products. A nonlinear matrix population model for Tribolium was developed by Cushing and collaborators, see for example Cushing et al. (2002). The Tribolium model contains three stages: larvae (L), pupae (P), and adults (D). It is nonlinear because, in addition to feeding on flour, Tribolium adults cannibalize eggs and pupae; larvae in turn cannibalize eggs. These nonlinearities lead to a plethora of interesting bifurcations, attractors, and transient and asymptotic dynamics, which have been studied extensively (Costantino et al., , 1997(Costantino et al., , 2005Henson et al., 2002;Cushing et al., 2002;Edmunds et al., 2003).
We use the Tribolium model as the basis for a genetic model by including one-locus, two allele (A and a) Mendelian genetics.
We define effective adult densities D* and D † , and an effective larval density L † , as linear combinations of the genotype densities and use these to express the density effects of pupa cannibalism by adults (23) and egg cannibalism by adults (24) and larva (25). Each genotype i has a survival and transition matrix, where μ i and ν i are the larval and adult mortality rates, respectively. In addition, each genotype i has a fertility matrix, where β i is the fecundity at low densities. This parameterization permits selection to operate on any of the life-history characteristics; i.e. stagespecific viability, fertility, and/or cannibalism rates. Genotype × stage dynamics. As an example of genotype × stage dynamics, we use parameters estimated from a laboratory population of Tribolium by Dennis et al. (1995) (Table 2). We introduce a hypothetical allele with additive effects on fecundity (parameter β). The simulation was initialized with a population of AA individuals with low fecundity at the stage distribution of the homozygote equilibrium. After 50 iterations, one larval heterozygote is introduced into the population with a fecundity exactly in between the two homozygotes. The invading allele with larger birthrate gradually takes over the population and becomes fixed. As the genetic composition of the population changes ( Fig. 1C and D), the population structure changes ( Fig. 1A and B). Eventually, as the population approaches the aa boundary, the dynamics bifurcate from a stable equilibrium to a two-point cycle. Matlab code and parameters used for Fig. 1 are in the Online Supplementary Materials.

The outcome of density-dependent selection: conditions for genotype coexistence
The most basic question about selection, density-dependent or Table 2 Parameter estimates for malathion resistant and susceptible strains at 3 ppm malathion, from Table 2 in Cheung (2002). Heterozygote parameters are assumed to be identical to the resistant homozygote. Eco-evolutionary dynamics of the invasion of an allele with a higher fecundity β (β AA = 5.68, β Aa = 8.68, β aa = 11.68). All other parameters are taken from Table 1 in Dennis et al. (1995) and are set equal for all three genotypes: Abundance of larvae and pupae. c. Frequencies of the A and a alleles. d. Frequencies of the three genotypes.
C. de Vries, et al. Ecological Modelling 416 (2020) 108875 otherwise, is the question of whether genotypes coexist, so that the population retains some degree of genetic diversity, or whether one allele becomes fixed. The question becomes more complicated, but no less basic, when demographic structure and nonlinearity are included.
In this section, we present a general condition for determining the outcome of selection. The dynamics of ñ take place in an ωg-dimensional space defined by combinations of ω stages and g genotypes (with g = 3 in the present context). In the absence of mutation, the ω-dimensional boundary subspaces corresponding to the homozygous genotypes AA and aa are invariant under the dynamics specified by A ñ [˜]; these dynamics are given by the nonlinear projection matrices for the homozygous genotypes (Fig. 2).
We assume the existence of a single equilibrium on each boundary, and that this equilibrium is stable with respect to perturbations in the boundary subspace. Coexistence of the two alleles in a protected polymorphism (Levene, 1953;Prout, 1968;see Nagylaki, 1992, Chap. 6) results when the boundary subspaces are both unstable to perturbations into the interior. That is, if allele A can invade a population of aa individuals and allele a can invade a population of AA individuals, then a homozygote state can never be reached again once both alleles are present in the population (since both alleles grow when rare). Mutual invasibility therefore leads to a protected genetic polymorphism.
In general, the dynamics in the interior are unknown, and could include multiple equilibria, strange attractors, cycles, etc. Provided both boundary subspaces are unstable to perturbations into the interior, we refer to any of the (possibly exotic) dynamics in the interior as a protected polymorphism.
In general, nonlinear models could possess multiple invariant sets (equilibria, cycles, strange attractors) on the boundary. We restrict our discussion to models with a unique equilibrium on the boundary. Extending the analysis to include a k-cycle on the boundary would be accomplished by transforming the k-point cycle into an equilibrium of the k-fold map, but we do not consider this here.

Stability of the homozygote boundaries
The stability of a boundary equilibrium is determined by the spectral radius (i.e. the magnitude of the eigenvalue with the largest magnitude) of the Jacobian matrix at the equilibrium. We denote this eigenvalue by ζ AA and ζ aa for the AA and aa boundaries, respectively. If |ζ| < 1 the equilibrium is stable; if |ζ| > 1 it is unstable. We assume that the boundary equilibria are locally stable to perturbations within the boundary subspace, so if the equilibrium is unstable, the associated eigenvector must point into the interior, which implies that the invading allele increases when rare.
The Jacobian matrix, is obtained by differentiating Eq.
(3) and evaluating the resulting derivative at the boundary equilibrium. Here we give the Jacobian at the AA boundary; the expression at the aa boundary can be derived afterwards by symmetry. The Jacobian matrix at the AA boundary is see Appendix A for a derivation. The linearization reflects the two sources of nonlinearity in the model: those due to ecological density dependence and those due to the genetic frequency-dependence. The matrix M is a block structured matrix with blocks corresponding to genotypes and entries within the blocks corresponding to stages within genotypes, The block M 11 represents the contribution of perturbations in the AA direction to growth or decline of perturbations in the AA direction, block M 12 represents the contribution of perturbations in the Aa direction to growth or decline of perturbations in the AA direction, etc. All of the blocks in the Jacobian are given, with their derivation, in Appendix A.
Evaluated at the equilibrium on the AA boundary, the Jacobian matrix M is block upper triangular, with (See Eq. (A.35) in Appendix A). Thus the spectral radius of M depends on the eigenvalues of the diagonal blocks. Block M 33 = U aa projects perturbations in the aa direction, and since ρ(U aa ) < 1 this direction is always stable. The block M 11 projects perturbations within the AA boundary. By assumption, the boundary equilibrium is stable to such perturbations, so the spectral radius of M 11 must be less than one. The stability of the AA boundary equilibrium therefore depends on the where ρ(·) denotes the spectral radius of the matrix. By symmetry, the dominant eigenvalue of the Jacobian matrix at the aa boundary, denoted by ζ aa , is obtained by replacing AA by aa in (32) We note that the conditions for polymorphism in Eqs. (35) and (36) are a function of the nonlinear demographic rates of both the invading heterozygote and the resident homozygote (through the matrices U i and F i ) and of the structure of the homozygote equilibrium, n AA or n aa .

Tribolium revisited: density-dependent selection and pesticide resistance
Armed with these coexistence conditions, we return to Tribolium to study the spread of malathion resistance in the red flour beetle, Tribolium castaneum. Malathion was a commonly used pesticide in grain storage in the 1950s, and malathion resistance has since become widespread in Tribolium castaneum. The evolution of pesticide resistance is usually assumed to involve a fitness trade-off, in which resistant genotypes are at a disadvantage in the absence of the pesticide. However, there appears to be no fitness trade-off related to malathion resistance in T. castaneum. The resistant strain appears to have higher fitness even in the absence of malathion (Haubruge and Arnaud, 2001;Cheung, 2002;Arnaud et al., 2005). Arnaud et al. (2005) suggest that the higher fitness of malathion resistant genes may be the result of posterior modification of the insect genome after resistance became prevalent.
The genetics of malathion resistance varies among strains; sometimes resistance is found to be polygenic, while in other strains, it is due to a dominant allele at a single autosomal locus (Wool et al., 1982). Cheung (2002) studied a Tribolium strain in which resistance is primarily controlled by a single, dominant allele or closely linked set of alleles.
We denote the resistant allele with r and the susceptible allele with s; the genotypes are ss, rs, and rr. Since the resistant allele is almost completely dominant, Cheung assumed that the demographic rates of the rs genotype are identical to the demographic rates of the rr genotype; i.e.
for any population vector ñ . Cheung estimated U ss , F ss , U rr , and U rr in the laboratory under three levels of malathion exposure (0 ppm, 1.5 ppm, and 3 ppm). Experimental populations were kept in 120 ml Wheaton vials containing 20 grams of media (92.5% bleached white flour, 5% dry brewer's yeast, 1.5% ground fumigation, and 1% sunflower oil). The media containing 1.5 ppm and 3 ppm malathion was prepared by diluting malathion into the sunflower oil before adding to the flour. Each malathion treatment was initiated with 250 larva, 5 pupa, and 250 adults.
Three treatment groups were established, with different initial allele frequencies and four replicate populations in each treatment group. Populations in the two homozygous boundary treatments were initiated with all ss and rr insects. In the "evolving treatment" each population was initiated with all ss larvae and pupae, and with 245 ss and 5 rs adults, resulting in an initial r allele frequency of 0.01 among adults. All life stages in each population were counted once every two weeks for 80 weeks (t = 0, 1, 2, …, 40). Adult mortality was set at 50% by removing half the number of adults counted in the previous census minus the number of dead adults found in the current census. One of the replicate populations in the evolving treatment group was lost to disease at week 56 (t = 28). Parameters for rr and ss genotypes were estimated from homozygous populations with a maximum likelihood procedure, which is described in Section 2 of Cheung (2002). The parameter values are given in Table 2.  Table 2. Under these conditions, a homozygous ss population goes extinct (Fig. 3a), and a population of homozygous rr individuals persists at a stable equilibrium (Fig. 3c). When a small number of heterozygote rs individuals are introduced into a susceptible population, evolutionary rescue saves the population from extinction (Fig. 3b). 1 MATLAB code and parameters used for Fig. 3 are in the Online Supplementary Materials. Cheung (2002) estimated allele frequencies at the end of the experiment (t = 40) in the three surviving populations of the evolving treatment group by sampling and isolating 45-50 large larva or pupa, mating them with ss genotypes, and testing the survival of their offspring in malathion media to determine the genotypes of the sampled insects. This yielded estimates of the frequency of the r allele of 0.857, 0.948, and 0.889, with a mean and standard error of 0.898 ± 0.027. In comparison, at t = 40, the model predicts a frequency of the r allele of 0.862, which is not statistically different from the mean experimental value (p = 0.31).

Boundary stability
The stability of a population of the ss homozygote when invaded by the resistant allele is determined by the dominant eigenvalue ζ ss of the Jacobian. The final term in (35) simplifies significantly because the Tribolium model contains only one reproducing stage. If we denote the stage abundances at the susceptible equilibrium by L ss , P ss , D ss , then At 3 ppm malathion, ζ ss = 1.8450, so the resistant allele is able to invade. In the absence of malathion, ζ ss = 1.0265, so the resistant allele is superior, even without the advantage of the presence of malathion.
Since the model assumes complete dominance, the rs and rr individuals have identical parameters. This implies that ζ rr = 1, and that linearization fails to show the stability of the equilibrium. This is a well known phenomenon in models for selection against recessive alleles (e.g., Nagylaki, 1992, Sect. 4.2). Because recessive homozygotes are so rare close to the dominant equilibrium, elimination of the recessive is extremely slow. Simulations, however, show that the resistant boundary is stable to invasion by the susceptible allele.
Evolutionary stability analysis. The analytical expression that determines stability of the homozygote boundaries allows us to study the effect of parameter changes on stability; i.e. we can perform an evolutionary stability analysis. To demonstrate such an analysis, we explore the effects on stability of the degree of dominance of the resistant allele. This analysis was inspired by Beeman (1983) who found incomplete dominance of the resistant allele in the Rmal strain of Tribolium at high concentrations of malathion.
We analyze the effect of incomplete dominance by writing the heterozygote parameters as a convex combination of the parameters of the homozygotes; e.g., for mortality, and likewise for all other parameters in the model. The parameter x weighs the relative effect of the two alleles on all vital rates in heterozygotes. When x = 1, the s allele is dominant and when x = 0 the r allele is dominant. Fig. 4 shows the dominant eigenvalue ζ, evaluated at both boundary equilibria, as a function of x. The ss boundary is unstable in the range [0, 1); the rr boundary is stable in the range (0, 1]. Fig. 5 shows the effect of dominance on the fixation time, defined as the time required to reduce the frequency of the susceptible allele from 0.99 to 0.01 in the larvae. Alternative definitions of fixation time involving sums of stages rather than the larval stage alone result in the same qualitative shape. Using a different threshold to define allele fixation does effect the shape of the curve. The initial increase of the resistant allele is faster at x = 0 than at x = 0.5, but the final approach to the boundary is faster at x = 0.5. If we had defined competitive exclusion as the susceptible allele having a frequency less than 0.2 in larvae, the line in Fig. 5 would be monotonically increasing with x.

Discussion
Demographic models classify individuals into i-states, such as age classes or developmental stages. The i-state captures all the relevant information about an individual such that the fate of an individual depends only on its current i-state and the environment (Metz, 1977;Caswell and John, 1992;Caswell, 2001). In this paper, we extended the i-state to include the individual genotype. Treating genotype as a demographic state variable makes the powerful mathematical machinery  of matrix population models available for the study genotype × stage dynamics ( Figs. 1 and 3). Using genotype × stage models makes the equally powerful mathematical machinery of matrix calculus available for stability analysis via linearization (Sections 4 and 5). Using such a linearization, we obtained conditions for the coexistence of two alleles at one locus for a general density-dependent demographic model with age-or stage-structure. This opens new possibilities for eco-evolutionary analysis of life history traits. It complements, but is not intended to replace, the quantitative genetics approaches recently developed by Coulson et al. (2010).
We applied the model to study the genotype × stage dynamics of pesticide resistance in Tribolium castaneum. The model does an excellent job of describing experimental populations exposed to malathion, and successfully predicts the outcome of selection (fixation of the resistant allele). Fixation would occur more quickly for intermediate dominance than for complete dominance (Fig. 5).
The invasion speed of a pesticide-resistant allele in an agricultural pest is an important quantity for agriculturalists and/or policy makers. Pesticide resistance in insects is often determined by only one or two loci (Roush and McKenzie, 1987;Ffrench-Constant et al., 2004). For example, a number of single gene mutations are known to result in DDT resistance in Drosophila (Pittendrigh et al., 1997;Joussen et al., 2008). Similarly, a single gene was found in houseflies that influences the rates of penetration of DDT and dieldrin by Hoyer and Plapp (1968). There are many more examples, see reviews by Georghiou (1969), Roush andMcKenzie (1987), or Ffrench-Constant et al. (2004). The one locus, two allele model presented in this paper could therefore be applicable to many cases of insecticide resistance.
Maximization principles. Early work on density-dependent selection focused on the maximization of equilibrium population size (hence "Kselection"; MacArthur (1962); MacArthur and Wilson (1967); Roughgarden (1971)). Charlesworth (1994) extended this result to agestructured populations and showed that density-dependent selection maximizes the abundance of the age class that is exerting the densitydependent pressure on the population. For example, if adults cannibalize juveniles and juvenile mortality is a nonlinear function of adult density, then a successful invasion will always lead to a higher adult density. In deriving this result, Charlesworth assumes males and females have identical demographic rates; i.e. there is no sexual dimorphism in any life-history characteristics.
Because one cannot, in general, find expressions for equilibria in nonlinear models, we do not consider the maximization of density here. Moreover, we note that in a structured model there is no reason to assume that "density" is a single scalar quantity (e.g., Caswell et al., 2004). Each stage may have its own effects on other stages; a total density obtained by adding together tiny seedlings and big trees has little biological meaning. Equally, each vital rate may be influenced by a different set of stages, as in the Tribolium model where fertility is a function of larval and adult densities, but pupal survival depends only on adult density. Thus, maximization of "density" will not be as simple a concept in structured models as it is in unstructured models.
Extensions. The model presented here can be extended by relaxing a number of assumptions. Our assumption of female dominance leads to the assumption that male and female population vectors are proportional, so that the female population vector can be used to calculate gene frequencies in the (male) mating population. de Vries and Caswell (2019a) show that this assumption is met provided males and females have equal survival and transition probabilities, and are born at equal proportions. Both these assumptions could be relaxed in a two-sex model.
The model could also be extended to include more ecological interaction, such as time-dependent demographic rates, interactions among species, or dependence on environmental resources. The genetic component of the model can be expanded to include nonrandom mating, more than two alleles, or mutations.
At the cost of additional mathematics, our results could be extended to dynamics on the boundaries that are more exotic than fixed points. A k-point cycle, for example, can be transformed into an equilibrium by studying the k-fold map, i.e. by applying the population projection matrix k times.
We did not discuss frequency-dependent selection in this paper, in which the vital rates of an individual are a function of the frequencies of the genotypes in the population. Negative frequency-dependent selection, in which the fitness of a genotype declines as its frequency increases, is often invoked as a mechanisms for maintaining polymorphisms. In a structured population, the demographic rates could depend on the genotype frequencies in some subset of stages rather than in the entire population. In species with alternative mating strategies, the mating strategies often show negative frequency dependence (e.g., studies of salmon by Gross, 1985;Berejikian et al., 2010). Incorporating frequency dependence in our model would permit analysis of such traits.
The model construction introduced in this paper remains unchanged for frequency-dependent models, or for models with both frequency dependence and density dependence. However, a model with only frequency dependence becomes linear on the boundaries. Therefore the population is either exponentially growing or shrinking rather than at equilibrium on the boundary. To calculate the vulnerability of the homozygote population to invasion by the other allele in a frequencydependent model, it would be necessary to renormalize the population vector and to project the resulting frequency vector instead, as is discussed in de Vries and Caswell (2019b).
Conclusion. We began this paper by emphasizing how demographic and genetic processes combine to determine eco-evolutionary dynamics. At this point we have shown several examples of how Mendelian genetics and density-dependent vital rates, combined into a stage × genotype-structured matrix model, can be used to this end. The genotypes determine survival, transitions, and fertility. The rules of mating and segregation determine the genetic composition of offspring. The model can project dynamics of joint stage × genotype distributions, for hypothetical ( Fig. 1) and experimental (Fig. 3) situations. The eventual genetic coexistence is determined by stability analysis of the homozygous boundary equilibria.

Declaration of interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Derivation of Jacobian matrix, general nonlinear model
The local stability of a homozygote boundary equilibrium is determined by the largest eigenvalue (in absolute value) of the Jacobian matrix of the nonlinear matrix model evaluated at that equilibrium. We denote this eigenvalue at the AA boundary by ζ AA , and at the aa boundary by ζ aa . If the magnitude of the dominant eigenvalue of the Jacobian matrix is larger than one when evaluated at an equilibrium, then this equilibrium is unstable. The Jacobian matrix, is obtained by differentiating Eq. (3), and evaluating the resulting derivative at the boundary equilibrium. This requires a long series of matrix calculus operations (Magnus and Neudecker, 1985;Caswell, 2008). The analysis repeatedly takes advantage of the fact that n at the AA boundary contains zeros for the blocks corresponding to Aa and aa genotypes. Differentiate Eq. (A.2) to obtain where the explicit dependence of Ã on ñ has been omitted to avoid a cluttering of brackets. Multiply the second term by an ωg × ωg identity matrix, and apply the vec operator to both sides, remembering that as ñ is a vector, = n n vec˜˜, Next, apply Roth's theorem (Roth, 1934), )vec , to replace the vec operator with the Kronecker product: The matrix Ã can be decomposed into nine ω × ω block matrices, which are created by adding Eqs. (22) and (20): The blocks are denoted by A ij , so that, for example, C. de Vries, et al. Ecological Modelling 416 (2020) 108875 Armed with this expression for A vec˜, we analyze the term T t n I A (˜( ) )dvecg in (A.6). Replace the derivative of A vec˜with Eq. (A.14), such that where the differentials on both sides are evaluated at the boundary equilibrium. Finally, the First Identification Theorem Neudecker, 1985, 1988) and the chain rule together give the following formula for the Jacobian, This is the block-structured Jacobian matrix that appears in Eq. (30).