Evolutionary invasion analysis for structured populations: a synthesis

Evolutionary invasion analysis for structured populations: a synthesis Ryosuke Iritani1 and Troy Day2 1RIKEN Center for Interdisciplinary Theoretical and Mathematical Sciences, Saitama, Japan 2Department of Mathematics and Statistics, Queen’s University, Canada March 24, 2026 Abstract1 Natural populations exhibit complex class structures that profoundly shape evolutionary trajectories. While evolutionary2 demography provides a formal framework to predict adaptation using invasion fitness, the high mathematical dimensionality3 of these models often precludes analytical solutions, obscuring biological interpretation and hindering the analysis of long-4 term evolutionary outcomes. Because current reduction techniques remain fragmented, a unifying theoretical foundation is5 critically needed. Here, we introduce “structural evolutionary invasion analysis, ” a systematic framework that integrates two6 complementary tools to simplify complex life cycles. First, we formulate the “invasion determinant, ” an algebraic method that7 yields a direct scalar condition for mutant invasion. Second, we develop the Projected Next-Generation Matrix (PNGM), which8 structurally compresses life-cycle graphs by eliminating secondary classes. W e demonstrate that this reduction is mathemat-9 ically equivalent to separating dynamical timescales, explicitly preserving Fisher’s reproductive values for the retained focal10 classes. Crucially, under the standard assumption of weak selection, our synthesized framework guarantees that all properties11 of evolutionary singularities—including their location, convergence stability, and evolutionary stability—are strictly identical to12 those derived from the full, unreduced model. Illustrated with diverse ecological examples, this framework provides modellers13 with a rigorous and tractable toolkit for decoding state-dependent selection in high-dimensional populations.14 1 Introduction15 Natural populations exhibit complex structures, defined by state variables such as sex, age, body size, or habitat quality. This16 class structure profoundly shapes individual life histories, generating substantial heterogeneity in the selection pressures ex-17 perienced across different classes. Theoretical studies have revealed that such state-dependent selection is the primary driver18 behind a vast array of evolutionary phenomena, from conditional sex allocation and differential dispersal to facultative social19 behaviours and pathogen defences (e.g., Rodrigues & Gardner 2013; Rodrigues & Gardner 2015; Úbeda & Jansen 2016; Iritani20 & Cheptou 2017; Boots & Best 2018; Iritani et al. 2019; Kuijper & Johnstone 2019; Buckingham et al. 2023). Crucially, em-21 pirical research has increasingly corroborated these theoretical predictions, explicitly quantifying how selection gradients and22 life-history trade-offs vary dynamically across ages, physical conditions, or environmental patches in wild populations (e.g.,23 Coulson et al. 2005; Pelletier et al. 2007; Kruuk & Hill 2008). For example, specific trade-off structures have been shown to24 drive the evolution of antagonistic pleiotropy, with differential selection operating across distinct life-history stages (Guillaume25 & Otto 2012). Therefore, elucidating how class structure modulates natural selection is not merely a descriptive task, but a26 fundamental prerequisite for predicting evolutionary trajectories in realistic biological systems.27 T o mathematically capture this complex state-dependent selection, evolutionary dynamics in class-structured populations28 are predominantly analyzed through “evolutionary demography” , a powerful synthesis of adaptive dynamics (evolutionary in-29 vasion analysis) and matrix population models (Metz et al. 1992; Takada & Nakajima 1992; Dieckmann & Law 1996; Geritz30 et al. 1998; Caswell 2001; Rees & Ellner 2016). This analytical framework translates the demographic transitions of individuals31 across different classes into the invasion fitness (i.e., the asymptotic growth rate) of a rare mutant lineage. Crucially, success-32 ful invasion in structured populations depends not simply on immediate offspring production, but on the long-term, relative33 demographic contribution of each state, rigorously quantified as class-specific reproductive values (Fisher 1930; Taylor 1990;34 1 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint Caswell 2001; Rousset 2004). By explicitly linking demographic mechanics to these reproductive values, the evolutionary de-35 mography framework provides the essential theoretical machinery to decode how adaptive evolution proceeds amidst inherent36 population heterogeneity.37 One technical challenge is to obtain an analytically tractable expression for the invasion growth rate. This calculation typi-38 cally involves determining the eigenvalues of a Jacobian matrix, which describes the short-term growth of the mutant popula-39 tion. Y et, eigenvalues often lack closed-form expressions, and even when they are obtainable, their biological interpretation is40 frequently obscure. T o address this challenge, Hurford et al. (2010) developed the next-generation matrix approach building41 on earlier advances in mathematical epidemiology (Diekmann et al. 1990; van den Driessche & W atmough 2002). The core42 idea is to transform the Jacobian into an alternative matrix (the ‘next-generation matrix’) that captures the lifetime reproduc-43 tive output of individuals in each class. This next-generation matrix allows us to calculate the expected number of offspring44 produced over an individual’s lifetime and enables a systematic characterization of invasion conditions in the form of R0 > 1,45 where unity represents the threshold value for reproductive success. Because the invasion condition determines the long-term46 consequence of evolution, evaluating the form of R0 can yield crucial biological insights.47 Nevertheless, the next-generation approach faces its own limitations, particularly in systems with higher dimensionality.48 Even in relatively simple cases, the biological interpretation of the invasion condition remains obscure. For example, Massol &49 Débarre (2015) analyzed the next-generation matrix of evolutionary dynamics of dispersal in a metapopulation composed of50 two habitat types (good and poor quality). They derived the invasion condition as a solution to a quadratic equation, yet the51 biological interpretation of the solution remained obscure. Such problems become more pronounced in systems with greater52 dimensionality, where obtaining simple or interpretable invasion conditions is generally elusive.53 The mathematical dimensionality of these models scales directly with the number of population classes. Several studies54 have proposed methods to reduce the dimension of class-structured population models, although primarily in ecological or55 epidemiological contexts rather than evolutionary settings (Caswell 2001; de Camino Beck & Lewis 2007; Rueffler & Metz56 2013; Lewis et al. 2019). Rueffler et al. (2012) applied a linear algebraic approach to life-history evolution, illustrating the57 relevance of the optimization principle (Charnov 1976). However, their work did not offer a systematic way to derive biologically58 interpretable invasion conditions. Similarly, Roberts & Heesterbeek (2003) proposed the type-reproduction number approach,59 which allows for focusing on the spread of disease in a subset of compartments (classes) in the host population and thereby60 reducing the model dimension. However, despite its general utility, its application to evolutionary invasion analysis remains61 scarce.62 In this article, we formulate a unifying framework, which we term structural evolutionary invasion analysis, for deriving the63 invasion condition of a mutant in class-structured populations. W e present two complementary methods—the invasion deter-64 minant and the Projected Next-Generation Matrix (PNGM)—that are applicable to any class-structured model. W e highlight65 the pros and cons of these methods in the context of Fisher’ s theory of reproductive value (Fisher 1930; Taylor 1990; Caswell66 2001; Rousset 2004). Following Caswell (2001) and Rueffler et al. (2012), we also use life-cycle graphs to highlight the biological67 interpretation of these methods. Notably, these proposed methods can be seamlessly combined. W e illustrate their utility by68 applying them to four examples drawn from the literature. Although we assume weak selection (i.e., the phenotypic effect of69 mutation is small) throughout the manuscript, this synthesized framework provides modellers with powerful tools to simplify70 complex analyses and derive meaningful biological predictions.71 2 Overview: Evolutionary invasion analysis72 W e first review the next-generation matrix (NGM) methodology for evolutionary invasion analysis within the adaptive dy-73 namics framework (Hofbauer & Sigmund 1990; Dieckmann & Law 1996; Geritz et al. 1998; Hurford et al. 2010). In general,74 we begin by modelling nonlinear ecological dynamics of a class-structured, monomorphic population. Models may be cast75 in continuous time (ordinary differential equations) or discrete time (recursions). For now, we focus on the continuous-time76 dynamics and examine the discrete-time case later in the Examples section. The generic notation is summarized in Table 1.77 2 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint Table 1: Summary of notation used in the main text. Notation Definition Note K number of classes k, m class labels ∈ {1, 2, … , K} ∈ membership z evolving trait z for resident; z′ for mutant ′ mutant (superscript) ∘ resident/neutrality (superscript) δ phenotypic effect of mutation = z′ − z; assumed small Nk density in class k Nk for resident; Nk for mutants ⊤ transpose # ‌N vector of densities = (N1, … , NK)⊤ A rate of change in densities matrix with elements ak,m C or D continuous-time or discrete-time used for Jacobian label J Jacobian of mutant dynamics J D for discrete-time model σ largest real part of the eigenvalues spectral abscissa U inflow rate matrix (typically birth) subcomponent of J V remaining terms in J U − V = J F birth matrix need not be ‘birth’ S state-transition matrix need not be ‘transition’ G next-generation matrix (NGM) defined as UV −1 (or similar) gk,m (k, m)-element ofG gene transfer to class k from class m R0 reproduction number (generic) ρ(G); R0 > 1 implies invasion ρ largest absolute value of the eigenvalues spectral radius 𝕏 primary group of classes 𝕐 secondary group of classes 𝕏 ∪ 𝕐 = {1, 2, … , K}, 𝕏 ∩ 𝕐 = ∅ # ‌N𝕏, # ‌N𝕐 density vector of primary/secondary group 𝒫𝕏(G) Projected NGM (PNGM) NGM projected onto 𝕏 W e follow a standard adaptive dynamics approach and assume that the rate of mutation is small enough that the ecolog-78 ical dynamics reach an equilibrium while monomorphic (referred to as resident equilibrium). Biologically, we separate the79 timescales of ecological and evolutionary dynamics. W e also assume that the effect of mutations is small, meaning that the80 strength of natural selection is weak. Specifically, we suppose the resident equilibrium is monomorphic with trait value z; this81 trait is assumed to be a quantitative trait, possibly multidimensional. W e then consider a slightly different mutant phenotype82 z′ = z + δ, where δ is the small phenotypic effect of mutation.83 W e primarily focus on evaluating the invasion condition. This focus is grounded in the “invasion implies substitution”84 principle (Geritz 2005; Priklopil & Lehmann 2020): under the assumption of rare mutations, a mutant that satisfies the invasion85 condition is expected to fix in the population, driving the system to a new resident equilibrium. By iteratively assessing whether86 a mutant can increase in frequency when rare, we can therefore characterize the long-term evolutionary process of successive87 allelic substitutions without tracking the full frequency dynamics.88 3 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint Ecological dynamics89 Consider a population structured intoK discrete classes. These might be discrete habitats, ages, stages, sizes, or epidemiological90 states. For any pair of classes k and m, we assume that an individual of class m can contribute to class k; this contribution may91 be direct (via reproduction or class-transition) or indirect (via other classes). Such class structure is called irreducible (Caswell92 2001; Stott et al. 2010).93 T o understand the model structure, we use a life-cycle graph (Caswell 2001), where vertices (nodes) represent classes and94 directed edges (paths) represent positive contribution via class-transition or births. The irreducibility assumption then means95 that the life-cycle graph is strongly connected: for any pair of classes, (k, m), there exists at least one finite-length path linking96 the vertices k and m.97 The rate of change of class densities is modelled using a system of ordinary differential equations (ODEs):98 d # ‌N dt = A(z || # ‌N , z )# ‌N , (1) where # ‌N = (N1, … ,NK)⊤ collects the class densities. For any given resident trait value z, we assume there exists a unique,99 globally stable ecological equilibrium # ‌N ∘ = (N∘ 1 , … ,N ∘ K)⊤ with N ∘ k > 0 a for all classes, k = 1, 2, … ,K and A(z || # ‌N ∘, z ) = O100 (zero matrix). The potential existence of alternative equilibria or attractors does not affect our analysis, since the weak selection101 assumption prevents the population from converging to an alternative attractor (see Geritz et al. 2002; Priklopil & Lehmann102 2020). Throughout this article, we assume that the resident population dynamics are regulated by sufficient negative density103 dependence at high densities, and that the population can grow at low densities. Under these standard ecological conditions,104 the resident dynamics admit a strictly positive, globally stable equilibrium (Cushing 1998). W e refer to this globally stable105 monomorphic equilibrium as the resident equilibrium.106 Dynamics of a rare mutant107 Suppose a mutation arises with phenotype z′ = z + δ (where |δ| ≪ 1). The density of this mutant in class k, denoted by N ′ k,108 changes according to the following ordinary differential equations (ODEs):109 dN ′ k dt = K ∑ m=1 ak,m(z′ || # ‌N , # ‌N ′, z, z′ )N ′ m, k = 1, … ,K, (2) where ak,m represents the per-capita rate at which class-m mutants contribute to the class-k mutant population. This rate110 depends on the environmental feedback generated by the densities and traits of both the resident and the mutant. Accordingly,111 the resident dynamics (Equation (1)) must also be augmented to account for the mutant density.112 Assuming the mutant is initially rare, we can linearize the full joint system around the resident equilibrium, where the113 mutant is absent ( # ‌N ′ = # ‌0 ). In this regime, the ability of the mutant to invade is governed by the linearized system:114 d # ‌N ′ dt = J # ‌N ′. (3) Here, the Jacobian matrix J has entries jk,m ≔ ak,m(z′ || # ‌N , # ‌0 , z, z′ ), which can be interpreted as the net rate at which class- m115 mutants produce class-k mutants in the resident environment. The invasion condition is determined by the spectral abscissa of116 J, denoted by σ(J)(i.e., the largest real part among its eigenvalues). The mutant can successfully invade if σ(J)> 0, whereas it117 faces extinction if σ(J)< 0.118 The primary difficulty with this direct approach is that the spectral abscissa,σ(J), is often impossible to compute analytically.119 Although one might occasionally derive invasion conditions through ad-hoc algebraic manipulations, such as evaluating the120 sign of the determinant ofJ, these methods are typically highly model-specific. Furthermore, even when a closed-form invasion121 condition can be obtained, it frequently lacks a clear biological interpretation. Consequently, alternative, more systematic122 4 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint approaches have been developed.123 Most alternatives to directly analyzingσ(J)involve shifting the timescale of the evolutionary dynamics. Rather than tracking124 the mutant population over absolute time (e.g., days or years), these approaches evaluate population change over an organismal125 generation. This concept is best illustrated using a simple unstructured population. Suppose the density of a mutant,n′, changes126 according to the scalar ODE dn′ dt = un′ − vn′, where u and v are the constant per-capita birth and death rates, respectively, and127 time is measured in days. The spectral abscissa is simply σ = u − v, which has units of days −1. Invasion occurs if u − v > 0,128 meaning the daily fecundity exceeds the daily mortality. However, this inequality can be algebraically rearranged as uv−1 > 1.129 The left-hand side is now a dimensionless quantity, rendering it independent of the absolute time units. Biologically, because a130 mutant individual experiences a constant mortality rate v, its expected lifespan is v−1 days. Multiplying this expected lifespan131 by the daily offspring production rate, u, yields the dimensionless ratio uv−1, which represents the expected total number of132 offspring a mutant individual will produce over its entire lifetime. If this lifetime reproductive output exceeds unity, the mutant133 population will grow and successfully invade.134 While this unstructured example is simple enough that both invasion criteria (σ> 0 and uv−1 > 1) can be readily computed135 and biologically interpreted, introducing class structure significantly complicates the analysis. In multi-dimensional systems,136 the transition from continuous-time rates to generational reproductive output is far less straightforward, necessitating the for-137 mal mathematical machinery of the next-generation matrix.138 Next-generation matrix (NGM)139 The NGM method generalizes the concept of an individual’s expected lifetime reproductive output to structured populations.140 This is achieved by constructing a matrix,G, whose entries gk,m represent the expected number of class-koffspring produced by141 a single class-m individual over its lifetime. T o construct this matrix, we partition the continuous-time Jacobian asJ = U − V.142 Here, U is a non-negative matrix capturing the inflow of new individuals (e.g., “births”) into each class, while V captures143 all mortality and transition rates among the classes. From a mathematical standpoint, the biological labels assigned to these144 transitions are somewhat arbitrary, meaning the partitioning J = U − V is not necessarily unique. Although we will explore145 the implications of this non-uniqueness later, any valid partition must satisfy three fundamental mathematical requirements:146 (i) U ≥ 0, ensuring non-negative inflows; (ii) σ(−V) < 0, ensuring that in the absence of new inflows (U = O), the population147 densities asymptotically approach zero; and (iii) V −1 exists and is non-negative.148 Under these conditions, the next-generation matrix is defined as G ≔ UV −1. Hurford et al. (2010) formally demonstrated149 that150 σ(J) ⋛ 0 ⟺ ρ(G) ⋛ 1, (4) where ρ(G)denotes the spectral radius of the NGM (i.e., its largest eigenvalue in modulus), and the notation ⋛ indicates that151 the inequalities strictly correspond across the equivalence. Consequently, the invasion condition can be evaluated equivalently152 using eitherσ(J)or ρ(G). It is straightforward to verify that the earlier scalar ODE example represents a special, one-dimensional153 case of this general result.154 The primary advantage of the NGM framework is its capacity to yield simple, biologically interpretable invasion conditions155 that are often intractable when analyzing σ(J)directly. Nevertheless, in high-dimensional systems or models with complex156 life cycles, deriving clear analytical results from the standard NGM remains challenging. While alternative simplification tech-157 niques exist—such as graph-theoretical decompositions (Rueffler et al. 2012; Rueffler & Metz 2013)—they are often mathemati-158 cally demanding, which can limit their accessibility for applied modellers. A more direct strategy involves strategically choosing159 a non-standard decomposition J = U − V that still satisfies the mathematical requirements but yields a simpler NGM, thereby160 facilitating biological interpretation. However, this strategy currently lacks a unifying theory, typically forcing researchers to161 rely on ad hoc, model-specific derivations.162 In the following sections, we synthesize these disparate approaches into a unified framework that we term “structural evo-163 lutionary invasion analysis. ” This framework integrates two complementary tools: a generalized analysis of the invasion de-164 5 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint terminant of G, and a systematic methodology for decomposing the Jacobian J. W e call this approach “structural” because it165 exploits the algebraic architecture of the life cycle to characterize both short- and long-term evolutionary dynamics, entirely166 bypassing computationally demanding eigen-analyses. Ultimately, this framework serves as a systematic toolkit, empowering167 applied modellers to address complex evolutionary questions within the adaptive dynamics theory. T o apply this framework,168 we rely on three standard assumptions typical of adaptive dynamics: (i) natural selection is weak; (ii) the resident equilibrium169 is globally asymptotically stable; and (iii) the resulting NGM is irreducible.170 3 Synthesis: structural evolutionary invasion analysis171 The structural evolutionary invasion analysis synthesizes several distinct approaches. Specifically, we integrate two complemen-172 tary methods. The first method uses an algebraic approach (Appendix A): we calculate the determinant of matrices to capture173 the overall growth potential of a rare mutant lineage. This approach has remained underutilized in the literature, although it174 has appeared in previous works (e.g., Taylor & Bulmer 1980; Courteau & Lessard 2000; Day & Burns 2003; Metz & Leimar175 2011; Rueffler et al. 2012). W e refer to this quantity as the invasion determinant, highlighting both its mathematical nature as176 a determinant and its biological role in characterizing the invasion condition.177 Second, to systematically reduce the dimension, we apply the Nonnegative Decomposition Theorem (proven in Appendix178 B; Box 1). This theorem guarantees that any biologically motivated decomposition—provided it satisfies mild non-negativity179 conditions—yields an NGM that preserves the correct invasion threshold. Based on this theoretical foundation, we introduce180 the Projected NGM (PNGM) approach (Appendix C). The PNGM builds on the theory of the type-reproduction number devel-181 oped in theoretical epidemiology (Roberts & Heesterbeek 2003; Inaba & Nishiura 2008), focusing on the reproductive success182 of individuals in only a subset of classes while eliminating others. This reduces the model dimension and facilitates the calcula-183 tion of the invasion condition. While algebraically this amounts to choosing a specific decomposition J = U − V, biologically184 it can be interpreted as approximating the mutant dynamics by a separation of timescales: ‘slow’ classes are retained while ‘fast’185 classes are removed. From a life-cycle graph viewpoint, this effectively eliminates nodes from the original graph, as illustrated186 in Box 3. Crucially, this reduction preserves not only the invasion threshold but also Fisher’ s (1930) reproductive values of the187 retained classes, thus facilitating biological interpretation (Appendix C).188 A further key property is that these methods preserve the (i) position of the evolutionary equilibrium, (ii) convergence189 stability towards it (Appendix D), and (iii) evolutionary stability against invasion of alternative mutants (Appendix D). These190 properties provide useful tools for predicting long-term evolutionary outcomes, extending beyond the analysis of instantaneous191 invasion conditions.192 Key component 1: Invasion determinant193 Throughout our analysis, we assume that selection is weak, consistent with the adaptive dynamics framework. As alluded to194 above, it is possible to obtain a condition on the determinant of a matrix that provides an invasion threshold as follows:195 ρ(G)⋛ 1 ⟺ − det(I − G)⋛ 0 (5) (Appendix A). W e refer to − det(I − G)as the invasion determinant. The invasion determinant provides an algebraically196 closed-form expression for the invasion condition, which provides a tractable and general tool for arbitrary, class-structured,197 populations. Although the basic idea underlying this method appears in earlier studies (Taylor & Bulmer 1980; Courteau &198 Lessard 2000; Day & Burns 2003; Metz & Leimar 2011; Rueffler et al. 2012), it has received limited attention in the literature199 and has been rarely used. T o promote its practical utility, we here assign it a formal name and advocate its use as a systematic200 component of evolutionary invasion analysis.201 6 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint Incidentally, we can also derive the invasion conditions using the original, continuous-time Jacobian:202 ρ(G)⋛ 1 ⟺ − det(−J)⋛ 0 (6) (Appendix A), which has been also used in some studies (e.g., Day & Burns 2003; Hastings & Botsford 2006b; Boots & Best203 2018).204 One limitation of the invasion determinant approach is that it does not always yield biologically interpretable expressions.205 In particular, it does not account for Fisher’s (1930) reproductive value – a crucial quantity that helps biological interpretations206 of outcomes of evolutionary invasion analysis. Despite this limitation, the method significantly simplifies the invasion condi-207 tion and can thus facilitate theoretical analysis for evolutionary demography, with the advantage being more pronounced for208 complex or high-dimensional systems. Metz & Leimar (2011) further investigate the applicability of the invasion determinant209 and determine the condition for which we can relax the weak-selection assumption. This suggests that the invasion determinant210 can serve as a robust tool of deriving invasion condition beyond the traditional limit of adaptive dynamics theory.211 NGM: G = ( G𝕏,𝕏 G𝕏,𝕐 G𝕐,𝕏 G𝕐,𝕐 ) 𝒫𝕏 ⟼ 𝒫𝕏(G)= G𝕏,𝕏 + G𝕏,𝕐 (I − G𝕐,𝕐) −1 G𝕐,𝕏 ≅ ≅ life-cycle graph: Y2 X1 X2 Y1 Y3 Y4 𝒫𝕏 ⟼ Y2 X1 X2 Y1 Y3 Y4 ODE: d dt ( # ‌N ′ 𝕏# ‌N ′ 𝕐 ) d # ‌N ′ 𝕏 dt Figure 1: A schematic illustration of PNGM, with 𝕏 = {X1, X2} for primary and 𝕐 = {Y1, Y2, Y3, Y4} for secondary groups. The colors: red for the pathway [𝕏 ← 𝕏], orange for [𝕏 ← 𝕐], blue for [𝕐 ← 𝕏], and black for [𝕐 ← 𝕐]. Top panels depict linear- algebraic forms of the original next-generation matrix and PNGM. Bottom panels show the corresponding life-cycle graph (left) and its projection (right). Key component 2: Projected Next-Generation Matrix (PNGM)212 We refer to the next method as the projected NGM (PNGM) approach. We start with the above NGM G and partition the213 classes into two mutually exclusive groups: a ‘primary’ group and a ‘secondary’ group (Figure 1); for example, four classes can214 be divided into 𝕏 = {1, 2, 3} and 𝕐 = {4} (where {}is used for set notation; Figure 1). Such a grouping can be based on any215 biological feature of the classes.216 Symbolically, we block-partition the NGM as:217 G = ( G𝕏,𝕏 G𝕏,𝕐 G𝕐,𝕏 G𝕐,𝕐 ). (7) The two matrices in the diagonal blocks are square, and those in the off-diagonal blocks have compatible sizes. For exam-218 ple, when G𝕏,𝕏 is three-by-three and G𝕐,𝕐 is two-by-two, G𝕏,𝕐 is three-by-two and G𝕐,𝕏 is two-by-three. Biologically, each219 block matrix represents transfer of genes (either by reproduction or transition) among the primary and secondary groups. For220 7 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint example, G𝕏,𝕐 represents the transfer of genes from the secondary 𝕐 to primary 𝕏 group.221 Next, we consider how mutant individuals in the primary group (primary-group mutant) produce offspring in the primary222 group. There are two possibilities: (i) direct production, in which the primary mutants directly produce primary offspring223 during their lifetime. This process is captured by a loop-path [𝕏 ← 𝕏], or algebraically by the submatrix G𝕏,𝕏. The other224 possibility is (ii) indirect production, in which primary mutants produce secondary offspring, and these secondary offspring225 then produce primary offspring at some point in the future. This occurs along a path from 𝕏 to 𝕏 via 𝕐, [𝕏 ← 𝕐 ← 𝕏] . This226 process is represented by the matrix product G𝕏,𝕐G𝕐,𝕏. Following this logic, if we write Pn for the matrix whose elements give227 the total expected number of primary offspring produced by primary individuals within n generations (n ≥ 1), we get:228 Pn = G𝕏,𝕏 + G𝕏,𝕐 n−2 ∑ i=0 (G𝕐,𝕐) i G𝕐,𝕏 (8) This interpretation suggests that whenever taking the limit n → ∞ is possible, the following reduced-dimension matrix might229 be used to analyze invasion condition:230 𝒫𝕏(G)≔ lim n→∞ Pn = G𝕏,𝕏 + G𝕏,𝕐 (I − G𝕐,𝕐) −1 G𝕐,𝕏. (9) Most notably,Pn and G𝕏,𝕏 have the same size, and thus 𝒫𝕏(G)is the projection of G onto a lower dimensional matrix whose231 classes are given by the primary group 𝕏; the reader may note that while ‘projection’ in population demography often refers to232 an operator of population dynamics, here we use the term in the mathematical sense of ‘dimensionality reduction’ . In life-cycle233 graph, the projection corresponds to eliminating secondary groups and reallocating the paths (reproductive success) through234 secondary classes into primary classes (Figure 1). In other words, the projection method compresses the life-cycle graph.235 Since the matrix 𝒫𝕏(G)has a lower dimension, it is appealing to utilize this matrix as an alternative to the original NGM.236 Indeed, in Appendix C, we show that the limit always exists if selection is weak, and consequently prove the following equiva-237 lency:238 ρ(G)⋛ 1 ⟺ ρ(𝒫𝕏(G))⋛ 1. (10) The above result shows that the spectral radius of the reduced matrix 𝒫𝕏(G)provides an alternative invasion condition. W e239 refer to 𝒫𝕏(G)as projected next-generation matrix, or PNGM of G. This is referred to as the type-reproduction number in240 theoretical epidemiology (Roberts & Heesterbeek 2003).241 The PNGM may be also derived using the original NGM framework. W e first decompose the Jacobian as J = U − V, with242 block matrices J𝕏,𝕏 = U𝕏,𝕏 − V𝕏,𝕏, J𝕏,𝕐 = U𝕏,𝕐, J𝕐,𝕏 = U𝕐,𝕏, and J𝕐,𝕐 = U𝕐,𝕐 − V𝕐,𝕐, and write G = UV −1 for the resulting243 NGM, where we write G𝕏,𝕐 = U𝕏,𝕐V −1 𝕐,𝕐 etc.; we note that this decomposition satisfies the hypotheses of the theorem of NGM.244 On the other hand, we consider another decomposition and the resulting NGM as follows:245 ˜U = ( U𝕏,𝕏 U𝕏,𝕐 O O ), ˜V = ( V𝕐,𝕐 O −U𝕐,𝕏 −(U𝕐,𝕐 − V𝕐,𝕐) ), ˜G = ˜U ˜V −1 = ( 𝒫𝕏(G) ∗ O O ) (11) where the asterisk is some matrix. As a result, ρ(˜G)= ρ(𝒫𝕏(G)). Hence, PNGM can be regarded as an outcome of NGM with246 different decomposition of the Jacobian. Box 1 describes more technical details to clarify the NGMs that have the same invasion247 condition.248 Properties of the framework249 Beyond the invasion condition itself, the present framework possesses several, biologically useful properties (summarized in250 Box 3). First, the dimensionality reduction via the PNGM is mathematically equivalent to a separation of timescales, where the251 eliminated (secondary) classes are assumed to equilibrate instantaneously relative to the retained (primary) classes (Appendix252 C). Second, despite this reduction, the reproductive values of the retained classes are preserved (but only to the first order effect253 8 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint of selection; Appendix C). Finally, the framework preserves not only the invasion threshold (i.e., whether the spectral radius254 exceeds one) but also local stability properties around the equilibrium, thereby enabling the analysis of long-term consequence255 of evolutionary dynamics (e.g., convergence and evolutionary stability; Appendix D and E).256 Summary of the framework257 W e have provided key two methods: the invasion determinant and decomposition analysis. By exploiting the flexible decom-258 posability of the NGM and integrating the determinant formula with PNGM, we are able to generalize invasion analysis to259 higher-dimensional class-structured systems. The invasion determinant is an algebraic tool of writing down the invasion con-260 dition in closed form. The decomposition analysis allows us to partition NGM into two nonnegative matrices and calculate an261 alternative NGM based on the decomposition. The most notable application of the decomposition analysis includes the PNGM262 to calculate the type-reproduction number. These methods can be also used to locate candidate evolutionary equilibrium and its263 stability conditions (convergence- and evolutionary stability). One may preferably use any of the unifying methods. The struc-264 tural evolutionary invasion analysis may allow for systematic evaluation and prediction of evolutionary dynamics in structured265 populations. Some useful techniques are encapsulated in Box 2.266 In practice, the invasion determinant and PNGM serve complementary purposes, and the choice between them is best267 guided by the modelling goal. If the primary objective is analytical manipulation, differentiation, or stability calculations (e.g.,268 selection gradients, convergence stability, or evolutionary stability), the invasion determinant is often preferable because it yields269 a scalar condition that can be differentiated and simplified systematically. In contrast, if the goal is biological interpretation,270 PNGM is typically more useful because it retains a clear life-cycle meaning and preserves reproductive-value weighting on271 the retained (primary) classes, thereby clarifying which classes and pathways contribute most to mutant success. The main272 practical challenge in using PNGM is choosing primary and secondary groups. W e recommend selecting primary groups to273 match the biological question (typically the birth or transmission classes of interest) and treating as secondary those classes274 that act as transient intermediates. When multiple admissible partitions exist, one can choose the partition that yields the most275 interpretable or tractable reduced matrix, since the invasion threshold is invariant under these reductions.276 4 Examples277 W e demonstrate our methodology with four examples. More detailed analyses are encapsulated in Appendix E. The first two278 models use continuous time dynamics, while the last two use discrete time dynamics. W e briefly explain the procedure of279 invasion analysis for discrete time models in the third example.280 Example 1: Two-class model281 W e first examine a two-class model since it provides the simplest, non-trivial example. W e consider a stage-structured popula-282 tion with semi-adult and adult stages, both reproductive. The mutant dynamics read:283 dN ′ 1 dt = b′ 1N ′ 1 + b′ 2N ′ 2 − (a′ 1 + m1)N′ 1 dN ′ 2 dt = a′ 1N ′ 1 − m2N ′ 2, (13) where we have written parameters (each with a subscript) b′ for fecundity, which may be subject to density-dependent reduc-284 tion; a′ for transition; and m for mortality. Note that we have already substituted the resident equilibrium into the model.285 This example is relevant to models such as epidemiological susceptible-infected dynamics with host evolution, or source-sink286 metapopulation dynamics. W e can use the primed parameters a’ s orb’s for evolving trait(s), where these traits are traded-off287 against each other.288 9 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint The Jacobian of the mutant dynamics can be partitioned as:289 J = (b′ 1 b′ 2 0 0 ) ⏟ ⎵ ⎵ ⏟ ⎵ ⎵ ⏟ =B + ( 0 0 a′ 1 0) ⏟ ⎵ ⏟ ⎵ ⏟ =Tin − (a′ 1 0 0 0 ) ⏟ ⎵ ⏟ ⎵ ⏟ =Tout − (m1 0 0 m2 ) ⏟ ⎵ ⎵ ⏟ ⎵ ⎵ ⏟ =M , (14) where from left to right, the matrices represent birth, transition-in, transition-out, and mortality. W e write, dropping the290 prime symbol for brevity, f1 = b′ 1/(a′ 1 + m1)and f2 = b′ 2/m2 for the expected reproductive outputs in stage 1 and 2, and291 s = a′ 1/(a′ 1 + m′ 1)(< 1)for survival probability. Using these symbols, we can either decompose the Jacobian as (i) U = B,292 where the inflow term includes birth only; or (ii) U = B + T in, where the inflow term includes transitions. T o demonstrate the293 generality of the framework, we take U = B + T in − T out and V = M to get:294 G = (f1 − s f 2 s 0) ⏟ ⎵ ⎵ ⏟ ⎵ ⎵ ⏟ =H1 ((1 0 0 1 ) − (s 0 0 0 )) −1 ⏟ ⎵ ⎵ ⎵ ⎵ ⎵ ⎵ ⏟ ⎵ ⎵ ⎵ ⎵ ⎵ ⎵ ⏟ =(I−H2)−1 = ( f1−s 1−s f2 s 1−s 0) . (15) This matrix is not nonnegative when f1 < s (which occurs for instance when the class-1 individual does not reproduce), and295 may not be biologically relevant. However, this matrix has nonnegative off-diagonals (called as ‘Metzler matrix’ or ‘essentially296 nonnegative matrix’), and has similar properties to those of nonnegative matrices (Berman & Plemmons 1987). Indeed, taking297 PNGM with 𝕏 = {1} yields ρ(𝒫𝕏(G))= w1 ≔ (f1 + sf2 − s)/(1− s), which we can rearrange as the invasion condition w2 ≔298 f1 + sf2 > 1. Overall, we have essentially two formulas for invasion fitness: w1 = (f1 + sf2 − s)/(1− s),andw2 = f1 + sf2.299 More generally, the invasion condition for two-class models can be further simplified to:300 ρ(G)⋛ 1 ⟺ g1,1 + g2,2 − (g1,1g2,2 − g1,2g2,1)= tr(G)− det(G)⋛ 1, (16) which is useful for practical analyses for two-class models (e.g., Rodrigues & Gardner 2013; Massol & Débarre 2015; Iritani301 et al. 2019). There does not exist analogous formula when the system is three or higher dimensional.302 W e now examine the selection gradient, defined as the first derivative of the invasion fitness with respect to the mutant303 phenotypic deviation δ. While both w1 and w2 yield identical stability properties, we proceed with w2 due to its algebraic sim-304 plicity. Applying the differentiation technique outlined in Box 3, assuming thatf1, s, and f2 all depend on the mutant phenotype305 z′ = z + δ, yields:306 𝜕w2 𝜕δ || |δ=0 = f ∘ 1 ⋅ f (1) 1 f1 + s∘f ∘ 2 ⋅ (s(1) s + f (1) 2 f2 ). (17) Recalling that the neutrality condition implies w∘ 2 = f ∘ 1 + s∘f ∘ 2 = 1, the selection gradient can be rewritten by substituting307 s∘f ∘ 2 = 1 − f ∘ 1 :308 𝜕w2 𝜕δ |||δ=0 = f ∘ 1 ⋅ f (1) 1 f ∘ 1 + (1 − f ∘ 1 )⋅ (s(1) s∘ + f (1) 2 f ∘ 2 ). (18) This expression offers a clear biological interpretation: the total selection gradient is a weighted sum of the marginal changes309 in life-history traits. Specifically, the weights f ∘ 1 and 1 − f ∘ 1 correspond to the relative reproductive contributions of the class 1310 (direct reproduction) and class 2 (survival and reproduction) pathways to the next generation at the resident equilibrium. As311 such, the selection gradient captures the effect of natural selection throughout individual life-cycle. An analogous analysis is312 possible for the alternative fitness function w1, although the resulting expression may look more complex at first sight.313 While we focus here on the selection gradient (convergence stability), the second derivative required for assessing evolu-314 tionary stability can be derived in an equally straightforward manner using w2. This avoids the complexity of differentiating315 the eigenvalues of the original matrix (see Appendix D for the theoretical guarantee).316 T o summarize, this two-class example demonstrates that the structural evolutionary invasion analysis yields the same in-317 10 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint vasion condition across various types of biologically motivated decompositions of NGM. It also highlights the generality of318 the decomposition theorem, in which NGM need not be nonnegative and can be used as an alternative NGM with the same319 invasion condition.320 Example 2: Evolutionary epidemiological model321 The second, continuous-time example is based on the model of sex-structured epidemiological dynamics studied in Mitchell322 et al. (2022). While the authors have considered the coevolution between the virulence of vertically transmitting pathogen and323 sex-specific host recovery, we focus on host evolution for illustration.324 W e describe the mutant dynamics as:325 dS′ f dt = θffb′ f ⋅ S∘ m + I∘ m N ∘ T ⋅ (S′ f + (1− v)I′ f ) +θfmb′ m ⋅ S′ m + I′ m N ∘ T ⋅ (S∘ f + (1− v)I∘ f ) +γ′ fI′ f − hS′ f − μS′ f dI′ f dt = θffb′ f ⋅ S∘ m + I∘ m N ∘ T ⋅ vI′ f + θfmb′ m ⋅ S′ m + I′ m N ∘ T ⋅ vI∘ f + hS′ f − (μ + γ′ f + αf)I′ f dS′ m dt = θmfb′ f ⋅ S∘ m + I∘ m N ∘ T ⋅ (S′ f + (1− v)I′ f ) +θmmb′ m ⋅ S′ m + I′ m NT ⋅ (S∘ f + (1− v)I∘ f ) +γ′ mI′ m − hS′ m − μS′ m dI′ m dt = θmfb′ f ⋅ S∘ m + I∘ m N ∘ T ⋅ vI′ f + θmmb′ m ⋅ S′ m + I′ m N ∘ T ⋅ vI∘ f + hS′ m − (μ + γ′ m + αm)I′ m. (19) where S∘ f , I∘ f , S∘ m, and I∘ m are the resident equilibrium densities of susceptible female, infected female, susceptible male, and326 infected male, respectively. In this system of equations, we write b for birth (subscript f for maternity and m for paternity), as a327 function of recovery rate (fertility cost of recovery) and also the densities of the resident (density-dependent death and effective328 population size); θ for gene transmission between sexes (for example, all θs are one-half for diploid organisms); v for vertical329 transmission; γ for recovery; h = β(I∘ f + I∘ m)for the force of infection (with β transmission rate); μ for natural death; and α for330 virulence. Parameters are sex-specific except for (i) vertical transmission ( v) which occurs only from infected female and (ii)331 natural death rate μ which is common for both sexes.332 The Jacobian is a four-by-four matrix, with the (4,1)-st element (creation ofI′ m by S′ f) being zero. W e put only birth terms333 into U while all the others into V (including infection, recovery and mortality):334 J = ( Uff Ufm Umf Umm ) − ( Vf O O V m ), G = ( Uff Ufm Umf Umm ) ( V −1 f O O V −1 m ) = ( θffGf θfmGm θmfGf θmmGm ) (20) where Gf and Gm are fitness components of female and male, respectively, ignoring the genetic inheritance system. W e can335 immediately see that G has rank two (Mitchell et al. 2022, Eqn 10b), i.e., the model is effectively two-dimensional.336 W e now explore the formula of invasion fitness for diploidy and haplodiploigy separately. W e find that, for diploidy:337 G = 1 2 ( Gf Gm Gf Gm ). (21) Using the method in Box 1, we get the simple invasion condition ρ(Gf/2 + Gm/2) > 1. That is, essentially the NGM is given338 by the arithmetic mean of sex-specific NGMs. Note that Gf and Gm are both two-by-two, and so the two-class model re-339 sults are readily applicable to obtain the invasion condition. Meanwhile, the simple arithmetic-mean formula breaks down for340 haplodiploidy, leading to the invasion condition ρ(Gf/2 + GmGf/2) > 1 (Appendix E). Biologically, this multiplicative term341 (GmGf) explicitly captures the asymmetric nature of gene transmission: a male’ s reproductive success is strictly bottlenecked342 by the fitness of females in the previous generation. Consequently, this structural derivation immediately predicts a funda-343 mental asymmetry in intralocus sexual conflict; traits that benefit males at the expense of females (i.e., high Gm but low Gf)344 are strongly selected against. Our framework thus transparently demonstrates how genetic systems inherently bias sex-specific345 11 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint life-history evolution (such as sex-specific recovery or immunity) toward female optima, without the need to solve complex346 multi-generation recursions.347 Example 3: Discrete-time Lefkovitch model348 The above examples focus on continuous-time population dynamics and thus do not cover discrete-time models. Applications349 of discrete-time evolutionary dynamics include numerous traits – for example, dispersal (Hamilton & May 1977; Taylor 1988;350 Iritani & Iwasa 2014), mating system in plants (Lande & Schemske 1985; Iritani & Cheptou 2017), sex ratios (Hamilton 1967;351 Bulmer & Taylor 1980; Frank 1987; Iritani et al. 2021), and social behaviors (Taylor 1992; Rodrigues & Gardner 2012; Iritani352 2020). For discrete-time models, the population dynamics of the mutant linearized around the resident equilibrium is given353 by:354 # ‌N ′(t + 1)= J D # ‌N ′(t), (22) where J D represents the Jacobian (D for ‘discrete-time’), and also satisfies the mathematical conditions of a NGM. The invasion355 condition is thus given by ρ(J D)= ρ(G)> 1.356 W e use the result above to analyze Lefkovitch’s (1965) model, which helps us understand the biological meaning of PNGM.357 W e consider one juvenile stage (state 1) and three adult stages (states 2, 3 and 4), in which all adults are reproductively mature.358 W e first construct the NGM apply PNGM by partitioning the classes as𝕏 = {1, 2, 3} and 𝕐 = {4} (Figure 2):359 G ≔ ⎛ ⎜ ⎜ ⎜ ⎝ g1,1 g1,2 g1,3 g1,4 g2,1 g2,2 0 0 0 g3,2 g3,3 0 0 0 g4,3 g4,4 ⎞ ⎟ ⎟ ⎟ ⎠ ⟼ 𝒫 𝕏(G)= ⎛ ⎜⎜ ⎝ g1,1 g1,2 ˆg1,3 g2,1 g2,2 0 0 g3,2 g3,3 ⎞ ⎟⎟ ⎠ . (23) W e can immediately see that in this example, PNGM modifies the(1,3)-rd element asˆg1,3 = g1,3 + g1,4(1− g4,4)−1g4,3, where the360 second additional term represents the reproductive success of an individual in class 3 via class 4, assuming that the population361 dynamics in class-4 is at quasi-equilibrium (Figure 2). In general, PNGM allows the removal of intermediate classes from life-362 cycle graphs by appropriately redistributing their per capita contributions to other reproductive pathways (e.g., 3 to 1). This363 effectively compresses the life-cycle graph (i.e., matrix) while preserving the overall lifetime reproductive success.364 In Appendix E, we show that the invasion condition reads:365 ∑ m=1 fm m ∏ k=1 sk > 1, (24) where sk represents the survival probability of an individual in class k. Equation (24) represents the canonical form of repro-366 ductive success: fitness is survival times fecundity. Rueffler & Metz (2013) investigate the condition for which the reproductive367 success may be expressed in this canonical form (see Theorem 3 therein).368 T o further examine the long-term evolutionary outcomes, one may differentiate the left hand side of this invasion condition,369 using the techniques given in Box 2. The resulting formula can be interpreted as Fisher’ s (1941) average excess and thus yields370 a biologically transparent expression linking the effects of selection and pleiotropy (mutation influencing multiple phenotypes371 expressed at different stages).372 For ecological or epidemiological models, Hastings & Botsford (2006a), Rueffler et al. (2012), and Lewis et al. (2019) have373 rationalized this calculation by decomposing life-cycle graphs into sub-structures with distinct pathways of reproductive out-374 puts. These authors’ methods, however, require a careful treatment of signs (plus or minus. See Hastings & Botsford 2006a, Eqn375 2). On the other hand, applying the PNGM or invasion determinant automatically solves this problem because these calcula-376 tions only require solving linear equations or equivalently computing the determinant, which can be easily coded in standard377 computational software.378 12 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint G ≔ ⎛ ⎜ ⎜ ⎝ g1,1 g1,2 g1,3 g1,4 g2,1 g2,2 0 0 0 g3,2 g3,3 0 0 0 g4,3 g4,4 ⎞ ⎟ ⎟ ⎠ ≅ 4 1 2 3 g1,4 g2,1 g3,2g4,3 g1,2 g1,3 𝒫𝕏(G)= ( g1,1 g1,2 ˆg1,3 g2,1 g2,2 0 0 g3,2 g3,3 ) ≅ 4 1 2 3 g2,1 g3,2 g1,2 ˆg1,3 ( ˆg1,3 = g1,3 + g1,4 ⋅ (1 − g4,4) −1 ⋅ g4,3⏟ ⎵ ⎵ ⎵ ⎵ ⎵ ⏟ ⎵ ⎵ ⎵ ⎵ ⎵ ⏟ [1←4←3] ) Figure 2: Lefkovitch model. Top: the original NGM G and corresponding life-cycle graph. Self-loops are not shown for simplicity. Bottom: PNGM and resulting life-cycle graph. We see g1,3 is changed to g1,3 + g1,4 (1 − g4,4) −1 g4,3, which means that the reproductive pathway from class-3 to class-1 now accounts for the reproductive pathway via class-4. Example 4: Dispersal in stable habitats379 W e analyze the evolutionary model of dispersal in a saturated metapopulation (Hamilton & May 1977). The metapopulation380 consists of an infinitely many number of patches each occupied by the same number, K, of reproductive adults, following381 Wright’ s (1931) island-model of dispersal. For demography, we use the ‘birth-death Moran process’: all individuals produce the382 same number of offspring, and one individual is replaced by the newborn after natal dispersal of the offspring, where dispersal383 rate is d with mortality cost to dispersal c. Such replacement is based on random competition among offspring present in the384 patch to choose one individual, forming the next generation. All other offspring die.385 W e now consider the evolution of dispersal rate and introduce a mutant dispersald′ = d + δ under weak selection (i.e., δ is386 small). Writing k = 1, … ,K for the number of mutants in a patch, the transition probability S′ m,k from class k to m is given by:387 S′ m,k = ⎧⎪ ⎨⎪ ⎩ ϕ′ k = K − k K ⋅ k(1− d′) k(1− d′) + (K − k)(1− d) +K(1− c)d if m = k + 1, k ≠ 0; μ′ k = k K ⋅ (K− k)(1− d) +K(1− c)d k(1− d′) + (K− k)(1− d) +K(1− c)d if m = k − 1; 1 − ϕ′ k − μ′ k otherwise (25) (Ewens 2004), with conventionally ϕ′ 0 ≔ (1 − c)d′/(K(1− cd)), the probability that a patch with no mutant individual becomes388 a patch with 1 mutant. Using the transition probabilities in Equation (25) gives the discrete-time stochastic process, of:389 q′ 1(t + 1)= q′ 1(t)+ ϕ′ 0 K ∑ m=1 mq′ m(t)+ μ′ 2q′ 2(t)− (ϕ′ 1 + μ′ 1)q′ 1(t), ⋮ q′ k(t + 1)= q′ k(t)+ μ′ k+1q′ k+1(t)+ ϕ′ k−1q′ k−1(t)− (ϕ′ k + μ′ k)q′ k(t), ⋮ q′ K(t + 1)= q′ K(t)+ ϕ′ K−1q′ K−1(t)− μ′ Kq′ K(t), (26) 13 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint where q′ k(t)represents the proportion of patches with k mutants, conditioned that at least one mutant is present, at time t,390 k = 1, 2, … ,K. The invasion condition is given by ρ(G)= ∑K m=1 m ∏m k=1(ϕ′ k−1/μ′ k)(Mullon & Lehmann 2014), as is consistent391 with the metapopulation-fitness concept (Metz & Gyllenberg 2001; Ajar 2003; Massol et al. 2009; Mullon et al. 2016).392 W e recover the invasion condition by recursively applying the timescale separation to Equation (26). W e see that the dy-393 namics of class-K depends only on those of K and K − 1. W e therefore first solve q′ K(t + 1) = q′ K(t)to eliminate q′ K and get394 q′ K = q′ K−1ϕ′ K−1/μ′ K. Next, we solve q′ K−1(t + 1) = q′ K−1(t)to eliminate q′ K−1. Combining these results, we eliminate q′ K and395 q′ K−1. Recursively, we get q′ m = q′ 1 ∏m−1 k=1 (ϕ′ k/μ′ k+1)for m = 2, … ,K. Substituting these into the first line of Equation (26) yields:396 q′ 1(t + 1)− q′ 1(t)= ϕ′ 0 K ∑ m=1 mq′ 1(t) m−1 ∏ k=1 ϕ′ k μ′ k+1 − μ′ 1q′ 1(t)= ( K ∑ m=1 m m ∏ k=1 ϕ′ k−1 μ′ k − 1)μ′ 1q′ 1(t). (27) Thus the invasion condition is given by ∑K m=1 m ∏m k=1(ϕ′ k−1/μ′ k) > 1 as desired. This invasion condition is similar to Equa-397 tion (24) in the Lefkovitch model, due to the product term representing the sojourn time in each state. W e can also derive the398 explicit formula of the selection gradient by evaluating the number of generations a patch spends possessing k mutants until399 the extinction of the mutant lineage (Mullon et al. 2016).400 This example demonstrates that recursive operations of the quasi-equilibrium approximation may simplify the problem401 in some cases. The quasi-equilibrium approximation is especially useful when NGM is tri-diagonal, which occurs for stage-402 structured models or Moran processes. However, using other models (say Wright-Fisher process) would make the analysis403 more complex (Ajar 2003). Y et, it often suffices to examine the first or second derivative of the invasion condition. For this404 purpose, we can make use of properties of matrix derivatives. More details are presented in Box 2 and Appendix D.405 (A) Unidirectional (B) Independence N1 N2 N3N3 N1 N2 N3 N4 Figure 3: Two typical examples of reducible life-cycle graphs (i.e., non-irreducible life-cycles). Arrows represent transitions between classes, and gray bands highlight subsets of classes that form irreducible components — subgraphs in which every node is reachable from any other within the same component. In both examples, the invasion-implies-substitution principle fails: a mutation arising in class 3, for instance, cannot spread to fix in the entire population in both panels: (A) Classes 1 and 2 are interconnected, allowing gene flow between them. However, no gene flow from class 3 to class 2 is possible. Thus a mutant emerging in class 3 will never get fixated. (B) Classes 1 and 2 form one interconnected group, and classes 3 and 4 form another. There is no gene flow between these two groups, resulting in two isolated irreducible components. Again, fixation is impossible. 5 Discussion406 We have synthesized a general, unifying framework of evolutionary invasion analysis for class-structured populations. This407 synthesis integrates two main approaches, the invasion determinant and the projected next-generation matrix (PNGM). The408 invasion determinant provides a determinant-based formula for invasion fitness (Taylor & Bulmer 1980). We also show that the409 PNGM method has biologically clear interpretations (Box 3): (i) Fisher’s (1930) reproductive values remain invariant under410 this transformation, and (ii) PNGM results from a separation of dynamical timescales among classes. These advantages are411 illustrated in four examples. Beyond invasion conditions, these methods also allow us to derive long-term, stability properties412 of evolutionary equilibrium (Box 3). Illustrative examples demonstrate the utility of this framework, with practical protocols for413 two-class models (Example 1). Example 2 of evolutionary epidemiology highlights that even when the full Jacobian is complex414 (four-by-four), the structural approach can decompose the invasion fitness into interpretable sex-specific components. The415 14 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint framework is readily applicable to both continuous-time and discrete-time dynamics (Example 3), including spatially structured416 model (Example 4). Collectively, the structural evolutionary invasion analysis provides a systematic tool for analyzing invasion417 conditions within the adaptive dynamics theory.418 The invasion determinant, originally formulated by Taylor & Bulmer (1980) (see also Courteau & Lessard 2000; Rueffler419 et al. 2012), provides a general algebraic tool to derive invasion conditions. Our derivation relies crucially on the assumption of420 weak selection, where the phenotypic effect of a mutation is small. This assumption, however, is not highly restrictive; it often421 provides a robust approximation for evolutionary dynamics (Mullon & Lehmann 2014) and can be formally relaxed, particu-422 larly for continuous quantitative traits (for precise details, see Metz & Leimar 2011). Furthermore, the invasion determinant423 facilitates the derivation of higher-order properties, such as the stability of evolutionary equilibria, and serves as the mathe-424 matical foundation for the negative decomposition theorem (which yields the projected next-generation matrix). Thus, while425 the bare determinant might sometimes lack immediate biological transparency, it constitutes the analytical core of structural426 evolutionary invasion analysis.427 The projected next-generation matrix (PNGM) systematically simplifies life-cycle graphs by eliminating secondary classes.428 While this method builds upon the type-reproduction number theory from epidemiology (Roberts & Heesterbeek 2003),429 our synthesis highlights two crucial biological insights for evolutionary demography (Box 3). First, constructing a PNGM430 is mathematically equivalent to separating the timescales of class dynamics. This equivalence is particularly advantageous in431 high-dimensional systems—as demonstrated in the metapopulation model (Example 4)—because it reduces complex matrix432 problems to tractable scalar or low-dimensional conditions. It allows researchers to place any subset of classes into a quasi-433 equilibrium state without requiring an actual, physiological separation of timescales. As illustrated by the Lefkovitch model434 (Example 3), removing intermediate stages transparently exposes the net reproductive pathways. Second, the PNGM preserves435 Fisher’ s (1930) reproductive values for the retained classes. In structured populations, variations in how different classes trans-436 mit alleles (nonheritable class transmission) complicate the evaluation of allele-frequency changes (Taylor 1990; Frank 1998;437 Rousset 2004; Priklopil & Lehmann 2020, 2024). By employing this structural reduction, we demonstrate that tracking gene-438 frequency changes strictly within the focal subset is sufficient to determine mutant invasibility. This implies that transient class439 transmission dynamics can be safely bypassed when evaluating initial invasion. Therefore, applying the PNGM to evolutionary440 dynamics provides conceptual clarity, not merely a mathematical shortcut.441 The theory of evolutionary demography aims to understand how complex demographic structures shape patterns of adapta-442 tion (Takada & Nakajima 1992; Cushing 1998; Takada & Nakajima 1998; Easterling et al. 2000; Caswell 2001; T uljapurkar et al.443 2003). Over the past decade, integral projection models (IPMs) have emerged as powerful tools in this field (Ellner et al. 2016;444 Rees & Ellner 2016). IPMs incorporate continuous traits, such as size, weight, or age, into the adaptive dynamics framework,445 allowing for a realistic representation of individual heterogeneity. These models typically use an integral operator rather than a446 finite matrix for the NGM, resulting in an infinite-dimensional system. Nevertheless, the core logic of the PNGM framework447 remains robust for two reasons. First, the corresponding task of the PNGM in IPMs is simply to partition the physiological448 space into two regions and compute an effective operator for the primary region; this concept has already been integrated into449 epidemiological studies (Inaba & Nishiura 2008). Second, in practice, IPMs are numerically implemented by discretizing the450 continuous state space into a finite number of classes. This numerical scheme immediately reduces the problem to the stan-451 dard, matrix-based domain. Thus, our approach offers a principled way to extend evolutionary invasion analysis into the IPM452 framework, potentially improving the tractability and interpretability of extended IPMs.453 Beyond lifetime reproductive success, the present framework facilitates the determination of long-term evolutionary con-454 sequences, namely convergence stability and evolutionary stability (Maynard Smith & Price 1973; Maynard Smith 1982; Chris-455 tiansen 1991; Takada & Kigami 1991; Eshel 1996). Such stability analyses require the first and second partial derivatives of456 the invasion fitness with respect to phenotypes. Difficulty typically arises when the population comprises many classes, as the457 analytical differentiation of the leading eigenvalue becomes computationally prohibitive and requires explicit information on458 the reproductive values of all classes (e.g., Avila & Mullon 2023). In this context, our framework offers two complementary459 pathways to overcome this high-dimensional hurdle. For algebraic or computational purposes, one can directly differentiate460 15 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint the invasion determinant (Appendix D). This approach is often more straightforward than differentiating the spectral radius it-461 self, as it avoids the explicit calculation of eigenvectors. Alternatively, when interpretability is paramount, the PNGM approach462 allows one to reduce the dimensionality prior to differentiation. This structurally restricts the necessary reproductive value463 calculations to only the focal classes (Appendix C), thereby simplifying the stability conditions while retaining the biological464 intuition of the life cycle. By choosing between these algebraic and reductionist approaches, researchers can flexibly tackle the465 stability analysis of complex models.466 The present framework is also applicable to ecological models, particularly for analyzing the spread of invasive species.467 For example, Hastings & Botsford (2006a,b) used analogous equations to investigate how spatial structure influences inva-468 sion dynamics. More recently, building on the target reproduction number developed by Lewis et al. (2019), Harrington et al.469 (2023) examined source-sink metapopulation dynamics in sea lice (Lepeophtheirus salmonis), a parasitic species threatening470 the salmon farming industry. The next-generation matrix approach has proven highly effective for evaluating persistence in471 metapopulations with well-characterized dispersal networks. In this context, the PNGM offers a mathematically grounded472 strategy to compress complex network structures into interpretable forms, enhancing both biological understanding and man-473 agement decisions. However, this framework applies strictly when the invader is rare. Mathematically, the geometric series474 used to compute the PNGM must converge, which is guaranteed only when the spectral radius of the secondary sub-process475 is strictly less than unity. Once the species becomes widespread and density-dependent effects alter the transition rates, these476 initial approximations no longer hold, emphasizing the need for early detection in invasion management.477 Having presented multiple fitness proxies that share the same invasion condition, a practical question arises: which method478 should one use? In practice, the invasion determinant and the PNGM serve complementary purposes. If the primary objective479 is algebraic manipulation or stability calculations (e.g., deriving selection gradients), the invasion determinant is often prefer-480 able. It yields a scalar criterion that can be differentiated systematically without requiring eigenvectors or an explicit partitioning481 of classes. Conversely, if the goal is biological interpretation, the PNGM is typically superior. It retains a clear life-cycle interpre-482 tation and preserves the reproductive-value weighting on the retained classes, explicitly clarifying which pathways contribute483 most to mutant success. The main practical challenge when using the PNGM is selecting the primary and secondary groups.484 W e recommend defining the primary group based on the specific biological question (e.g., birth or key transmission events)485 and treating the intermediate classes as secondary. When multiple admissible partitions are available, one should select the486 partition that yields the most interpretable reduced proxy, as the fundamental invasion threshold remains invariant across all487 choices.488 Finally, we outline the underlying ecological assumptions for applying this framework. First, we assume that the resident489 population dynamics are regulated by sufficient density dependence to reach a unique, strictly positive, and globally stable equi-490 librium. A standard local stability analysis is typically sufficient to ensure that the population does not depart from the resident491 equilibrium once established. Formally, ecological equilibria often change only slightly following recurrent gene substitutions492 (Geritz et al. 2002). While bistability or multiple attractors can complicate the analysis, such systems can generally be handled493 by analyzing the invasibility of mutants at each stable equilibrium separately.494 Second, we assume that the transition matrix of the population model is irreducible (Stott et al. 2010). Biologically, this495 means that for any pair of distinct classes (k, m), there is at least one sequence of transitions from m to k. While seemingly re-496 strictive, this assumption aligns naturally with the core premises of adaptive dynamics. In reducible models, gene flow between497 certain classes is permanently blocked, which contradicts the invasion-implies-substitution principle (Figure 3). For instance,498 if a system has a one-way transition from class 2 to class 3 with no return path (Figure 3A), a mutant arising in class 3 has499 zero chance of fixation regardless of its selective advantage. Similarly, a population perfectly split into two disconnected sub-500 groups (Figure 3C) functions as two independent evolutionary entities and must be analyzed as such. Thus, the irreducibility501 assumption simply ensures that the population forms a single, cohesive evolutionary unit.502 T o conclude, we have unified the methodology of evolutionary invasion analysis for class-structured populations under the503 framework of structural evolutionary invasion analysis. By integrating the invasion determinant, the arbitrary decomposability504 of the next-generation matrix, and the projected next-generation matrix, we provide researchers with flexible tools that can505 16 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint be used interchangeably or in combination. Furthermore, this framework naturally complements integral projection models,506 paving the way for the rigorous analysis of infinite-dimensional evolutionary demography.507 References508 Ajar, É. (2003). Analysis of disruptive selection in subdivided populations. BMC Evol. Biol. , 3.1, p. 22. DOI: 10.1186/1471-2148-3-22.509 Avila, P. & Mullon, C. (2023). Evolutionary game theory and the adaptive dynamics approach: adaptation where individuals interact. Phil.510 Trans. R. Soc. B, 378.1876. DOI: 10.1098/rstb.2021.0502.511 Berman, A. & Plemmons, R. J. (1987). Nonnegative matrices in the mathematical sciences. Classics in applied mathematics 9. Society for Industrial512 and Applied Mathematics.513 Boots, M. & Best, A. (2018). The evolution of constitutive and induced defences to infectious disease. Proc. R. Soc. B, 285.1883. DOI: 10.1098/514 rspb.2018.0658.515 Buckingham, L. J., Bruns, E. L., & Ashby, B. (2023). The evolution of age-specific resistance to infectious disease. Proc. R. Soc. B, 290.1991. DOI:516 10.1098/rspb.2022.2000.517 Bulmer, M. G. & Taylor, P. D. (1980). Dispersal and the sex ratio. Nature, 284.5755, pp. 448–449. DOI: 10.1038/284448a0.518 Caswell, H. (2001). Matrix population models. Wiley Online Library.519 Charnov, E. L. (1976). Optimal foraging, the marginal value theorem. Theor. Popul. Biol., 9.2, pp. 129–136. DOI: 10.1016/0040-5809(76)90040-520 x.521 Christiansen, F. B. (1991). On conditions for evolutionary stability for a continuously varying character. Am. Nat., pp. 37–50. DOI: 10.1086/522 285203.523 Coulson, T., Benton, T., Lundberg, P., Dall, S., Kendall, B., & Gaillard, J.-M. (2005). Estimating individual contributions to population growth:524 evolutionary fitness in ecological time. Proc. R. Soc. B, 273.1586, pp. 547–55. DOI: 10.1098/rspb.2005.3357.525 Courteau, J. & Lessard, S. (2000). Optimal sex ratios in structured populations. J. Theor. Biol. , 207.2, pp. 159–175. DOI: 10.1006/jtbi.2000.526 2160.527 Cushing, J. M. (1998). An Introduction to Structured Population Dynamics . Society for Industrial and Applied Mathematics. DOI: 10.1137/1.528 9781611970005.529 Day, T. & Burns, J. G. (2003). A consideration of patterns of virulence arising from host-parasite coevolution. Evolution, 57.3, pp. 671–676.530 DOI: 10.1111/j.0014-3820.2003.tb01558.x.531 de Camino Beck, T. & Lewis, M. A. (2007). A new method for calculating net reproductive rate from graph reduction with applications to the532 control of invasive species. Bull. Math. Biol., 69.4, pp. 1341–1354. DOI: 10.1007/s11538-006-9162-0.533 Dieckmann, U. & Law, R. (1996). The dynamical theory of coevolution: a derivation from stochastic ecological processes. J. Math. Biol., 34.5-6,534 pp. 579–612. DOI: 10.1007/s002850050022.535 Diekmann, O., Heesterbeek, J. A. P., & Metz, J. A. J. (1990). On the definition and the computation of the basic reproduction ratio R0 in models536 for infectious diseases in heterogeneous populations. J. Math. Biol., 28.4. DOI: 10.1007/bf00178324.537 Easterling, M. R., Ellner, S. P., & Dixon, P. M. (2000). Size-specific sensitivity: applying a new structured population model. Ecology, 81.3,538 pp. 694–708. DOI: 10.1890/0012-9658(2000)081[0694:sssaan]2.0.co;2.539 Ellner, S. P., Childs, D. Z., & Rees, M. (2016). Data-driven modelling of structured populations . Springer.540 Eshel, I. (1996). On the changing concept of evolutionary population stability as a reflection of a changing point of view in the quantitative541 theory of evolution. J. Math. Biol., 34.5–6, pp. 485–510. DOI: 10.1007/bf02409747.542 Ewens, W. J. (2004). Mathematical Population Genetics: I. Theoretical Introduction. 2nd ed. Interdisciplinary Applied Mathematics 27. Springer-543 Verlag New York.544 Fisher, R. A. (1930). The genetical theory of natural selection . Clarendon Press. DOI: 10.5962/bhl.title.27468.545 — (1941). Average excess and average effect of a gene substitution. Ann. Eugen., 11.1, pp. 53–63. DOI: 10 . 1111 / j . 1469 - 1809 . 1941 .546 tb02272.x.547 Frank, S. A. (1998). Foundations of social evolution. Princeton University Press.548 — (1987). Demography and sex ratio in social spiders. Evolution, 41.6, pp. 1267–1281. DOI: 10.1111/j.1558-5646.1987.tb02465.x.549 Geritz, S. A. H., Kisdi, É., Meszéna, G., & Metz, J. (1998). Evolutionarily singular strategies and the adaptive growth and branching of the550 evolutionary tree. Evol. Ecol., 12.1, pp. 35–57. DOI: 10.1023/a:1006554906681.551 17 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint Geritz, S. A. H. (2005). Resident-invader dynamics and the coexistence of similar strategies. J. Math. Biol. , 50.1, pp. 67–82. DOI: 10 . 1007 /552 s00285-004-0280-8.553 Geritz, S. A. H., Gyllenberg, M., Jacobs, F. J., & Parvinen, K. (2002). Invasion dynamics and attractor inheritance. J. Math. Biol., 44.6, pp. 548–560.554 DOI: 10.1007/s002850100136.555 Guillaume, F. & Otto, S. P. (2012). Gene functional trade-offs and the evolution of pleiotropy. Genetics, 192.4, pp. 1389–1409. DOI: 10.1534/556 genetics.112.143214.557 Hamilton, W. D. (1967). Extraordinary sex ratios. Science, 156.3774, pp. 477–488. DOI: 10.1126/science.156.3774.477.558 Hamilton, W. D. & May, R. M. (1977). Dispersal in stable habitats. Nature, 269.5629, pp. 578–581. DOI: 10.1038/269578a0.559 Harrington, P. D., Cantrell, D. L., & Lewis, M. A. (2023). Next-generation matrices for marine metapopulations: The case of sea lice on salmon560 farms. Ecol. Evol., 13.4. DOI: 10.1002/ece3.10027.561 Hastings, A. & Botsford, L. W. (2006a). A simple persistence condition for structured populations. Ecol. Lett., 9.7, pp. 846–852. DOI: 10.1111/562 j.1461-0248.2006.00940.x.563 — (2006b). Persistence of spatial populations depends on returning home. Proc. Natl. Acad. Sci. USA, 103.15, pp. 6067–6072. DOI: 10.1073/564 pnas.0506651103.565 Hofbauer, J. & Sigmund, K. (1990). Adaptive dynamics and evolutionary stability. Appl. Math. Lett., 3.4, pp. 75–79. DOI: 10.1016/0893-9659(90)566 90051-c.567 Hurford, A., Cownden, D., & Day, T. (2010). Next-generation tools for evolutionary invasion analyses. J. R. Soc. Interface , 7.45, pp. 561–571.568 DOI: 10.1098/rsif.2009.0448.569 Inaba, H. & Nishiura, H. (2008). The state-reproduction number for a multistate class age structured epidemic system and its application to570 the asymptomatic transmission model. Math. Biosci., 216.1, pp. 77–89. DOI: 10.1016/j.mbs.2008.08.005.571 Iritani, R. & Cheptou, P.-O. (2017). Joint evolution of differential seed dispersal and self-fertilization. J. Evol. Biol. , 30.8, pp. 1526–1543. DOI:572 10.1111/jeb.13120.573 Iritani, R. (2020). Gametophytic competition games among relatives: when does spatial structure select for facilitativeness or competitiveness574 in pollination? J. Ecol., 108.1, pp. 1–13. DOI: 10.1111/1365-2745.13282.575 Iritani, R. & Iwasa, Y. (2014). Parasite infection drives the evolution of state-dependent dispersal of the host. Theor. Popul. Biol. , 92, pp. 1–13.576 DOI: 10.1016/j.tpb.2013.10.005.577 Iritani, R., Visher, E., & Boots, M. (2019). The evolution of stage-specific virulence: Differential selection of parasites in juveniles. Evol. Lett., 3.2,578 pp. 162–172. DOI: 10.1002/evl3.105.579 Iritani, R., West, S. A., & Abe, J. (2021). Cooperative interactions among females can lead to even more extraordinary sex ratios. Evol. Lett., 5.4,580 pp. 370–384. DOI: 10.1002/evl3.217.581 Kemeny, J. G., Snell, J., & Knapp, A. (1966). Markov chains. New York, NY: University Series in Higher Mathematics, Van Nostrand.582 Kruuk, L. & Hill, W. (2008). Evolutionary dynamics of wild populations: the use of long-term pedigree data. Proc. R. Soc. B, 275.1635, pp. 593–583 596. DOI: 10.1098/rspb.2007.1689.584 Kuijper, B. & Johnstone, R. A. (2019). The evolution of early-life effects on social behaviour—why should social adversity carry over to the585 future? Phil. Trans. R. Soc. B, 374.1770, p. 20180111. DOI: 10.1098/rstb.2018.0111.586 Lande, R. & Schemske, D. W. (1985). The evolution of self-fertilization and inbreeding depression in plants. I. Genetic models. Evolution, 39.1,587 pp. 24–40. DOI: 10.1111/j.1558-5646.1985.tb04077.x.588 Lefkovitch, L. P. (1965). The study of population growth in organisms grouped by stages. Biometrics, 21.1, p. 1. DOI: 10.2307/2528348.589 Lewis, M. A., Shuai, Z., & van den Driessche, P. (2019). A general theory for target reproduction numbers with applications to ecology and590 epidemiology. J. Math. Biol., DOI: 10.1007/s00285-019-01345-4.591 Li, C.-K. & Schneider, H. (2002). Applications of Perron-Frobenius theory to population dynamics. J. Math. Biol. , 44.5, pp. 450–462. DOI: 10.592 1007/s002850100132.593 Lion, S. & Metz, J. A. J. J. (2018). Beyond R0 maximisation: on pathogen evolution and environmental dimensions. Trends Ecol. Evol., 33, pp. 75–594 90. DOI: 10.1016/j.tree.2018.02.004.595 Massol, F., Calcagno, V., & Massol, J. (2009). The metapopulation fitness criterion: Proof and perspectives. Theor. Popul. Biol., 75.2–3, pp. 183–596 200. DOI: http://dx.doi.org/10.1016/j.tpb.2009.02.005.597 Massol, F. & Débarre, F. (2015). Evolution of dispersal in spatially and temporally variable environments: The importance of life cycles. Evo-598 lution, 69.7, pp. 1925–1937. DOI: 10.1111/evo.12699.599 Maynard Smith, J. (1982). Evolution and the Theory of Games. Cambridge university press.600 18 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint Maynard Smith, J. & Price, G. R. (1973). The logic of animal conflict. Nature, 246, p. 15. DOI: 10.1038/246015a0.601 Metz, J. A. J. & Leimar, O. (2011). A simple fitness proxy for structured populations with continuous traits, with case studies on the evolution602 of haplo-diploids and genetic dimorphisms. J. Biol. Dynam., 5.2, pp. 163–190. DOI: 10.1080/17513758.2010.502256.603 Metz, J. A. & Gyllenberg, M. (2001). How should we define fitness in structured metapopulation models? Including an application to the604 calculation of evolutionarily stable dispersal strategies. Proc. R. Soc. B, 268.1466, pp. 499–508.605 Metz, J. A., Nisbet, R., & Geritz, S. (1992). How should we define ‘fitness’ for general ecological scenarios? Trends Ecol. Evol., 7.6, pp. 198–202.606 Mitchell, E., Graham, A. L., Úbeda, F., & Wild, G. (2022). On maternity and the stronger immune response in women. Nat. Commun., 13.1. DOI:607 10.1038/s41467-022-32569-6.608 Mullon, C., Keller, L., & Lehmann, L. (2016). Evolutionary stability of jointly evolving traits in subdivided populations. Am. Nat., 188.2, pp. 175–609 195. DOI: 10.1086/686900.610 Mullon, C. & Lehmann, L. (2014). The robustness of the weak selection approximation for the evolution of altruism against strong selection.611 J. Evol. Biol. , 27.10, pp. 2272–2282. DOI: 10.1111/jeb.12462.612 Pelletier, F., Clutton-Brock, T., Pemberton, J., Tuljapurkar, S., & Coulson, T. (2007). The Evolutionary Demography of Ecological Change: Linking613 Trait Variation and Population Growth. Science, 315.5818, pp. 1571–1574. DOI: 10.1126/science.1139024.614 Priklopil, T. & Lehmann, L. (2020). Invasion implies substitution in ecological communities with class-structured populations. Theor. Popul.615 Biol., 134, pp. 36–52. DOI: 10.1016/j.tpb.2020.04.004.616 — (2024). On the interpretation of the operation of natural selection in class-structured populations. Am. Nat., 203.2, pp. 292–304. DOI:617 10.1086/727970.618 Rees, M. & Ellner, S. P. (2016). Evolving integral projection models: evolutionary demography meets eco-evolutionary dynamics. Methods Ecol.619 Evol., 7.2, pp. 157–170. DOI: 10.1111/2041-210x.12487.620 Roberts, M. G. & Heesterbeek, J. A. P. (2003). A new method for estimating the effort required to control an infectious disease. Proc. R. Soc.621 B, 270.1522, pp. 1359–1364. DOI: 10.1098/rspb.2003.2339.622 Rodrigues, A. M. M. & Gardner, A. (2015). Simultaneous failure of two sex-allocation invariants: implications for sex-ratio variation within and623 between populations. Proc. R. Soc. B, 282.1810, p. 20150570. DOI: 10.1098/rspb.2015.0570.624 Rodrigues, A. M. M. & Gardner, A. (2012). Evolution of helping and harming in heterogeneous populations. Evolution, 66.7, pp. 2065–2079.625 DOI: 10.1111/j.1558-5646.2012.01594.x.626 — (2013). Evolution of helping and harming in heterogeneous groups. Evolution, 67.8, pp. 2284–2298. DOI: 10.1111/j.1558-5646.2012.627 01594.x.628 Rousset, F. (2004). Genetic Structure and Selection in Subdivided Populations (MPB-40). Princeton University Press.629 Rueffler, C. & Metz, J. A. J. (2013). Necessary and sufficient conditions for R0 to be a sum of contributions of fertility loops. J. Math. Biol., 66.4-5,630 pp. 1099–1122. DOI: 10.1007/s00285-012-0575-0.631 Rueffler, C., Metz, J. A. J., & Van Dooren, T. J. M. (2012). What life cycle graphs can tell about the evolution of life histories. J. Math. Biol., 66.1-2,632 pp. 225–279. DOI: 10.1007/s00285-012-0509-x.633 Stott, I., Townley, S., Carslake, D., & Hodgson, D. J. (2010). On reducibility and ergodicity of population projection matrix models. Methods Ecol.634 Evol., 1.3, pp. 242–252. DOI: 10.1111/j.2041-210x.2010.00032.x.635 Takada, T. & Kigami, J. (1991). The dynamical attainability of ESS in evolutionary games. J. Math. Biol. , 29.6, pp. 513–529. DOI: 10 . 1007 /636 bf00164049.637 Takada, T. & Nakajima, H. (1992). An analysis of life history evolution in terms of the density-dependent Lefkovitch matrix model. Math. Biosci.,638 112.1, pp. 155–176. DOI: 10.1016/0025-5564(92)90091-a.639 — (1998). Theorems on the invasion process in stage-structured populations with density-dependent dynamics. J. Math. Biol., 36.5, pp. 497–640 514. DOI: 10.1007/s002850050111.641 Taylor, P. D. (1988). An inclusive fitness model for dispersal of offspring. J. Theor. Biol., 130.3, pp. 363–378. DOI: 10.1016/s0022-5193(88)642 80035-3.643 — (1990). Allele-frequency change in a class-structured population. Am. Nat., 135, pp. 95–106. DOI: 10.1086/285034.644 — (1992). Altruism in viscous populations—an inclusive fitness model. Evol. Ecol., 6.4, pp. 352–356. DOI: 10.1007/bf02270971.645 Taylor, P. D. & Bulmer, M. G. (1980). Local mate competition and the sex ratio. J. Theor. Biol., 86.3, pp. 409–419. DOI: 10 . 1016 / 0022 -646 5193(80)90342-2.647 Tuljapurkar, S., Horvitz, C. C., & Pascarella, J. B. (2003). The many growth rates and elasticities of populations in eandom environments. Am.648 Nat., 162.4, pp. 489–502. DOI: 10.1086/378648.649 19 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint Úbeda, F. & Jansen, V. A. A. (2016). The evolution of sex-specific virulence in infectious diseases. Nat. Commun., 7, p. 13849. DOI: 10.1038/650 ncomms13849.651 van den Driessche, P. & Watmough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of652 disease transmission. Math. Biosci., 180.1, pp. 29–48. DOI: 10.1016/s0025-5564(02)00108-6.653 Wright, S. (1931). Evolution in Mendelian populations. Genetics, 16.2, p. 97. DOI: 10.1007/bf02459575.654 20 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint Box 1: Classifying NGMs based on the invasion threshold In the main text, we have seen that some NGMs have the same invasion condition. It is therefore useful to clas- sify typical formulas of the NGMs with shared invasion threshold. In this box, we say that two NGMs G1 and G2 are invadability-equivalent if ρ(G1)⋛ 1 ⇔ ρ(G2)⋛ 1. In Appendix C, we show that the following NGMs are invadability-equivalent, whenever F and S are both nonnegative regardless of their actual biological meanings: F + S, (I − S)−1 F, F (I − S)−1 and (F S F S )(Box-Fig 1). The last NGM has twice as large size as the original, and may be practically relevant when the class transition is unbiased (Rueffler et al. 2012; Lion & Metz 2018; Iritani et al. 2019, also see Example 2). We duplicated the focal class into two replicas with the same reproductive value. The life-cycle graph is expanded by creating a replica of the duplicated classes. Whenever the inversion exists and F and S are nonnegative, these NGMs are invadability-equivalent to each other. This result is therefore a simple generalization of the so-called ‘fundamental matrix’ in population demography (Cushing 1998; Li & Schneider 2002). One may wonder which NGM to use. One may preferably use the ‘snapshot’ fitness F + S for its numerical stability over the others involving matrix inversions. From an algebraic perspective, determining the stability properties of evolutionary equilibrium involves partial derivatives; again, the snapshot form ( F + S) may seem preferable, but at the cost of higher dimensionality than the inverse formulas such as F (I − S)−1 if F contains many zeros or redundancies (said to be ‘rank deficient’ in linear algebra). G = F + S ; (I − S)−1 F = F + SF + S2F + ⋯ ; F (I − S)−1 = F + FS + FS2 + ⋯ ; ( F S F S ); ≅ ≅ ≅ ≅ ### S F ### ### ⋮ F S S ⋮ ### ### ⋮ F S F S ### ### F S repl. repl. Box-Fig 1:The invadability-equivalent NGMs (top) and the corresponding structures of NGMs (bottom). The colors correspond to the highlighted terms in the equation. Closed circle: focal individual; open cicles: offspring of focal individual; Gray circle: focal-individual replica. Gray pathway should not be counted in to avoid the double-counting of fitness. “repl.” for mirroerd replica. 21 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint Box 2: Calculation techniques In practice, the invasion condition may involve fractional terms such as w ≔ ϕ(z + δ, z)/μ(z + δ, z), where z is resident phenotype, with δ small. The derivative has a useful formula: 1 w ⋅ 𝜕w 𝜕δ = (ϕ(1) ϕ − μ(1) μ ), where (1) represents the first derivative with respect to δ. This is a form of the marginal value theorem linking the trade- off relationship between the marginal increases in the numerator (typically, fecundity) and the denominator (typically, mortality) (Charnov 1976). In population genetics, this expression is called as the average excess of gene substitution (Fisher 1941). Invasion conditions often contain the cumulative survival probability terms like s′ ≔ ∏k s′ k where s′ k represents a survival or transition probability from class k and ∏k represents the product operator over k. Its derivative also has a useful formula: 1 s ⋅ 𝜕s 𝜕δ = ∑ k s(1) k sk . This expression can be also viewed as selection leading to the balanced, marginal increases in survival probabilities. The common technique behind these derivatives is logarithmic differentiation to compute marginal changes in re- productive success. The following, matrix-analogue of the logarithmic derivative of the sojourn time is also useful: 𝜕 𝜕δ (I − S)−1 = (I − S)−1 S(1) (I − S)−1 . (12) The matrix (I − S)−1, called as the fundamental matrix (Kemeny et al. 1966), represents the expected number of visits to a class starting from another class. A note on interpretability: while symbolic computation software is powerful, automated simplification often ob- scures biological insight by expanding meaningful compound terms (e.g., life-history trade-offs). Differentiation should be viewed not merely as an algebraic operation but as a mapping that preserves the biological structure of fitness. Therefore, we recommend rearranging the derived equations to maintain biological interpretability, and prioritizing meaningful groupings over mathematical brevity. 22 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint Box 3: Properties of the structural evolutionary invasion analysis Timescale separation First, we obtain the PNGM as a result of separation of timescales between primary (which we assume is slow) and secondary (fast) groups. Biologically, this means that we can focus on only an arbitrary subset of classes to capture the overall gene-frequency changes (Appendix C). Note that one does not need any a priori knowledge or assumption of the actual differences in timescales of the classes. The utility of this method is clearly illustrated in the dispersal model (Example 3), where the class is defined as the number of mutant in a patch and the timescale is separated based on the purpose of deriving fitness, regardless of biological assumption on the timescales. Fisher’s reproductive values Second, Fisher’s (1930) reproductive values of the primary group do not change with PNGM (Appendix C). This property reflects the reproductive values’ nature of measuring the asymptotically long- term contributions to the gene pool (Taylor 1990; Caswell 2001; Rousset 2004). Note that more precisely, the invariance in the reproductive values hold only for neutral mutants (first order effect; Taylor 1990) but not for selectively non-neutral mutants. Inheriting stability characteristics Practically, the ultimate goal of applying the adaptive dynamics framework is to lo- cate and characterize evolutionary equilibrium (EE), namely the direction of selection and possibility of disruptive selection. These characteristics of EE can be analyzed through: (i) the first derivative of invasion fitness (i.e., selection gradient capturing the location of EE); (ii) how the direction of selection changes with the resident trait value (i.e., convergence stability or attainability); and (iii) the second derivative of invasion fitness (i.e., whether selection is disruptive or stabilizing around EE; May- nard Smith & Price 1973; Maynard Smith 1982; Christiansen 1991; Takada & Kigami 1991; Eshel 1996). Regardless of using the original NGM, invasion determinant, or PNGM, the resulting invasion conditions give the same stability properties (i)–(iii). Thus, one can characterize both the location and stability of EE using the structural evolutionary invasion analysis. See Appendix D for details. 23 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint APPENDIX655 A Invasion determinant formula656 Resident NGM657 W e here show that the resident NGM has a unity spectral radius. This result is obvious for discrete-time models. W e thus derive658 the result for continuous-time models.659 T o see this, consider two types each with density # ‌N and # ‌N ′. Ignoring mutation, we can describe their dynamics as:660 d # ‌N dt = A(z | E )# ‌N , d # ‌N ′ dt = A(z′ | E )# ‌N ′, (A-1) where E = (# ‌N , # ‌N ′, z, z′). Note that these are generally nonlinear because of E.661 At the resident equilibrium with no mutant,662 A(z | E∘ )# ‌N ∘ = # ‌0 , E∘ = (# ‌N ∘, # ‌0 , z, z). (A-2) When the mutant is rare, its dynamics can be linearized as:663 d # ‌N ′ dt = A(z′ | E∘ )# ‌N ′ = J(z′ | E∘ )# ‌N ′. (A-3) Thus, the Jacobian matrix evaluated at the resident equilibrium is exactly J(z′ | E∘ )= A(z′ | E∘ ).664 At neutrality (z′ = z), substituting this into the resident equilibrium condition yields:665 J ∘ # ‌N ∘ = # ‌0 , (A-4) where J ∘ ≔ J(z | E∘ ).666 Decomposing the Jacobian as J ∘ = U ∘ − V ∘ according to the NGM framework, we get:667 (U ∘ − V ∘)# ‌N ∘ = # ‌0 ⟹ U ∘ # ‌N ∘ = V ∘ # ‌N ∘. (A-5) Since the neutral NGM is defined as G∘ = U ∘ (V ∘)−1, we get:668 G∘ (V ∘ # ‌N ∘)= V ∘ # ‌N ∘. (A-6) Because the matrix U ∘ is non-negative and the equilibrium density # ‌N ∘ is strictly positive, the vector U ∘ # ‌N ∘ must be non-669 negative and non-zero. The equation G∘U ∘ # ‌N ∘ = U ∘ # ‌N ∘ indicates that U ∘ # ‌N ∘ is an eigenvector of G∘ associated with the unit670 eigenvalue. By the Perron-Frobenius theorem for non-negative irreducible matrices, an eigenvalue associated with a non-671 negative eigenvector must be the spectral radius. Therefore, ρ(G∘)= 1.672 Theorem on the invasion determinant673 W e here prove the invasion determinant.674 Theorem 1 (Invasion determinant formula). Let G be a next-generation matrix, and assume that the selection is weak. Then,675 ρ(G)⋛ 1 ⟺ − det(I − G)⋛ 0. (A-7) 24 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint Proof. Let ΓM(x) ≔ det(xI − M)be the monic characteristic polynomial of a matrixM, and G∘ be the mutant NGMG evaluated676 at δ = 0 (i.e., the neutral mutant matrix). Thus, we have ρ(G∘) = 1 and so, by continuity of the spectrum of G, we have677 ρ(G) = ρ(G∘)+ ε = 1 + ε where ε is small. Furthermore, the Perron-Frobenius theorem for irreducible nonnegative matrix678 implies that the spectral radius of G is a simple eigenvalue of G and so ΓG(ρ(G))= 0 and d ΓG(x)/dx|x=ρ(G) > 0. As a result,679 if δ is small we can expand the equation ΓG(1+ ε) = 0 in powers of ε to obtain ΓG(1)+ dΓG(x)/ dx|x=ρ(G)ε + o(ε)= 0. Given680 ε = ρ(G) −1, for small enough δ and thus ε this can be re-written as sign(ρ(G) −1) = sign(−ΓG(1)). Finally, from the definition681 of ΓG(x)this can be re-written as sign(ρ(G) −1) = sign(−det(I− G))from which the theorem follows.682 The corresponding result for the original continuous-time Jacobian reads:683 ρ(G)⋛ 1 ⇔ − det(−J)⋛ 0. (A-8) T o prove this we decomposeJ = U − V into the components that satisfy the hypotheses of NGM. With some rearrangement,684 we get:685 det(−J)= det(V − U)= det(I − G)det(V). (A-9) W e notice thatV has nonnpositive elements at off-diagonal components while positive element at diagonal components, and686 V −1 is nonnegative. W e can thus, using the theory of M-matrix, we have det(V)> 0 (Berman & Plemmons 1987). Therefore,687 det(−J)and det(I − G)have the same sign, and thus both can be used for invasion condition.688 B Nonnegative decomposition theorem689 W e here show that any NGM can be split into nonnegative and nonzero matrices and converted to an alternative.690 Theorem 2(Nonnegative decomposition theorem). Let G be a (irreducible) next-generation matrix, and assume that the selec-691 tion is weak. Put G = H1 + H2 where H1 and H2 are both nonnegative and nonzero for values ofδ within some neighbourhood692 of δ = 0. Then,693 1. for small enough δ, the matrix I − H2 is invertible and thus ˆG ≔ H1 (I − H2)−1 is well-defined, and694 2. if ˆG is irreducible, the following equivalency holds:695 ρ(G)⋛ 1 ⟺ ρ(ˆG)⋛ 1. (B-1) Proof. W e first show that I − H2 is invertible whenever H1 is nonzero. For this, we use Corollary 2.1.5 of Berman & Plemmons696 (1987, p. 27), which is stated as follows:697 Lemma 1 (Berman and Plemmons 1994). Let A and B be nonnegative matrices. If B − A is nonnegative with B − A ≠ O, and698 if A + B is irreducible, then ρ(A)< ρ(B).699 W e takeA = H∘ 2 and B = G∘. W e haveB − A = H∘ 1 , which is nonnegative and nonzero by our assumption. Furthermore,700 because B = G∘ is irreducible by the standard assumption of the NGM, and O ≤ A ≤ B, the sum A + B shares the exact701 same zero-nonzero pattern as B, making A + B irreducible as well. Thus, the hypotheses of Lemma 1 are satisfied, yielding702 ρ(H∘ 2 )< ρ(G∘)= 1. By the continuity of the spectral radius with respect to matrix entries, ρ(H2)is continuous in δ. Therefore,703 ρ(H2)< 1 under weak selection (i.e., for small enough |δ|), which implies that I − H2 is invertible and completes the proof of704 the first statement.705 T o prove the second part, we use the invasion determinant. By definition of ˆG = H1 (I − H2)−1, we have:706 − det(I − ˆG)= − det(I − H1 − H2) det(I − H2) = − det(I − G) det(I − H2). (B-2) 25 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint Because ρ(H2)707 0. Therefore, − det(I − ˆG)and − det(I − G)share the exact same sign.708 Because G is irreducible, Theorem 1 guarantees thatρ(G)⋛ 1 ⟺ − det(I − G)⋛ 0. Likewise, because the reduced matrix709 ˆG is explicitly assumed to be irreducible in the theorem statement, applying Theorem 1 to ˆG yields − det(I − ˆG) ⋛ 0 ⟺710 ρ(ˆG)⋛ 1. Connecting these equivalences through the shared sign of their determinants completes the proof.711 The nonnegative decomposition theorem is based on a matrix splitting method known as ‘regular splitting’ for numerical712 analysis (V arga 1963).713 C PNGM714 Derivation of PNGM715 W e here prove the following theorem.716 Theorem 3 (Projected Next-Generation Theorem). Let G be a next-generation matrix of a mutant, with the following block-717 partition structure:718 G = ( G𝕏,𝕏 G𝕏,𝕐 G𝕐,𝕏 G𝕐,𝕐 ). (C-1) Then, under weak selection, the following results hold true:719 1. ρ(G𝕐,𝕐)< 1.720 2. ρ(G)⋛ 1 ⟺ ρ(𝒫𝕏(G))⋛ 1, (C-2) where721 𝒫𝕏(G)≔ G𝕏,𝕏 + G𝕏,𝕐 (I − G𝕐,𝕐) −1 G𝕐,𝕏. (C-3) The first part can be proved in the same way as the first part of Theorem 1. The second part is also a direct implication of722 Theorem 2: choose the decomposition as H1 ≔ Π𝕏G with H2 = G − H1 = Π𝕐G, where Π𝕏 and Π𝕐 are projection matrices723 (Roberts & Heesterbeek 2003); note that a matrix Π is called as projection if Π2 = Π.724 W e here provide a similar, alternative proof. Specifically, we use the Schur complement and Schur’ s determinant identity725 (Crabtree & Haynsworth 1969).726 Definition 1 (Schur complement). Let M be a square matrix block-partitioned as:727 M = ( M𝕏,𝕏 M𝕏,𝕐 M𝕐,𝕏 M𝕐,𝕐 ) (C-4) Suppose M𝕐,𝕐 is invertible. W e define the Schur complement ofM to M𝕐,𝕐, denoted by 𝒮𝕏(M), as728 𝒮𝕏(M)≔ M𝕏,𝕏 − M𝕏,𝕐M−1 𝕐,𝕐M𝕐,𝕏. (C-5) By definition, 𝒫𝕏(M)= I − 𝒮𝕏(I − M).729 Lemma 2 (Schur). Let M be a square matrix satisfying the assumptions of Definition 1. Then, the following identity holds true:730 det(𝒮𝕏(M))= det(M) det(M𝕐,𝕐). (C-6) 26 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint (Proof of Lemma 2). Consider the identity:731 (M𝕏,𝕏 M𝕏,𝕐 M𝕐,𝕏 M𝕐,𝕐 ) ( I O −M−1 𝕐,𝕐M𝕐,𝕏 I ) = (𝒮𝕏(M) M𝕏,𝕐 O M 𝕐,𝕐 ), (C-7) taking the determinant of which yields:732 det(M)⋅ 1 = det(𝒮𝕏(M))⋅ det(M𝕐,𝕐), (C-8) which yields the desired result.733 (Proof of Theorem 3). In Lemma 2, put M = I − G. By definition, we have:734 𝒮𝕏(M)= I − 𝒫𝕏(G)= I − G𝕏,𝕏 − G𝕏,𝕐 (I − G𝕐,𝕐) −1 G𝕐,𝕏. (C-9) Applying Lemma 2, we have:735 det(I − 𝒫𝕏(G))= det(I − G) det(I − G𝕐,𝕐). (C-10) Because ρ(G𝕐,𝕐) 0. Hence the736 LHS and the numerator of the RHS have the same sign. This completes the proof.737 Interpretation of PNGM738 W e consider a mutant individual residing in one of the primary group classes and assess how many offspring it produces to the739 primary group. Specifically, we write Pn for a matrix of the total expected number of the mutant’s offspring produced to the740 primary group over n subsequent generations. W e then get:741 Pn = G𝕏,𝕏 + G𝕏,𝕐G𝕐,𝕏 + G𝕏,𝕐G𝕐,𝕐G𝕐,𝕏 + ⋯ + G𝕏,𝕐(G𝕐,𝕐)n−2G𝕐,𝕏 = G𝕏,𝕏 + n−2 ∑ i=0 G𝕏,𝕐(G𝕐,𝕐)iG𝕐,𝕏 = G𝕏,𝕏 + G𝕏,𝕐 (I − (G𝕐,𝕐)n−1) (I − G𝕐,𝕐) −1 G𝕐,𝕏 (C-11) where the first term represents the number of primary group-offspring directly produced by the mutant in one generation; the742 second term represents the number of primary group-grand-offspring whose parents are in the secondary group (the remaining743 terms can be read analogously). Because ρ(G𝕐,𝕐)< 1, we have lim n→∞ Pn = 𝒫𝕏(G)(Horn & Johnson 1985, Theorem 5.6.12).744 Hence, the PNGM counts the total number of offspring in the primary classes produced by the mutant from the primary classes745 throughout the lifespan.746 Embedding of PNGM747 Here we examine the relation between PNGM and nonnegative decomposition theorem. W e specifically consider embedding748 G = H1 + H2 into a high dimensional space as:749 ˜G ≔ ( O I H1 H2 ), (C-12) with:750 𝒫𝕏(˜G)= (I − H2)−1 H1 𝒫𝕐(˜G)= G = H1 + H2 (C-13) 27 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint X1 X2 ⋮ XK Y1 Y2 ⋮ YK ≅ ˜G = ( O I H1 H2 ) ⟼ { 𝒫𝕏(˜G) = (I − H2)−1 H1 𝒫𝕐(˜G) = H1 + H2 Figure 4: Embedding of G and its projection. Here, Y s are the orignal classes and X are their ‘replicas’. (Figure 4). On the other hand, we also have751 (I − H2)−1 H1 = H1 + H2 (I − H2)−1 H1 = 𝒫𝕏 ( H1 H2 H1 H2 ). (C-14) Since H1 (I − H2)−1 and H2 (I − H1)−1 have the same spectral radius, the following matrices share the same invasion conditions:752 H1 + H2, H1 (I − H2)−1 , (I − H2)−1 H1 ( H1 H2 H1 H2 ) . (C-15) The relation to separation of timescales753 We here show that the separation of timescales among classes produces the PNGM.754 Theorem 4 (Constructing PNGM from continuous-time mutant dynamics). Consider the following, K-dimensional mutant755 dynamics:756 d # ‌N ′ dt = J # ‌N ′, (C-16) where # ‌N ′ represents a vector collecting the densities in the classes, and J represents the Jacobian, block-partitioned as Equa-757 tion (C-1). Consider the Schur-complement 𝒮𝕏(J)= J𝕏,𝕏 − J𝕏,𝕐J −1 𝕐,𝕐J𝕐,𝕏. Then,758 (1) taking the quasi-equilibrium approximation d # ‌N ′ 𝕐/ dt ≈ #‌0 leads to d # ‌N ′ 𝕏/ dt = 𝒮𝕏(J)# ‌N𝕏;759 (2) there exists a decomposition of J = U − V that satisfies the hypotheses of NGM theorem, V𝕏,𝕏 nonsingular, and 𝒮𝕏(J)=760 𝒫𝕏(UV −1)V𝕏,𝕏 − V𝕏,𝕏. Thus for such U, V with UV −1 = G, the resulting NGM is 𝒫𝕏(G)V𝕏,𝕏 (V𝕏,𝕏) −1 = 𝒫𝕏(G).761 Proof. (1) This statement results from calculating the Schur-complement for the original K-dimensional Jacobian: consider762 the quasi-equilibrium approximation given by:763 d dt # ‌N ′ 𝕐 = J𝕐,𝕏 # ‌N ′ 𝕏 + J𝕐,𝕐 # ‌N ′ 𝕐 ≈ #‌0 , (C-17) which gives # ‌N ′ 𝕐 ≈ −J −1 𝕐,𝕐J𝕐,𝕏 # ‌N ′ 𝕏, and thus764 d # ‌N ′ 𝕏 dt ≈ J𝕏,𝕏 # ‌N ′ 𝕏 + J𝕏,𝕐 # ‌N ′ 𝕐 ≈ (J𝕏,𝕏 − J𝕏,𝕐J −1 𝕐,𝕐J𝕐,𝕏)# ‌N ′ 𝕏 = 𝒮𝕏(J)# ‌N ′ 𝕏. (C-18) (2) It is easy to check that if we choose a block-diagonal, invertible, and inverse-nonnegative matrix as V, the theorem imme-765 diately follows.766 28 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint This theorem guarantees the equivalence of quasi-equilibrium approximation (sepataion of timescales) and PNGM: the767 primary group dynamics are taken as ‘slow’ while the secondary group dynamics are ‘fast. ’768 W e can consequently derive the discrete-time analogue of the separation of timescales.769 Proposition 1 (Constructing PNGM from discrete-time mutant dynamics). Suppose the mutant dynamics:770 # ‌N ′(t + 1)= G # ‌N ′(t), (C-19) where G represents a next-generation matrix of the system, block-partitioned as Equation (C-1). Then, the PNGM obtains by771 taking the quasi-equilibrium approximation of fast-primary group and slow-secondary group.772 Proof. W e explicitly solve G𝕐,𝕏 # ‌N𝕏 + G𝕐,𝕐 # ‌N𝕐 ≈ # ‌N𝕐 to get # ‌N𝕐 ≈ (I − G𝕐,𝕐) −1 G𝕐,𝕏 # ‌N𝕏. Substituting this into the dynamics of773 # ‌N𝕏 yields774 # ‌N𝕏(t + 1)= G𝕏,𝕏 # ‌N𝕏 + G𝕏,𝕐 # ‌N𝕐 ≈ (G𝕏,𝕏 + G𝕏,𝕐 (I − G𝕐,𝕐) −1 G𝕐,𝕏) # ‌N𝕏 = 𝒫𝕏(G)# ‌N𝕏(t), (C-20) which proves the proposition.775 Reproductive values776 Here we show that the eigenvectors for the primary group at neutrality are unchanged with the projection. Let # ‌r denote the777 right eigenvector of G∘ (at neutrality) associated with ρ(G∘)= 1:778 G∘ # ‌r = # ‌r . (C-21) Now suppose that the vector is partitioned so the dimension is consistent with the block-partitioning of G∘; namely:779 ( #‌r𝕏 #‌r𝕐 ) = G∘ ( #‌r𝕏 #‌r𝕐 ) = (G∘ 𝕏,𝕏 G∘ 𝕏,𝕐 G∘ 𝕐,𝕏 G∘ 𝕐,𝕐 ) ( #‌r𝕏 #‌r𝕐 ) = (G∘ 𝕏,𝕏 #‌r𝕏 + G∘ 𝕏,𝕐 #‌r𝕐 G∘ 𝕐,𝕏 #‌r𝕏 + G∘ 𝕐,𝕐 #‌r𝕐 ). (C-22) By solving for #‌r𝕐 the linear equation at the second block, we have:780 #‌r𝕐 = (I − G∘ 𝕐,𝕐) −1 G∘ 𝕐,𝕏 #‌r𝕏, (C-23) substituting which into the first block gives:781 #‌r𝕏 = (G∘ 𝕏,𝕏 + G∘ 𝕏,𝕐 (I − G∘ 𝕐,𝕐) −1 G∘ 𝕐,𝕏) #‌r𝕏 = 𝒫𝕏(G∘)#‌r𝕏, (C-24) as desired. The same argument applies to the left eigenvector associated with the spectral radius ρ(G∘)= 1.782 D Stability analysis783 W e have presented various invadability-equivalent NGMs withκ1(δ, z)≔ ρ(G1)− 1 ⋛ 0 if and only if κ2(δ, z)≔ ρ(G2)− 1 ⋛ 0.784 From this expression, we can derive three properties of long-term evolutionary dynamics, namely (i) the location of evolutionary785 equilibrium (EE), (ii) convergence stability, and (iii) evolutionary stability.786 W e write sign()for the sign function.787 Location of evolutionary equilibrium and convergence stability788 W e expandκi(δ, z)to the first order in δ as:789 κi(δ, z)= δκ(1) i (z)+ O(δ2), (D-1) 29 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint where κ(1) i (z)≔ 𝜕κi 𝜕δ || |δ=0 is the selection gradient, and the zero-th order term vanishes due to neutrality (κ i(0, z)= 0).790 For any given resident trait z, κ1(δ, z)and κ2(δ, z)share the same sign for sufficiently small |δ| > 0. For this point-wise791 property to hold, their leading-order coefficients in the Taylor expansion evaluated at z must also share the same sign. Thus,792 we have:793 sign(κ(1) 1 (z))= sign(κ(1) 2 (z)). (D-2) From this equivalence, two properties immediately follow: (i) κ(1) 1 (z∗)= 0 if and only if κ(1) 2 (z∗)= 0, which proves that the794 location of the evolutionary equilibrium (EE) is identical. (ii) In the neighborhood of the EE, the sign of the selection gradient795 is identical for both measures. This implies that the direction of trait substitution is exactly the same, thereby ensuring that796 their convergence stability agrees.797 Evolutionary stability798 Fix the resident trait at an evolutionary equilibrium, z = z∗. W e expand the fitness functions with respect toδ up to the second799 order:800 κi(δ, z∗)= κi(0, z∗)+ δκ(1) i (z∗)+ 1 2 δ2κ(2) i (z∗)+ O(δ3). (D-3) The zero-th order term vanishes due to neutrality (κi(0, z∗)= 0), and the first-order term vanishes becausez∗ is an evolutionary801 equilibrium (κ(1) i (z∗)= 0).802 Because κ1(δ, z∗)and κ2(δ, z∗)share the same sign for any sufficiently small |δ| > 0, their leading non-zero terms in the803 expansion must also share the same sign. Thus, we conclude that804 sign(κ(2) 1 (z∗))= sign(κ(2) 2 (z∗)), (D-4) which ensures that the condition for evolutionary stability is identical between any two invadability-equivalent NGMs.805 Reproductive value-weighted selection gradient806 W e show that the selection gradients derived either by the original NGM and PNGM agree. Let us derive the selection gradient807 using the reproductive values: we write # ‌ℓ ⊤ or # ‌r for the left or right eigenvector of neutral NGMG∘ associated with the spectral808 radius 1, respectively (Taylor & Frank 1996; Frank 1998; Avila & Mullon 2023). Given PNGM, we partition the eigenvectors809 according to the same partition-size:810 # ‌ℓ = ( # ‌ℓ𝕏 # ‌ℓ𝕐 ), # ‌r = ( #‌r𝕏 #‌r𝕐 ). (D-5) The first derivative of the PNGM, with the following lines being evaluated at neutrality z′ = z unless otherwise stated, is811 given by:812 𝜕𝒫(G(z′, z)) 𝜕δ = 𝜕 𝜕δ (G𝕏,𝕏 + G𝕏,𝕐 (I − G𝕐,𝕐) −1 G𝕐,𝕏) = G(1) 𝕏,𝕏 + G(1) 𝕏,𝕐 (I − G∘ 𝕐,𝕐) −1 G∘ 𝕐,𝕏 + G∘ 𝕏,𝕐 (I − G∘ 𝕐,𝕐) −1 G(1) 𝕐,𝕐 (I − G∘ 𝕐,𝕐) −1 G∘ 𝕐,𝕏 + G∘ 𝕏,𝕐 (I − G∘ 𝕐,𝕐) −1 G(1) 𝕐,𝕏. (D-6) 30 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint Premultiplying the vector # ‌ ℓ⊤ 𝕏 and also postmultiplying #‌r𝕏 gives:813 # ‌ ℓ⊤ 𝕏 𝜕𝒫(G(z′, z)) 𝜕δ #‌r𝕏 = ( # ‌ ℓ⊤ 𝕏, # ‌ ℓ⊤ 𝕐)𝜕G(z′, z) 𝜕δ ( #‌r𝕏 #‌r𝕐 ) = 𝜕ρ(G) 𝜕δ , (D-7) because:814 (I − G∘ 𝕐,𝕐) −1 G∘ 𝕐,𝕏 #‌r𝕏 = #‌r𝕐, (D-8) # ‌ ℓ⊤ 𝕏G∘ 𝕏,𝕐 (I − G∘ 𝕐,𝕐) −1 = # ‌ ℓ⊤ 𝕐. (D-9) W e have therefore confirmed the equivalence between the convergence stability conditions for the original and projected next-815 generation matrix.816 E Analyses of the examples817 Example 1: Two-class model818 W e describe the result of NGM by taking the decomposition of J to be:819 U = (b′ 1 − a′ 1 b′ 2 a′ 1 0 ), V = (m1 0 0 m2 ). (E-1) W e writef1 ≔ b′ 1/(m1 + a′ 1),s ≔ a′ 1/(m1 + a′ 1)and f2 ≔ b′ 2/(m2)for the fecundity of class-1 individual, survival probability820 (0 ≤ s ≤ 1), and fecundity of class-2 individual, respectively. A direct calculation then gives:821 G = UV −1 = (b′ 1 − a′ 1 b′ 2 a′ 1 0 ) (m1 −1 0 0 m2 −1) = ( f1 − s 1 − s f2 s 1 − s 0 ) . (E-2) Using the invasion determinant yields the invasion condition, of:822 − det(I − G)= − 1 − f1 − sf2 1 − s > 0. (E-3) W e thus get the invasion condition in the main text. Note that we have not assumed anything about the sign off1 − s, and thus823 G need not be nonnegative.824 Example 2: Epidemiological model825 The mutant dynamics is given by:826 dS′ f dt = θffb′ f ⋅ S∘ m + I∘ m N ∘ T ⋅ (S′ f + (1− v)I′ f ) +θfmb′ m ⋅ S′ m + I′ m N ∘ T ⋅ (S∘ f + (1− v)I∘ f ) +γ′ fI′ f − hS′ f − μS′ f dI′ f dt = θffb′ f ⋅ S∘ m + I∘ m N ∘ T ⋅ vI′ f + θfmb′ m ⋅ S′ m + I′ m N ∘ T ⋅ vI∘ f + hS′ f − (μ + γ′ f + αf)I′ f dS′ m dt = θmfb′ f ⋅ S∘ m + I∘ m N ∘ T ⋅ (S′ f + (1− v)I′ f ) +θmmb′ m ⋅ S′ m + I′ m NT ⋅ (S∘ f + (1− v)I∘ f ) +γ′ mI′ m − hS′ m − μS′ m dI′ m dt = θmfb′ f ⋅ S∘ m + I∘ m N ∘ T ⋅ vI′ f + θmmb′ m ⋅ S′ m + I′ m N ∘ T ⋅ vI∘ f + hS′ m − (μ + γ′ m + αm)I′ m. (E-4) 31 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint The resulting Jacobian reads:827 J = ( Uff Ufm Umf Umm ) − ( Vf O O V mm ), G = ( Uff Ufm Umf Umm ) ( V −1 f O O V −1 mm ), (E-5) where, writing x∘ m ≔ (S∘ m + I∘ m)/(N∘ T)for the sex ratio, y∘ f ≔ If/NT for the proportion of infected females, and θk,m for the828 probability that an individual of sex k derives a gene from a parent of sex m (for diploid organisms, these are all 1/2):829 Uff = θff (η′ fx∘ m η′ f(1− v)x∘ m 0 η′ fvx∘ m ), Ufm = θfm (, η′ m(1− x∘ m − vy∘ f ) η′ m(1− x∘ m − vy∘ f ) η′ mvy∘ f η′ mvy∘ f ), Umf = θmf (η′ fx∘ m η′ f(1− v)x∘ m 0 η′ fvx∘ m ), Umm = θmm (η′ m(1− x∘ m − vy∘ f ) η′ m(1− x∘ m − vy∘ f ) η′ mvy∘ f η′ mvy∘ f ), Vf = (h + μ −γ′ f −h μ + γ′ f + αf ), Vm = (h + μ −γ′ m −h μ + γ′ m + αm ). (E-6) Although we can write down the Jacobian (which is four-by-four) as above, we here analyze the model more formally without830 doing so. This approach would help us more easily code numerical analyses. For this aim, we use tensor decomposition. Note831 that the idea underlying the tensor decomposition is symmetry; for example, the vertical transmission (from mother) to her832 offspring is independent of offspring sex (at random with no bias). This enables us to handle the four-by-four matrix in more833 compact, two-by-two form:834 J = ( (θff θfm θmf θmm )⏟ ⎵ ⏟ ⎵ ⏟ =Θ ⊗ (1 0 0 1 )) ⎛ ⎜⎜ ⎝ η′ f (1 0 0 0 )⏟ =Πf ⊗ (x∘m (1−v)x∘m 0 vx∘m )⏟ ⎵ ⎵ ⏟ ⎵ ⎵ ⏟ =X∘m +η′ m (0 0 0 1 )⏟ =Πm ⊗ ( 1−x∘m−vy∘ f 1−x∘m−vy∘ f vy∘ f vy∘ f )⏟ ⎵ ⎵ ⎵ ⎵ ⎵ ⏟ ⎵ ⎵ ⎵ ⎵ ⎵ ⏟ =X∘ f ⎞ ⎟⎟ ⎠ − (Vf O O V m ), (E-7) which, using the mixed-product rule (Graham 2018), we can rewrite as:835 J = (η′ fΘΠf)⊗ X∘ m + (η′ mΘΠm)⊗ X∘ f − (Vf O O V m ), (E-8) where Πs are projections:836 Πf = ( I O O O ), Πm = ( O O O I ). (E-9) W e invert the lastV-matrix and get NGM:837 G = ((η′ fΘΠf)⊗ (X∘ m)) (Vf O O V m ) −1 + ((η′ mΘΠm)⊗ (X∘ f )) (Vf O O V m ) −1 . (E-10) Using Πf and Πm, we can rewrite the inversion as:838 (Vf O O V m ) −1 = (V −1 f O O V −1m )= Πf ⊗ V −1 f + Πm ⊗ V −1 m . (E-11) 32 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint Again using the mixed product rule and ΠfΠm = O, Π2 f = Πf and Π2 m = Πm, we get:839 G = (η′ fΘΠf)⊗ (X∘ mV −1 f )+ (η′ mΘΠm)⊗ (X∘ f V −1 m ). (E-12) W e note that840 V −1 f = 1 1 − ιv′ f ( 1 h+μ v′ f ι 1 μ+γ′ f +αf ) , (E-13) V −1 m = 1 1 − ιv′m ( 1 h+μ v′ m ι 1 μ+γ′m+αm ) (E-14) with ι = h/(h+ μ),v′ f = γ′ f/(μ+ γ′ f + αf)and v′ m = γ′ m/(μ+ γ′ m + αm). Taken together,841 G = ( θffGf θfmGm θmfGf θmmGm ). (E-15) For diploidy,842 G = 1 2 ( Gf Gm Gf Gm ). (E-16) Since this takes the form of block-partitioned matrix in Box 1, the invasion condition reduces to:843 ρ(Gf + Gm 2 ) > 1. (E-17) For haplodiploidy,844 G = ( 1 2 Gf 1 2 Gm Gf O ) . (E-18) By applying the PNGM to eliminate the male class, we get:845 𝒫f(G)= 1 2 Gf + 1 2 GmGf. (E-19) The first component represents the direct transmission of genes from female to female, and the second component represents846 the indirect transmission from female to female via male. The invasion condition for haplodiploidy is thus847 ρ(1 2 Gf + 1 2 GmGf) > 1. (E-20) Example 3: Lefkovitch model848 W e consider the following NGM:849 G = ⎛ ⎜ ⎜ ⎜ ⎝ g1,1 g1,2 g1,3 g1,4 g2,1 g2,2 0 0 0 g3,2 g3,3 0 0 0 g4,3 g4,4 ⎞ ⎟ ⎟ ⎟ ⎠ . (E-21) 33 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint This NGM does not have a closed form of spectral radius. W e partition the classes into primary and secondary groups: 𝕏 =850 {1, 2, 3} and 𝕐 = {4}. The PNGM then reads:851 𝒫𝕏(G)= ⎛ ⎜⎜ ⎝ g1,1 g1,2 g1,3 g2,1 g2,2 0 0 g3,2 g3,3 ⎞ ⎟⎟ ⎠ + ⎛ ⎜⎜ ⎝ g1,4 0 0 ⎞ ⎟⎟ ⎠ 1 1 − g4,4 (0 0 g4,3). (E-22) The second matrix is zero except for the (1,3)-rd elementg1,4(1− g4,4)−1g4,3, which recovers the result of the main text.852 T o derive the ‘canonical form’ invasion condition∑k fk ∏m sm > 1, we first partition the NGM into fecundity and transition853 components G = F + S, with:854 F = ⎛ ⎜ ⎜ ⎜ ⎝ f1,1 f1,2 f1,3 f1,4 0 0 0 0 0 0 0 0 0 0 0 0 ⎞ ⎟ ⎟ ⎟ ⎠ , S = ⎛ ⎜ ⎜ ⎜ ⎝ s1,1 0 0 0 s2,1 s2,2 0 0 0 s3,2 s3,3 0 0 0 s4,3 s4,4 ⎞ ⎟ ⎟ ⎟ ⎠ . (E-23) Inverting I − S yields:855 (I − S)−1 m,k = 1 1 − sk,k m ∏ i=k+1 si,i−1 1 − si,i , m ≥ k. (E-24) Since F has elements only in the first row, what we actually need is856 (I − S)−1 m,1 = 1 1 − s1,1 m ∏ i=2 si,i−1 1 − si,i . (E-25) Using this we get:857 G = F (I − S)−1 = ⎛ ⎜ ⎜ ⎜ ⎝ g1,1 ∗ ∗ ∗ 0 0 0 0 0 0 0 0 0 0 0 0 ⎞ ⎟ ⎟ ⎟ ⎠ , g1,1 = K ∑ m=1 f1,m 1 − sm,m m−1 ∏ i=1 si+1,i 1 − si,i , (E-26) which thus says that ρ(G)= g1,1, leading to the canonical form given in the main text following the transformation:858 fm ≔ f1,m 1 − sm,m , si ≔ si+1,i 1 − si,i , (E-27) by abuse of notation.859 W e can derive the same canonical form by recursively applying𝒫{1,2,3}, 𝒫{1,2}, 𝒫{1}:860 G = ⎛ ⎜ ⎜ ⎜ ⎝ f1,1 + s1,1 f1,2 f1,3 f1,4 s2,1 s2,2 0 0 0 s3,2 s3,3 0 0 0 s4,3 s4,4 ⎞ ⎟ ⎟ ⎟ ⎠ 𝒫{1,2,3} ⟼ ⎛ ⎜⎜ ⎝ f1,1 + s1,1 f1,2 f1,3 s2,1 s2,2 0 0 s3,2 s3,3 ⎞ ⎟⎟ ⎠ + 1 1 − s4,4 ⎛ ⎜⎜ ⎝ 0 0 f1,4s4,3 0 0 0 0 0 0 ⎞ ⎟⎟ ⎠ 𝒫{1,2} ⟼ … 𝒫{1} ⟼ f1,1 + s1,1 + K ∑ m=2 f1,m 1 − sm,m m−1 ∏ i=1 si+1,i 1 − si,i . (E-28) Rearranging the last quantity> 1 yields the canonical form of invasion condition, with the abuse of notationfm ≔ f1,m/(1−sm,m)861 34 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint and si ≔ si+1,i/(1− si,i).862 Example 4: Dispersal in stable habitats863 Moran process864 The transition probabilities Sk,m are given by Equation (25) (e.g., Mullon & Lehmann 2014). Based on these probabilities, we865 can describe the dynamics of the number of mutants per patch as:866 q′ 1(t + 1)− q′ 1(t)= ϕ′ 0 K ∑ k=1 kq′ k(t)+ μ′ 2q′ 2(t)− (ϕ′ 1 + μ′ 1)q′ 1(t), ⋮ q′ m(t + 1)− q′ m(t)= μ′ m+1q′ m+1(t)+ ϕ′ m−1q′ m−1(t)− (ϕ′ m + μ′ m)q′ m(t), ⋮ q′ K(t + 1)− q′ K(t)= ϕ′ K−1q′ K−1(t)− μ′ Kq′ K(t). (E-29) W e now perform the separation of timescales. W e eliminate the state variables in descending order from q′ K, q′ K−1, to q′ 2.867 First, setting the left-hand side of the last line to zero yields q′ K ≈ ϕ′ K−1q′ K−1/μ′ K. Second, we consider868 q′ K(t + 1)+ q′ K−1(t + 1)− (q′ K(t)+ q′ K−1(t))= ϕ′ K−2q′ K−2 − μ′ K−1q′ K−1 (E-30) and set it to zero to obtain869 q′ K−1 ≈ ϕ′ K−2 μ′ K−1 q′ K−2. (E-31) By induction, we have870 q′ m ≈ ϕ′ m−1 μ′m q′ m−1, (E-32) and thus871 q′ m ≈ q′ 1 m−1 ∏ k=1 ϕ′ k μ′ k+1 , m = 1, … ,K. (E-33) Substituting these expressions into the dynamics of class 1 leads to:872 q′ 1(t + 1)− q′ 1(t)≈ ϕ′ 0 ( K ∑ m=1 m ( m−1 ∏ k=1 ϕ′ k μ′ k+1 ))q′ 1(t)− μ′ 1q′ 1(t). (E-34) The invasion condition is determined by requiring the right-hand side to be positive, which gives873 K ∑ m=1 m m ∏ k=1 ϕ′ k−1 μ′ k > 1, (E-35) as desired.874 The W right-Fisher model and the advantage of the structural approach875 The standard approach to evolutionary invasion analysis in subdivided populations, particularly under the Wright-Fisher pro-876 cess with local density dependence, heavily relies on tracking the probability distribution of local patch states (i.e., the number877 of mutants in a patch). As brilliantly demonstrated by Mullon et al. (2016, 2018), one can evaluate the invasion condition by878 differentiating the lineage fitness. However, because the Wright-Fisher process yields a dense transition matrix, analytically879 solving the recursions for the perturbed mutant distributions or identical-by-descent (IBD) probabilities requires remarkable880 35 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint mathematical effort and biological intuition. This complexity severely limits the generalizability of the traditional approach881 when the model incorporates additional class structures such as sex, age, or environmental heterogeneity.882 Here, we demonstrate how our structural framework, specifically the synthesis of the projected next-generation matrix883 (PNGM) and the invasion determinant, effortlessly bypasses this high-dimensional hurdle.884 W e start by constructing the fullK × K next-generation matrix G for a focal patch, where the state k ∈ {1, 2, … ,K} repre-885 sents the number of mutants in the patch. W e biologically decompose the NGM into two distinct pathways: local philopatric886 transitions S and global seeding of new patches F, such that887 G = F + S. (E-36) The matrix S describes the local demography (e.g., binomial sampling in the Wright-Fisher process) where mutants stay in the888 focal patch. Importantly, ρ(S)< 1 because without global seeding, the local mutant lineage is doomed to ultimate extinction889 under the infinite-island assumption.890 The matrix F describes the production of successful emigrants that colonize resident patches. In the infinite-island model,891 a dispersing mutant will almost surely land in a patch entirely occupied by residents, creating a new local lineage starting with892 exactly one mutant. This biological feature implies a rank-one structure for F:893 F = # ‌e1 # ‌ f ⊤, (E-37) where # ‌e1 = (1, 0, … ,0)⊤ is the unit vector representing the new state of exactly one mutant, and# ‌f = (f1, f2, … ,fK)⊤ is the vector894 whose k-th element represents the expected number of successful global seedings produced by a patch currently containing k895 mutants.896 Under the traditional method, evaluating the spectral radius ρ(G)> 1 directly is analytically prohibitive for a dense matrix.897 However, using the invasion determinant, the invasion condition is equivalent to − det(I − G)> 0. Substituting the rank-one898 structure into the determinant and applying the matrix determinant lemma (the Sherman-Morrison formula), we obtain:899 det(I − G)= det(I − S − # ‌e1 # ‌ f ⊤)= det(I − S)det(1 − # ‌ f ⊤ (I − S)−1 # ‌e1). (E-38) Since ρ(S) 0. Thus, the invasion condition collapses exactly to a single scalar inequality:900 W ≔ # ‌ f ⊤ (I − S)−1 # ‌e1 > 1. (E-39) This scalar W is precisely the 1 × 1 PNGM projected onto the state of a newly seeded patch (i.e., primary group 𝕏 = {1}), the901 so-called metapopulation fitness (Metz & Gyllenberg 2001; Ajar 2003; Massol et al. 2009; Mullon et al. 2016). W e have entirely902 eliminated the need to explicitly track the K × K mutant frequency dynamics. The matrix (I − S)−1 is the fundamental matrix,903 where the k-th entry of # ‌τ ≔ (I − T)−1 # ‌e1 simply represents the expected sojourn time in state k starting from a single founder.904 The true power of this structural approach shines under weak selection (first-order effects). Instead of solving complex905 recursions for the perturbed mutant distributions 𝜕qk/𝜕δ as required in previous studies (Mullon et al. 2016, 2018), we can906 mechanically Taylor-expand the scalar fitness W. Let δ be the small phenotypic deviation of the mutant. W e expand S ≈907 S∘ + δS(1) and # ‌f ≈ # ‌f ∘ + δ # ‌f (1). Using the standard matrix identity for the derivative of an inverse, (I − S)−1 ≈ (I − S∘)−1 +908 δ (I − S∘)−1 S(1) (I − S∘)−1, the selection gradient emerges automatically:909 𝜕W 𝜕δ || |δ=0 = # ‌ f ⊤(1) # ‌τ∘ + # ‌ ℓ⊤S(1) # ‌τ∘, (E-40) where # ‌τ∘ = (I − S∘)−1 # ‌e1 is the neutral sojourn time vector, and # ‌ ℓ⊤ = # ‌f ∘⊤ (I − S∘)−1 is proportional to (1,2, … ,K)⊤.910 In this formulation, the term # ‌ ℓ⊤S(1)∘ elegantly and automatically generates the inclusive fitness effects (kin selection). Under911 36 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint neutrality, all individuals are selectively equivalent and thus have the exact same individual reproductive value. Consequently,912 the total reproductive value of a patch (which corresponds to the class state in our NGM) is strictly proportional to the number913 of mutants it contains, enforcing # ‌ℓ ∝ (1, 2, … ,K)⊤. Thus, calculating # ‌ ℓ⊤S(1) # ‌f ∘ algebraically reduces to evaluating the neutral914 moments ∑k (# ‌f ∘) k and ∑k2 (# ‌f ∘) k , which are fundamentally linked to standard, readily available identity-by-descent (IBD)915 probabilities or relatedness coefficients.916 Under Wright-Fisher demography the local transition matrix is typically dense, so the convenient tri-diagonal recursions917 used above no longer apply; nevertheless, the metapopulation-fitness criterion still reduces to the scalar condition.918 References919 Ajar, É. (2003). Analysis of disruptive selection in subdivided populations. BMC Evol. Biol. , 3.1, p. 22. DOI: 10.1186/1471-2148-3-22.920 Avila, P. & Mullon, C. (2023). Evolutionary game theory and the adaptive dynamics approach: adaptation where individuals interact. Phil.921 Trans. R. Soc. B, 378.1876. DOI: 10.1098/rstb.2021.0502.922 Berman, A. & Plemmons, R. J. (1987). Nonnegative matrices in the mathematical sciences. Classics in applied mathematics 9. Society for Industrial923 and Applied Mathematics.924 Crabtree, D. E. & Haynsworth, E. V. (1969). An identity for the Schur complement of a matrix. Proc. Am. Math. Soc., 22.2, pp. 364–366. DOI:925 10.1090/s0002-9939-1969-0255573-1.926 Frank, S. A. (1998). Foundations of social evolution. Princeton University Press.927 Graham, A. (2018). Kronecker products and matrix calculus with applications . Courier Dover Publications.928 Horn, R. A. & Johnson, C. R. (1985). Matrix Analysis. 169. Cambridge University Press, Cambridge.929 Massol, F., Calcagno, V., & Massol, J. (2009). The metapopulation fitness criterion: Proof and perspectives. Theor. Popul. Biol., 75.2–3, pp. 183–930 200. DOI: http://dx.doi.org/10.1016/j.tpb.2009.02.005.931 Metz, J. A. & Gyllenberg, M. (2001). How should we define fitness in structured metapopulation models? Including an application to the932 calculation of evolutionarily stable dispersal strategies. Proc. R. Soc. B, 268.1466, pp. 499–508.933 Mullon, C., Keller, L., & Lehmann, L. (2016). Evolutionary stability of jointly evolving traits in subdivided populations. Am. Nat., 188.2, pp. 175–934 195. DOI: 10.1086/686900.935 — (2018). Social polymorphism is favoured by the co-evolution of dispersal with social behaviour. Nat. Ecol. Evol., 2.1, p. 132. DOI: 10.1038/936 s41559-017-0397-y.937 Mullon, C. & Lehmann, L. (2014). The robustness of the weak selection approximation for the evolution of altruism against strong selection.938 J. Evol. Biol. , 27.10, pp. 2272–2282. DOI: 10.1111/jeb.12462.939 Roberts, M. G. & Heesterbeek, J. A. P. (2003). A new method for estimating the effort required to control an infectious disease. Proc. R. Soc.940 B, 270.1522, pp. 1359–1364. DOI: 10.1098/rspb.2003.2339.941 Taylor, P. D. & Frank, S. A. (1996). How to make a kin selection model. J. Theor. Biol., 180.1, pp. 27–37. DOI: 10.1006/jtbi.1996.0075.942 Varga, R. S. (1963). Matrix iterative analysis . Graduate Texts in Mathematics. Prentice–Hall, Englewood Cliffs, NJ.943 37 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 25, 2026. ; https://doi.org/10.64898/2026.03.23.713828doi: bioRxiv preprint