Ecoevolutionary dynamics of clonal multicellular life cycles
Abstract
The evolution of multicellular life cycles is a central process in the course of the emergence of multicellularity. The simplest multicellular life cycle is comprised of the growth of the propagule into a colony and its fragmentation to give rise to new propagules. The majority of theoretical models assume selection among life cycles to be driven by internal properties of multicellular groups, resulting in growth competition. At the same time, the influence of interactions between groups on the evolution of life cycles is rarely even considered. Here, we present a model of colonial life cycle evolution taking into account group interactions. Our work shows that the outcome of evolution could be coexistence between multiple life cycles or that the outcome may depend on the initial state of the population – scenarios impossible without group interactions. At the same time, we found that some results of these simpler models remain relevant: evolutionary stable strategies in our model are restricted to binary fragmentation – the same class of life cycles that contains all evolutionarily optimal life cycles in the model without interactions. Our results demonstrate that while models neglecting interactions can capture shortterm dynamics, they fall short in predicting the populationscale picture of evolution.
Editor's evaluation
This article models the evolution of simple multicellular life cycles using evolutionary game theory. The authors discuss natural selection between different life cycles by modeling growth, fragmentation, and interactions between propagules, discovering conditions for selection of a single life cycle or coexistence of multiple ones. Overall, the model is biologically intuitive, the results are rigorous, and the implications for the evolution of multicellularity are interesting.
https://doi.org/10.7554/eLife.78822.sa0Introduction
Multicellular organisms are found everywhere. In all major branches of complex multicellularity (animals, plants, fungi, red and brown algae), organisms are formed by cells staying together after cell division – unlike unicellular species, in which cells part their ways before the next division occurs (MárquezZacarías et al., 2021; Herron et al., 2022). However, organisms have to reproduce as otherwise their species will eventually go extinct. For a multicellular organism, this means that some cells must depart in order to develop into an offspring individual. The combination of organism growth and reproduction constitutes a clonal life cycle. Emergence of clonal multicellular life cycles was the central innovation in the earlier stages of the evolution of multicellularity. There, traits, which do not even exist for unicellular species, become crucial for longterm success of even the most primitive colony of cells (Maynard Smith and Szathmáry, 1995, Michod, 2007). These include the number of cells in the colony, how often cells depart to give rise to new colonies, how large the released propagules are, and how many of them are produced. As the reproduction and, consequently, fitness of simple cell colonies are dependent on these traits, they immediately become subjected to natural selection, favoring some life cycles over others. Since complex multicellular life descends from those loose cell colonies, the understanding of the prior evolution of primitive life cycles is essential to our understanding of the later evolution of complex traits (Ratcliff et al., 2012; De Monte and Rainey, 2014; Hammerschmidt et al., 2014; Doulcier et al., 2020).
There are several theoretical approaches to the modeling of the evolution of multicellular life cycles. The mechanistically simplest class of models assumes that natural selection operates by means of growth competition. Colonies are born small, but due to cell divisions they increase in size and, eventually, fragment, so the number of colonies in the population increases. The life cycle maximizing the population growth rate has a selective advantage as it outgrows all competitors (Roze and Michod, 2001; Libby et al., 2014; Pichugin et al., 2017; Pichugin et al., 2019; Staps et al., 2019; Gao et al., 2019; Pichugin and Traulsen, 2020; Gao et al., 2021; Pichugin, 2022; Pichugin and Traulsen, 2022). For groups made of identical cells, growth competition models of evolution predict that some life cycles cannot be the winners of this growth competition under any conditions. For instance, if the fragmentation event is instantaneous and its execution does not cost anything to the group, only fragmentation into two pieces can evolve (Pichugin et al., 2017; Pichugin and Traulsen, 2020). And indeed, the division into two pieces is, by a large margin, the most common reproductive strategy among microscopic life forms.
However, these models, due to their conceptual simplicity, assume unconstrained (exponential) growth of the population, which cannot be sustained for a prolonged period of time, because resources and space are limited. Other models consider densitydependent growth (Rossetti et al., 2011; Tarnita et al., 2013; van Gestel and Nowak, 2016; Doulcier et al., 2020, Henriques et al., 2021), where the population growth decreases with the number of groups. A similar approach is the Moran birth–death process on the group level, where whenever a new group emerges, one other group dies (Traulsen and Nowak, 2006; Matias Rodrigues et al., 2012; Simon et al., 2013; Luo, 2014; Kaveh et al., 2016; Olejarz et al., 2018; Cooney, 2019; Cooney, 2020). While the population dynamics of densitydependent population growth is vastly different from the exponential explosion found in models of unconstrained growth, these two approaches lead to identical results for life cycle evolution: as shown in Pichugin and Traulsen, 2020, the dynamics of the fraction of a given life cycle in a population are identical in models with unconstrained and densitydependent growth. Therefore, even in models with densitydependent growth, the evolutionary success of the life cycle is still fully determined by the population growth rates.
Nevertheless, densitydependent growth is also a simplification as different groups may differ in their competitiveness. For instance, largecell colonies are able to block single cells from access to vital resources (Rainey and Travisano, 1998; Rainey and Rainey, 2003; Hammerschmidt et al., 2014), which may even lead to a complete extinction of solitary cells. Thus, the population dynamics of multicellular life cycles is not necessarily density dependent, but could be frequency dependent – the impact of resource limitation on the population growth depends on both the size and the composition of the population.
Hence, the evolution of multicellular life cycles cannot always be reduced to growth competition, but may arise from ecoevolutionary dynamics.
From a broader empirical perspective, frequencydependent dynamics is found to be common among microbial populations (Levin, 1988; Ribeck and Lenski, 2015; Healey et al., 2016; Friedman et al., 2017). From the perspective of the theoretical ecology, frequencydependent evolutionary dynamics arising from interactions between diverse population members has also been considered in detail (May, 1972; Wangersky, 1978; Bomze, 1983; Huang et al., 2015; Huang et al., 2017; Bunin, 2017; Barbier et al., 2018; Kotil and Vetsigian, 2018; Tarnita, 2018; Farahpour et al., 2018; Park et al., 2019; Park et al., 2020). The impact of interactions between individuals is recognized in the context of the emergence of aggregative multicellularity, where cells come together to form collectives (Garcia and De Monte, 2013; De Monte and Rainey, 2014; Garcia et al., 2015; Miele and De Monte, 2021). However, both empirical and theoretical ecology approaches tend to overlook frequency dependence in the context of clonal life cycles, where the organism’s growth in the course of the life cycle may cause a change in its role in interactions (but see an example in Tverskoi and Gavrilets, 2022 modeling evolution of germ–soma differentiation).
In this article, we developed a model of the evolution of clonal life cycles under frequencydependent dynamics, implemented in the form of frequencydependent colony death rate. We focus on three questions:
What is the population dynamics of a single life cycle?
What kind of evolutionary outcomes does frequencydependent selection bring?
Are there any patterns or constraints among possible evolutionary outcomes that are universal for multiple forms of frequency dependence?
We first address these questions in the context of simpler models with unconstrained growth. In these models, a population performing any life cycle grows exponentially, the competition between life cycles always results in a single one outcompeting all others (which life cycle will be the winner depends only on the growth/death rates but not on the initial composition of population), and finally, the winner always comes from a limited subset of life cycles (in the simplest version of the model with costless fragmentation – it must be a fragmentation mode producing exactly two offspring). We found that including frequencydependent interactions between organisms performing different life cycles and thus constraining the total population size changes the answers to these questions. First, the population dynamics leads to a stationary state with a finite population size. Second, we found that interactions between groups allow for situations with bistability or coexistence of multiple life cycles – scenarios impossible in the unconstrained growth model. Third, evolutionary stable strategies in our present model always belong to a limited subset of life cycles – the same one containing possible winners of growth competition in models without group interactions. Thus, we found that despite the fundamental differences between our present model and simpler models with unconstrained growth, some of their results have a direct analogy in a much more general ecoevolutionary context considered here.
Model
Population dynamics of a single life cycle
We consider a population consisting of cell groups that grow in size and fragment, giving rise to new groups. Cells within a group of size $i$ divide at rate $b}_{i$, thus a group of size $i$ grows at rate $i{b}_{i}$. Groups also die due to both external environmental factors and withinpopulation competition for resources or space. The death rate of groups of size $i$ due to external factors is $d}_{i$. Frequencydependent competition is modeled as the death of groups of size $i$ upon encounter with groups of size $j$ at rate ${K}_{i,j}$ (see Figure 1A).
Whenever a group of maturity size $m$ grows to $m+1$ cells, it immediately fragments. The fragmentation always occurs by the same pattern and determines the life cycle of a population. We represent a fragmentation pattern by $\kappa $ – a partition of a number $m+1$. For example, the fragmentation pattern of the unicellular life cycle, in which two daughter cells always go apart, is $\kappa =1+1$ (see Figure 1B). Other fragmentation patterns correspond to multicellular life cycles. The simplest of them are the life cycles in which groups grow up to two cells, but fragment upon reaching size three. Such a fragmentation can be performed in two ways: either detachment of a single cell, leading to the fragmentation pattern $\kappa =2+1$, or fission into three solitary cells, $\kappa =1+1+1$ (see Figure 1B). For simplicity, we assume that the cell number does not change during fragmentation (no cell loss), the sum of a fragmentation pattern $\kappa $ is equal to $m+1$.
If we denote the abundance of cell groups containing $i$ cells as $x}_{i$, then the dynamics of population is described by a system of differential equations
where the first two terms $(i1){b}_{i1}{x}_{i1}i{b}_{i}{x}_{i}$ describe the growth of groups – the positive term represents growth from size $i1$ to $i$ and the negative term represents growth from $i$ to $i+1$. The third term ${d}_{i}{x}_{i}$ is the environmentally caused death. The term $m{b}_{m}{\pi}_{i}(\kappa ){x}_{m}$ describes the birth of new groups of size $i$ via fragmentation of larger groups, where ${\pi}_{i}(\kappa )$ is the number of groups of size $i$ produced in the result of that fragmentation. Finally, ${x}_{i}{\sum}_{j=1}^{m}{K}_{i,j}{x}_{j}$ is the death of groups due to the competition between groups.
Summarizing the dynamics into matrix notation, the system (1) can be written as
Here, $\mathbf{x}$ is the column vector of group abundances
The linear term $\mathbf{Ax}$ represents the processes of group growth, fragmentation, and frequencyindependent death. The matrix $\mathbf{A}$ of size $m\times m$ is called a population projection matrix in the field of formal demography – in the sense of the projection of the current state into the future state. For an arbitrary life cycle, matrix $\mathbf{A}$ is given by
The elements of the population projection matrix ${\mathbf{A}}_{i,j}$ represent changes to number of groups of size $i$ due to processes occurring with groups of size $j$ (but not due to interactions). Hence, the population projection matrix has nonzero elements only on the main diagonal (group death and growth of groups to larger sizes), lower subdiagonal (growth of smaller groups), and rightmost column (fragmentation at the end of the life cycle). The elements of the competition matrix $\mathbf{K}$ are given by ${K}_{i,j}$ for $i,j=1,\mathrm{\dots},m$. The operation $\mathrm{diag}(\cdot )$ takes an input vector of length $m$ and transforms it into a diagonal matrix of size $m\times m$ with the entries of the input vector on the main diagonal.
Population dynamics of multiple life cycles
To investigate the ecoevolutionary dynamics of clonal life cycles, we consider a composition of subpopulations executing various life cycles: $\kappa}^{(1)},{\kappa}^{(2)},\dots {\kappa}^{(r)$. In this composite population, the cell growth ($b}_{i$), environmentally caused (constant) group death ($d}_{i$), and group fragmentation (${\pi}_{i}(\kappa )$) occur independently in each subpopulation. However, frequencydependent death due to competition entangles the dynamics of the subpopulations as groups with different life cycles growing together have to compete with each other. If we denote with $x}_{i}^{(j)$ the number of groups containing i cells in a subpopulation executing the life cycle $\kappa}^{(j)$, the dynamics of the composite population is described by the differential equations
where ${m}^{(j)}$ is the maturity size of the life cycle ${\kappa}^{(j)}$, and $n=\mathrm{max}({m}^{(1)},{m}^{(2)},\mathrm{\dots},{m}^{(r)})$ is the maximal maturity size in the composite population. The difference between the one life cycle system (1) and the system of multiple life cycles (5) is in the last term, where groups from every competing subpopulation contribute to frequencydependent death.
In vector form, the state of the composite population is described by a concatenation of vector states of each subpopulation:
where $\stackrel{~}{\mathbf{x}}$ is the column vector describing the state of the composite population, ${\mathbf{x}}^{(j)}$ is the column vector describing $j$th subpopulation in a form (3). Note that the last entries of any ${\mathbf{x}}^{{(j)}^{T}}$ will be zero if ${m}^{(j)}<n$. The dynamics of the composite population in Equation (5) can be represented in the vectorized form, similar to Equation (2):
Here, the composite population projection matrix representing the cell growth, environmentally caused (constant) group death, and group fragmentation is a diagonal block matrix
where ${\mathbf{A}}^{(i)}$ is the population projection matrix of the life cycle ${\kappa}^{(i)}$ extended to size $n\times n$ ($n$ is the maximal maturity size across all competing life cycles). If the maturity size of the life cycle $i$ is ${m}^{(i)}=n$, this matrix has a form exactly as in Equation (4). If the maturity size is smaller, ${m}^{(i)}<n$, then the topleft ${m}^{(i)}\times {m}^{(i)}$ has the form (4), while the remaining elements are nonzero at the main diagonal and the lower subdiagonal, as dictated by Equation (5).
The composite competition matrix $\stackrel{~}{\mathbf{K}}$ is constructed as
where each block $\mathbf{K}$ is a competition matrix.
Invasions from rare
In the general case, the investigation of the composite population dynamics given by Equation (7) is a very complex problem without a general solution. Hence, in our study we consider a specific class of initial conditions – invasion from rare, where the composite population contains only two subpopulations: the abundant ‘resident’ executing life cycle ${\kappa}^{(R)}$ and rare ‘invader’ executing different life cycle ${\kappa}^{(I)}$. This scenario represents an emergence of a mutant in previously stable population of the resident. The population changes if this mutant is capable of invading the resident – otherwise, the mutant goes extinct and the resident population remains the same.
In this scenario, the composite dynamics in Equation (7) can be disentangled into resident and invader components. Since the invader population is small, its contribution to frequencydependent competition is negligible. The members of the resident population compete predominantly between themselves, so the resident dynamics is effectively a singlespecies scenario,
where the vector ${\mathbf{x}}^{(R)}$ represents the composition of the resident population, ${\mathbf{A}}^{(R)}$ is the population projection matrix of the resident, ${\mathbf{x}}^{(R)*}$ is the equilibrium composition, and we introduced the selfinvasion population projection matrix ${\mathbf{A}}^{(R,R)}={\mathbf{A}}^{(R)}\mathrm{diag}({\mathbf{Kx}}^{(R)*})$. Since the resident is assumed to be at a stable equilibrium in the absence of invaders, the selfinvasion matrix ${\mathbf{A}}^{(R,R)}$ has an eigenvalue that is zero, and the equilibrium population composition ${\mathbf{x}}^{(R)}$ is given by the corresponding eigenvector. The resident population dynamics can be obtained by solving the nonlinear Equation (10), which in the general case can be performed only numerically.
The rare invader population also competes primarily with the resident and selfcompetition is negligible. Thus, its dynamics is given by
where vector ${\mathbf{x}}^{(I)}$ represents the composition of the invader population, and we introduced the invasion matrix
Unlike the resident dynamics, the dynamics of the invader population is linear – the invasion population projection matrix ${\mathbf{A}}^{(I,R)}$ is independent from the invader population state ${\mathbf{x}}^{(I)}$. The linear dynamics of clonal life cycles has been extensively studied in previous work (Pichugin et al., 2017). If the largest eigenvalue of the invasion matrix ${\mathbf{A}}^{(I,R)}$ is positive, then the invader population will increase in numbers, independently of its initial demography. Otherwise, the invader population goes extinct.
The assumption of a negligible impact of the invader population on competition limits the analysis to the early stages of invasion, when the invader population is small. Nevertheless, this makes it possible to investigate the stability of resident life cycles with respect to invasions.
Results
We first briefly recap the population dynamics and evolution under a more basic model with unconstrained growth (${K}_{ij}=0$) (Pichugin et al., 2017; Pichugin and Traulsen, 2020). This model has three main features. First, a population executing a single life cycle grows exponentially in the long run. The population growth rate is given by the leading eigenvalue of the population projection matrix ($\mathbf{A}$) and the demographic composition is given by the associated eigenvector. Second, selection always finds a single winner. In a composite population, where different subpopulations execute different life cycles, only one life cycle survives in the long run – the one that has the largest growth rate. This outcome is independent of the initial state of the population. Third, some life cycles cannot be optimal under any combination of cell division rates ($b}_{i$) and group death rates ($d}_{i$). In the simplest case of instant and costless group fragmentation, life cycles with more than two offspring cannot win the growth competition.
Next, we consider how these features transfer to a system taking into account competition between groups. We begin with the dynamics of a single life cycle (section ‘Dynamics of a single life cycle’) in a system with a population size constraint, which is very different from exponential growth. Then, we consider how the competition proceeds in our model: for two (section ‘Competition between pairs of life cycles may result in coexistence or bistability’) and multiple life cycles (section ‘Competition between multiple life cycles’), where a rich spectrum of possible stationary states is found. After that, we outline the restrictions imposed on evolutionary stable strategies (section ‘The set of possible evolutionary stable strategies is restricted’). We conclude with presenting scenarios of a special interest: interactions with killer and victim kernels, where the evolutionary dynamics is reduced to the competition for growth rate and carrying capacity respectively (section ‘Killer and victim kernels guarantee a dominance of a single life cycle’), and investigate the competition between unicellular and multicellular life cycles (section ‘Conditions promoting the evolution of multicellular life cycles’).
Dynamics of a single life cycle
For the simplest unicellular life cycle ($\kappa =1+1$), where population is composed only of solitary cells, the dynamics of our model given by Equation (2) reduces to logistic growth,
where the net growth rate is equal to ${b}_{1}{d}_{1}$ and the carrying capacity is $({b}_{1}{d}_{1})/{K}_{1,1}$. The characteristic feature of logistic growth is that the population approaches the carrying capacity with time, starting from either high or low populations size (see Figure 2A). The population dynamics of more complex life cycles also bears similarity to the logistic growth of the unicellular life cycle. If a population is small, the competition term is small, so the population grows exponentially. While population size increases, so does the competition term – group mortality rises and the overall population growth slows down. The growth stops when the population reaches a stationary state ${\mathbf{x}}^{*}$, where
Numerical simulations show that a population executing a single life cycle always comes to the same stationary state ${\mathbf{x}}^{*}$ from any initial distribution of group sizes (see Figure 2B and C).
Competition between pairs of life cycles may result in coexistence or bistability
A composite population containing several subpopulations executing different life cycles also reaches a stationary state. In the simplest case, only one life cycle survives, while others go extinct. However, we found that the stationary state may contain more than one life cycle (coexistence). Also, the stationary state may depend on the initial state of the population (multistability). Neither of these can occur in the linear model without competition.
To illustrate these effects, we focus on a pair of life cycles (${\kappa}^{(1)},{\kappa}^{(2)}$) with the special initial conditions, where one of these life cycles is abundant, while the other one is rare. The life cycle ${\kappa}^{(1)}$ can invade from rare into the abundant ${\kappa}^{(2)}$ if the largest eigenvalue of the invasion matrix ${\mathbf{A}}^{(1,2)}$ is positive; see Equation (12). Otherwise, the rare life cycle ${\kappa}^{(1)}$ goes extinct. For a pair of life cycles, there are four possible scenarios of invasion from rare (see the outline in Table 1 and corresponding dynamics in Figure 3).
Life cycle ${\kappa}^{(1)}$ dominates life cycle ${\kappa}^{(2)}$ if ${\kappa}^{(1)}$ spreads from rare, but ${\kappa}^{(2)}$ does not. This is equivalent to the leading eigenvalue of the invasion matrix ${\mathbf{A}}^{(1,2)}$ being positive, while the leading eigenvalue of ${\mathbf{A}}^{(2,1)}$ is negative. Then, independently of the initial conditions, only the life cycle ${\kappa}^{(1)}$ survives in the long run (see Figure 3B). The opposite signs result in the dominance of ${\kappa}^{(2)}$ over ${\kappa}^{(1)}$ (see Figure 3C).
Beyond a dominance of one life cycle over another, it is possible that each life cycle is able to spread in the abundance of another. This happens when the leading eigenvalues of both invasion matrices ${\mathbf{A}}^{(1,2)}$ and ${\mathbf{A}}^{(2,1)}$ are positive. There, the result of interactions between life cycles is a coexistence of both – an outcome impossible in the model without competition (see Figure 3D).
Finally, the leading eigenvalues of both invasion matrices could be negative – then neither of the life cycles can invade into another. Then, the result is a bistability between life cycles, where the outcome of interactions depends on the initial conditions – another result impossible in the model without competition (see Figure 3A).
The competition between groups plays a major role in determining which of the four invasion patterns occurs. For instance, it is possible that a life cycle having an advantage in the raw growth rate (i.e., dominating in the unconstrained growth model) is dominated by the result of the competition (see Figure 3C), where the life cycle 1 + 1 has a faster growth but 2 + 1 still dominates due to the advantage of multicellular groups in competition.
Competition between multiple life cycles
Extending our analysis beyond just a pair of life cycles, we considered the evolutionary dynamics in a population with multiple of them. We numerically investigated the evolutionary dynamics in a population containing all life cycles in which groups do not exceed a size of three cells. There are seven such life cycles: unicellular (1 + 1), two with bicellular groups (2 + 1 and 1 + 1 + 1), and four with groups of three cells (3 + 1, 2 + 2, 2 + 1 + 1, and 1 + 1 + 1 + 1) (see Figure 4A). We generated a set of 13,000 randomly drawn combinations of growth and competition rates from an exponential distribution with unit rate parameter. For each set, we simulated 100 independent replicates of population dynamics that differ in the initial conditions (abundances and demographic composition of each of the seven life cycles).
These runs were classified by the state of population at the end of simulation (see Figure 4B and Table 2). The majority of simulations (about 75%) resulted in the survival of a single life cycle across all 100 replicates. The next most common outcome is a coexistence between two or more life cycles – found in about 17% of simulations. Also, a multistability between two or three life cycles was observed in about 6.5% of simulations. Here, coexistence describes a situation where the stationary state of the population is composed of the same set of at least two life cycles in every replicate. Multistability describes a situation, where in every replicate the stationary state contained only a single life cycle, but there were different stationary states among the replicates. More complex outcomes were also observed – these were the compositions of multistability and coexistences, which contributed to only about 1.5% of simulations. The most common of these composite situations is a minimal combination of bistability and coexistence: there are two possible stationary states, one is a single life cycle and another is a coexistence of two other life cycles (0.7% of simulations).
The numerical investigation of evolutionary dynamics of multiple life cycles revealed that a diverse range of outcomes are possible, including multistability and coexistence of life cycles, as well as their combinations. At the same time, the most common result is still the dominance of a single life cycle.
The set of possible evolutionary stable strategies is restricted
Between simple dominance and multistability, about 80% of evolutionary simulations ended with the survival of a single life cycle. This happens if a life cycle is an evolutionary stable strategy (ESS), for example, if it is abundant, it cannot be invaded by any other life cycle. In the basic model without competition, the life cycle with the maximal growth rate also satisfies the definition of an ESS. Here, we found that the set of evolutionary stable strategies is similarly restricted – it also contains only fragmentations into exactly two pieces.
To show this restriction, we consider a triplet of life cycles ${\kappa}^{(1)}$, ${\kappa}^{(2)}$, ${\kappa}^{(3)}$. If the life cycle ${\kappa}^{(1)}$ is a resident, there are four variants of its stability against invasions: (i) either it is stable against invasions from both ${\kappa}^{(2)}$ and ${\kappa}^{(3)}$ (then ${\kappa}^{(1)}$ is an ESS), or (ii) stable against only ${\kappa}^{(2)}$, or (iii) stable against only ${\kappa}^{(3)}$, or (iv) both ${\kappa}^{(2)}$ and ${\kappa}^{(3)}$ can invade. Similar four variants exist for the two other life cycles. As a result, for the whole triplet, there are ${4}^{3}=64$ possible pairwise invasion patterns, which could feature 0, 1, 2, or 3 evolutionary stable strategies.
Numerical simulations show that all 64 patterns can be expressed for the same triplet of life cycles for some combination of cell birth, group death, and competition rates (see Figure 5A). We generated a set of 40,000 randomly drawn combinations of these rates from an exponential distribution with unit rate parameter and analyzed the pairwise invasion patterns for each. The six most frequent patterns, comprising 77% of the generated dataset, feature a hierarchical dominance, where the life cycles can be ordered in a way that higherorder life cycle dominates (always invade) lowerorder life cycles. These six patterns are all possible hierarchical dominance triplets as there are exactly six ways how three items can be placed in order. If we use the same analysis for the basic, linear model with unrestricted growth, we will only observe hierarchical dominance as larger growth rate results in domination there. On the opposite side of the frequency spectrum, the two most rare patterns feature cyclic dominance, together comprising only 0.015% of the dataset. There, in any pair of life cycles one dominates another but the whole triplet follows a ‘rock–paper–scissors’ dynamics with no evolutionary stable strategies present (cf. Park et al., 2020).
While an arbitrary triplet of life cycles may demonstrate up to 64 invasion patterns, some triplets, which we will call ‘constrained,’ feature much smaller diversity of patterns. A triplet is constrained if the fragmentation rule of one (constrained) life cycle ${\kappa}^{(M)}$ can be represented as a combination of fragmentation rules of two other (constraining) life cycles ${\kappa}^{(C1)}$ and ${\kappa}^{(C2)}$. The simplest example is the triplet ${\kappa}^{(C1)}=2+1$, ${\kappa}^{(C2)}=1+1$, ${\kappa}^{(M)}=1+1+1$, where the fragmentation of a threecelled group into three solitary cells ($3\to 1+1+1$) can be presented as a combination of the detachment of a single cell ($3\to 2+1$) immediately followed by the dissolution of the twocell group ($2\to 1+1$). A lot of constrained triplets can be constructed, for example, in our illustrations we use the triplet ${\kappa}^{(C1)}=2+2,{\kappa}^{(C2)}=4+4,{\kappa}^{(M)}=4+2+2$. Originally, constrained triplets emerged in the model with unconstrained growth (Pichugin et al., 2017), where the growth rate of the constrained life cycle ${\kappa}^{(M)}$ was found to be always between growth rates of the constraining life cycles ${\kappa}^{(C1)}$ and ${\kappa}^{(C2)}$. It follows that in the present model with competition (${K}_{ij}\ne 0$), each of two constraining life cycles (${\kappa}^{(C1)}$ and ${\kappa}^{(C2)}$) must be either stable against two others or unstable against both. The constrained life cycle ${\kappa}^{(M)}$ in turn is always invaded by one of constraining life cycles and is stable against the other (see Appendix 1). Hence, the number of possible pairwise invasion patterns for such a triplet is limited to $2\cdot 2\cdot 2=8$ (see Figure 5B). Among these eight patterns, two feature hierarchical dominance, where the life cycles can be ordered in a way that a higherorder life cycle dominates the lowerorder life cycles (with the constrained life cycle ${\kappa}^{(M)}$ always being in the middle of hierarchy) (see Figure 6A). In a larger dataset (66,000 entries) with random birth, death, and competition rates, hierarchical dominance was observed in about 78% of entries. Two patterns feature bistability between constraining life cycles ${\kappa}^{(C1)}$ and ${\kappa}^{(C2)}$ (see Figure 6B). These patterns appear in about 11% of the dataset. Two more patterns feature a coexistence between all three life cycles (see Figure 6C). Note that unlike a pair of life cycles with unique coexistence equilibrium considered in section ‘Competition between pairs of life cycles may result in coexistence or bistability,’ the triplet features a range of stable coexistence states. The coexistence pattern is similarly frequent, observed in 11% of the dataset. Finally, the least frequent two patterns are nonhierarchical dominance, where one constraining life cycle dominates another, but in the abundance of the constrained life cycle, the invasion pattern is inverse (see Figure 6D). They appear with three orders of magnitude lower frequency, smaller than 0.01% of cases.
The fundamental feature of constrained triplets is that any constrained life cycle (${\kappa}^{(M)}$) can always be invaded by exactly one constraining life cycles (${\kappa}^{(C1)}$ and ${\kappa}^{(C2)}$). Hence, any life cycle, which can be constrained by two others, cannot be an ESS. We found that any life cycle with more than two offspring can be constrained (see Appendix 1 for the proof). As a result, only binary fragmentation life cycles can be evolutionary stable strategies.
Killer and victim kernels guarantee a dominance of a single life cycle
As shown above, the evolutionary dynamics of interacting life cycles can be quite complex. In this context, a special interest attracts these cases, where the complexity of dynamics is limited. In this section, we present two forms of interaction matrices ($\mathbf{K}$), which guarantee that the evolutionary outcome is a straightforward domination of a single life cycle.
The first special case is the killer kernel, defined as
There, the probability of a group to die in an encounter depends only on the size of the opponent group ($j$), hence the name.
For an arbitrary killer kernel, a single life cycle has the same demographic composition in the stationary state as it has in the nointeractions model (${K}_{ij}=0$) (see Appendix 2 for the proof).
This feature of the demography leads to the result that the outcome of evolution under a killer kernel is also the same as in the nointeraction model. To show that, consider a composite population containing multiple life cycles; its dynamics is described by Equation (7) and depends on the composite projection ($\stackrel{~}{\mathbf{A}}$) and competition ($\stackrel{~}{\mathbf{K}}$) matrices. If the competition matrix $\mathbf{K}$ is a killer kernel, the composite competition matrix $\stackrel{~}{\mathbf{K}}$ defined in Equation (9) is a killer kernel as well. Since the population dynamics of both single and multiple life cycles are governed by equations with the same structure (Equation (2) and (7), respectively), the results of Appendix 2 carry over to a composite population. Specifically, the stationary state of the composite population is proportional to the eigenvector of the composite population projection matrix $\stackrel{~}{\mathbf{A}}$ corresponding to its leading eigenvalue. The composite population projection matrix $\stackrel{~}{\mathbf{A}}$ defined in Equation (8) is a block diagonal matrix, composed of population projection matrices of the life cycles constituting the composite population. Thus, the leading eigenvalue of $\stackrel{~}{\mathbf{A}}$ is the largest among the eigenvalues of all population projection matrices comprising $\stackrel{~}{\mathbf{A}}$. Additionally, the corresponding eigenvector has nonzero components only at the positions in $\stackrel{~}{\mathbf{x}}$ associated with the block having the maximal leading eigenvalue – that is, only that life cycle constitutes the stationary state. This rule is equivalent to the choice of the fastest growing life cycle in the linear model. Thus, the evolution of life cycles competing by a killer kernel (${K}_{ij}={k}_{j}$) can be reduced to growth competition, which results in the survival of only one life cycle.
Another special case is the victim kernel defined as
There, the chance of a group to die depends only on the size of that group ($i$). For an arbitrary victim kernel, the carrying capacity of a single life cycle can be explicitly found. The total number of groups at the stationary state ${N}^{*}={\sum}_{i}{x}_{i}^{*}$ is equal to the growth rate of this life cycle in the nointeraction model (${K}_{ij}=0$) with modified cell division rates ${b}_{i}^{\prime}={b}_{i}/{k}_{i}$ and group death rates ${d}_{i}^{\prime}={d}_{i}/{k}_{i}$ (see Appendix 3 for the proof).
We found that in the composite population with many life cycles interacting by a victim kernel, only one survives in the long run (see the proof in Appendix 3). In such a scenario, each life cycle grows in numbers if the total population size is smaller than the carrying capacity of that life cycle (derived in Appendix 3). If the whole population size exceeds this value, the life cycle gradually dies out. Hence, selection favors the life cycle with the largest carrying capacity because it can grow in a dense population, when all other life cycles die from intense competition.
In both killer and victim kernels, the evolutionary dynamics is reduced to optimization of a single trait: growth rate for the killer kernel and carrying capacity for the victim kernel. Thus, the result of selection in these scenarios is a dominance of a single life cycle over all suboptimal competitors.
Conditions promoting the evolution of multicellular life cycles
We conclude our findings with one of the most biologically interesting cases – the evolution of multicellular life cycles in a population dominated by unicellular organisms. This is also one of the simplest scenarios as it deals with the simplest, unicellular life cycle 1 + 1.
If a unicellular population is abundant in the system, its stationary state can be explicitly found (see Appendix 4). Hence, for an arbitrary invading multicellular life cycle ${\kappa}^{(M)}$, its invasion matrix (${\mathbf{A}}^{(M,1)}$) can be found explicitly. As a result, a multicellular invader can spread from rare in a population of unicellular residents if this life cycle has a positive growth rate in a model with unconstrained growth with modified group death rates,
The successful multicellular invader drives the unicellular life cycle to extinction if the unicellular life cycle cannot invade from rare. If the unicellular life cycle is an invader, its invasion matrix (${\mathbf{A}}^{(1,M)}$) has size $1\times 1$ and the invasion rate is just equal to its only element. The unicellular life cycle cannot invade into a multicellular resident life cycle if
where ${m}^{(M)}$ is the maximal group size in the life cycle ${\kappa}^{(M)}$, and ${x}_{i}^{(M)*}$ is the composition of the population when ${\kappa}^{(M)}$ is abundant. Derivations of both conditions are presented in Appendix 4.
Discussion
In this article, we have added an ecological dimension to the problem of life cycle evolution. Specifically, we are interested in three aspects of the ecoevolutionary dynamics: (i) What is a population dynamics of a single life cycle? (ii) What is the evolutionary dynamics of multiple life cycles? (iii) What are the possible constraints on the outcomes of that evolution? All three questions has been extensively studied for the model with unconstrained growth (Pichugin et al., 2017), and it is natural to contrast our current findings with these results.
First, the introduction of group competition completely changes the population dynamics of the life cycle growth. Instead of the exponential explosion of the population size occurring in the model with unconstrained population size, the growth in our current model slows down with the population size. Eventually, the population approaches a stationary state with the limited population size and constant composition (see Figure 2B).
Second, frequencydependent selection arising from competition allows diverse outcomes of evolution: the stationary state can feature a coexistence of multiple life cycles or multistability between several possible end states, as shown in Figure 3. By contrast, evolution in the model with unconstrained growth always results in a survival of a single life cycle, independently of the initial conditions.
Third, despite all the differences in the population and evolutionary dynamics, the restrictions on the possible evolutionary stable strategies are exactly the same in the models with and without population size constraints (see section ‘The set of possible evolutionary stable strategies is restricted’). In both models, an ESS may only feature a fragmentation into exactly two pieces. Any life cycle producing more than two offspring can always be invaded by another life cycle with a smaller number of offspring.
Beyond the scope of our model, life cycles with fragmentation into multiple pieces may become evolutionary stable strategies if fragmentation is costly (e.g., imposes a risk of cell loss). There, binary fragmentation is no longer special as a fragmentation into multiple pieces can win growth rate competition. Nevertheless, constrained triplets of life cycles still exist under a costly fragmentation scenario but the constraining condition is different from the one presented here. There, a life cycle containing two different subsets of offspring with the same combined size is constrained between two life cycles, which use either of these subsets twice (Pichugin and Traulsen, 2020). For instance, ${\kappa}^{(M)}=2+1+1$ is constrained as it has different offspring subsets 2 and $1+1$ with the same combined size. This life cycle is constrained between life cycles ${\kappa}^{(C1)}=2+2$ and ${\kappa}^{(C2)}=1+1+1+1$, which use one of these subsets twice. Since the scenario of costly fragmentation still contains constrained triplets, life cycles satisfying the rule above, such as 2 + 1 + 1, cannot be evolutionary stable strategies there. Note that this rule works in the present model too (so there are actually two ways to construct constraining triplets), but it is a weaker condition than the rule presented in section ‘The set of possible evolutionary stable strategies is restricted’ and does not allow to construct any additional constraining triplets.
In the broad context of the ecoevolutionary dynamics, our dynamical Equations (2) and (7) bear a similarity with the generalized Lotka–Volterra (GLV) equations widely used in ecology: both contain a linear growth term and a nonlinear competition term (typically of the second order) balancing out the linear growth. However, our equations are not equivalent to the GLV. In the GLV, individuals corresponding to different elements of the population vector reproduce independently, that is, in our terms, the population projection matrix $\mathbf{A}$ is diagonal. In our model, however, an individual group changes its state in the course of the life cycle and the population projection matrix $\mathbf{A}$ is not diagonal. Of course, if the population projection matrix $\mathbf{A}$ can be diagonalized, we can perform a linear transformation $\mathbf{x}\to \mathbf{Cy}$ ($\mathbf{C}$ is a matrix) to make the linear term in our model diagonal as in GLV. However, in this case, the interaction term will lose the GLV form of the modification of the growth rate (${x}_{i}{\sum}_{j}{K}_{ij}{x}_{j}$, see Equation (1)) and will become a general secondorder term instead (${\sum}_{jk}{K}_{ij}{x}_{j}{x}_{k}$). Given that our system is not a GLV in disguise, it is surprising how much of the analysis presented here has been performed using the approaches designed to analyze GLV systems.
In this article, we found that the competition plays a major role in the evolution of multicellular life cycles. Our choice of the competition matrix values was driven by theoretical aspects of this article: we either use randomly drawn values for numerical simulations or choose specific forms leading to analytical results. Yet, competition in natural populations is neither random nor finetuned to mathematically beautiful outcomes. What might empirical competition matrices in a stagestructured population look like? Unfortunately, the demographics of simple multicellular species is not sufficiently studied experimentally. Still, we can consider an example of emergence of Pseudomonas fluorescens colonies, where a competition plays a major role in the population dynamics and evolution (Rainey and Travisano, 1998; Rainey and Rainey, 2003; Hammerschmidt et al., 2014). In a still liquid media, these initially unicellular bacteria evolve a ‘glue’ production, which causes formation of multicellular aggregates. These aggregates float atop the media, gaining an exclusive access to oxygen. Once the entire surface is covered by continuous bacterial biofilm, the unicellular phenotype living in the body of the media is completely denied the oxygen access and dies out. In the framework of our study, the competition matrix $\mathbf{K}$ is determined by the capability to block oxygen and surface area access. Naturally, the more cells there are in the group, the more stress they apply to others and, at the same time, the more resistant to competition they are. In the limit of small population size, single cells have almost no impact on others (${K}_{i,1}\ll 1$) and are the most susceptible to oxygen denial (${K}_{1,i}\gg 1$). In opposite limit, an established mat of millions of cells just drives everything else to extinction (${K}_{i,j\gg 1}\gg 1$). For arbitrary competitors size, the terms ${K}_{ij}$ should increase with the size of an opponent group ($K}_{ij}<{K}_{i,j+1$) but decrease with the size of the focal group ($K}_{ij}>{K}_{i+1,j$). Similar competition matrices should arise in scenarios where being bigger is better.
Appendix 1
Invasion of life cycles and restrictions on ESS
In this section, we show that any resident life cycle with fragmentation into multiple parts can always be invaded by at least one life cycle with a smaller number of offspring.
Consider a resident life cycle ${\kappa}^{(R)}$, in which more than two offspring groups are produced as a result of fragmentation. The initial dynamics of any life cycle ${\kappa}^{(I)}$ invading from rare can be described by a linear model with death rates modified as
The invasion is successful if the leading eigenvalue of the corresponding population projection matrix ${\mathbf{A}}^{(I,R)}$, defined as in Equation (11), is positive.
In the analysis of the linear model (Pichugin et al., 2017), it was shown that if the fragmentation preserves the number of cells (no cell loss), then for any multiple fragmentation life cycle, there exist two constraining life cycles with a smaller number of offspring. For any combination of cell division rates (b_{i}) and group death rates (d_{i}), one of the constraining life cycles has a larger growth rate than the focal multiple fragmentation life cycle, while another has a smaller growth rate.
Now consider the invasion from rare in our present model with competition. The initial invasion rate computed to the leading eigenvalue of the matrix ${\mathbf{A}}^{(I,R)}$ is equal to the growth rate of the invader life cycle in an environment modified according to Equation (19). Therefore, for any resident population, the invasion rate of the constrained life cycle is always in between the invasion rates of the constraining life cycles. If the resident population is formed by the constrained life cycle itself, its selfinvasion rate is zero. Hence, one of the constraining life cycles has a larger invasion rate (i.e., it is positive), while another has a smaller invasion rate (negative). As a result, the constrained life cycle can always be invaded by exactly one of its constraining life cycles. No constrained life cycle can be an ESS. Since any life cycle with more than two offspring is constrained, only binary fragmentation can be an ESS.
To conclude Appendix 1, we consider the resident population formed by a constraining life cycle. If the constrained life cycle has a positive invasion rate, then another constraining life cycle must have a positive invasion rate as well. Alternatively, the constrained life cycle has a negative invasion rate, then another constraining life cycle also has a negative invasion rate. Thus, a constraining resident is either resistant to invasion from both other life cycles in a triplet or can be successfully invaded by both of them.
Population dynamics under the killer kernel
In this section, we show that the demographic distribution of populations in a stationary regime is identical in the linear model (${K}_{i,j}=0$) and under the killer kernel (${K}_{i,j}={k}_{j}$).
First, we consider a linear model, where the population dynamics is governed by the population projection matrix $\mathbf{A}$:
A number of useful properties of this dynamics comes from the Perron–Frobenius theorem for irreducible nonnegative matrices. Here, nonnegative matrix means a matrix with all elements greater or equal to zero. Matrix $\mathbf{M}$ is irreducible if its representation by a directed graph, where node $i$ has an edge to node $j$ only if ${M}_{ij}>0$, is a strongly connected graph, that is, there is a path from any node to any other node. Since the states $i$ and $j$ represent groups of different sizes, and the population projection matrix $\mathbf{A}$ describes the dynamics of the system, the irreducibility of $\mathbf{A}$ means that group of any given size $i$ can reach any other size $j$ through growth and fragmentation.
However, the population projection matrix $\mathbf{A}$ itself is not nonnegative as its main diagonal is negative (${A}_{ii}=i{b}_{i}{d}_{i}$). It is not always irreducible as well: if the fragmentation mode produces only multicellular offspring (e.g., $\kappa =2+2$), no group can ever reach unicellular state $j=1$. Both these limitations can be resolved. First, since the negative elements are located only on the main diagonal, the matrix ${\mathbf{A}}^{\prime}=\mathbf{A}+\mathbf{I}\cdot \mathrm{max}(i{b}_{i}+{d}_{i})$ is nonnegative. The matrix ${\mathbf{A}}^{\prime}$ has the same eigenvectors as $\mathbf{A}$, and its eigenvalues (${\lambda}^{\prime}$) relate to eigenvalues of the population projection matrix ($\lambda $) as ${\lambda}^{\prime}=\lambda +\mathrm{max}(i{b}_{i}+{d}_{i})$. Second, the irreducibility arises only when there are groups sizes smaller than the smallest produced offspring. The number of groups of these sizes will continuously decrease in time: they will grow into larger sizes and spontaneously die. Since fragmentation is not able to resupply these small groups, the population will eventually get rid of them, so these small sizes can be discarded from consideration of the longterm dynamics. Thus, without loss of generality, we can truncate the population projection matrix to take into account only the sizes actually emerging in a life cycle. The resulting modified matrix is nonnegative and irreducible, hence, the Perron–Frobenius theorem applies. From this theorem, it immediately follows that for arbitrary birth rates (b_{i}), death rates (d_{i}), and the life cycle executed ($\kappa $), the population projection matrix has a simple leading eigenvalue $\lambda $, its eigenspace is onedimensional, all components of the eigenvector are positive, and no other eigenvectors of this matrix have all their components positive.
Therefore, in the linear model, the stationary regime is an exponentially growing population (Pichugin et al., 2017)
where $\lambda $ is the leading eigenvalue of the matrix $\mathbf{A}$, and ${\mathbf{W}}^{*}$ is the corresponding eigenvector,
Under the killer kernel, the death rates due to competition are the same for groups of all sizes. Hence, following Equation (2), the population dynamics of a life cycle under a killer kernel is
where 1 is the vector of ones, and $\mathbf{I}$ is the identity matrix ($\mathrm{diag}(1)=\mathbf{I}$). It can be shown by a direct calculation that the vector
is the stationary state of the dynamics (23):
Note that while any other eigenvector ${\mathbf{W}}^{\prime}$ of the population projection matrix $\mathbf{A}$ can be used in Equation (25), only the leading eigenvector ${\mathbf{W}}^{*}$ can represent a biologically meaningful population, as all other eigenvectors have negative elements, which would mean negative number of groups of some size at the stationary state.
Since the stationary state ${\mathbf{x}}^{*}$ under the killer kernel is proportional to the vector describing the stationary distribution ${\mathbf{W}}^{*}$ in the linear model, the population compositions in both scenarios are the same.
Population dynamics under the victim kernel
Dynamics of a single life cycle
In this section, we show that the task of finding the equilibrium population size (${N}^{*}={\sum}_{i}{x}_{i}^{*}$) under the victim kernel (${K}_{i,j}={k}_{i}$) is mathematically equivalent to the task of finding the population growth rate in the linear model (${K}_{i,j}=0$) with modified cell birth rates (${b}_{i}\to {b}_{i}/{k}_{i}$) and group death rates (${d}_{i}\to {d}_{i}/{k}_{i}$).
In the linear model, the population growth rate is found as the leading eigenvalue of the population projection matrix determined by Pichugin et al., 2017
Under the victim kernel, the death rate due to the competition depends only on the size of the outcompeted group. Hence, following Equation (2), the stationary state of a life cycle under a victim kernel satisfies
where $\mathbf{k}$ is a vector constructed from any column of the competition matrix $\mathbf{K}$ (they are all identical), and ${N}^{*}={\sum}_{i}{x}_{i}^{*}$ is the equilibrium population size.
The last equality in Equation (27) implies that by Fredholm alternative one out of two conditions is satisfied (Hoffman, 1971):
Limiting ourselves to scenarios where the stationary state ${\mathbf{x}}^{*}$ is not an empty population, we can conclude that the population at the equilibrium satisfies
where in the last step, we divided the $i$th column by k_{i} for all $i$. Comparing the determinants in Equations (26) and (29), we find that they are identical after the substitution
Thus, the equilibrium population size (${N}^{*}$) under a victim kernel can be found as a population growth rate in the linear model with modified cell birth and group death rates.
Competition of multiple life cycles
In this section, we show that in a composite population, in which multiple ($r$) life cycles compete through a victim kernel (${K}_{i,j}={k}_{i}$), only one life cycle survives to the stationary state, and this is the same life cycle that has the maximal carrying capacity if grown in isolation, and the same life cycle that would be evolutionarily optimal in the linear model (${K}_{i,j}=0$) with modified cell division (${b}_{i}\to {b}_{i}/{k}_{i}$) and group death (${d}_{i}\to {d}_{i}/{k}_{i}$) rates.
If the competition matrix $\mathbf{K}$ constitutes a victim kernel, then the composite competition matrix $\stackrel{~}{\mathbf{K}}$ defined as in Equation (9) is a victim kernel as well. Hence, our results from section ‘Dynamics of a single life cycle’ hold also for a composite population. In particular, the population size at the stationary state (${\stackrel{~}{N}}^{*}={\sum}_{p=1}^{r}{\sum}_{i=1}^{n}{x}_{i}^{(p)}$) is determined by
where $\stackrel{~}{\mathbf{k}}$ is the vector made from a single column of the composite competition matrix $\stackrel{~}{\mathbf{K}}$ (it is a concatenation of $r$ vectors $\mathbf{k}$). In the second step of Equation (31), we used the property of blockdiagonal matrices: the determinant of a blockdiagonal matrix is the product of determinants of its blocks. Equation (31) is satisfied if one of the multipliers is equal to zero. The $p$th multiplier becomes zero when the composite population reaches a size equal to the carrying capacity of the $p$th life cycle; cf. Equation (28). At that moment, the $p$th life cycle can neither grow or decay, the life cycles with lower carrying capacity decay as they cannot keep up with competition caused by overcrowding, and only the life cycles with carrying capacity larger than the population size can grow in numbers. Equation (31) has $r$ possible solutions with respect to ${\stackrel{~}{N}}^{*}$ – one for each competing life cycle. However, only one solution can represent a stationary population, where no life cycle can grow in numbers – the one with the maximal population size. There, one life cycle is stationary, while all others decrease in numbers due to overcrowding. Thus, the outcome of life cycle competition under a victim kernel is the survival of a single life cycle, which has the maximal equilibrium population size among all competitors. According to section ‘Dynamics of a single life cycle,’ these population sizes are equal to the growth rates of the corresponding life cycles in a linear model with the modified cell division and group death rates. Consequently, the maximal population size corresponds to the fastest growing life cycle in that modified linear model. The population corresponding to the highest eigenvalue takes over and will dominate the system. This means that the life cycle with the largest population size in isolation dominates over all other life cycles in the competition through a victim kernel.
Invasion into unicellular resident and invasion of the unicellular invader
If the resident is unicellular (${\kappa}^{(R)}=1+1$), its steady state is given by the solution of
equal to
Then, for an arbitrary invader multicellular life cycle ${\kappa}^{(M)}$, the invasion matrix is given by
This is equivalent to the linear growth of the invader life cycle in an environment with modified death rates
If the invader is unicellular (${\kappa}^{(I)}=1+1$), then the invasion matrix is effectively a $1\times 1$ matrix since the invader contains only isolated cells. Even if ${\mathbf{A}}^{(I,M)}$ formally has a larger size, it is a block matrix, with the element ${\mathbf{A}}_{1,1}^{(I,M)}$ being a block, with a value that is equal to the growth rate of the unicellular life cycle. This value is
where ${x}_{i}^{(M)}$ is the number of groups of size $i$ in the resident population, and ${m}^{(M)}$ is the maximal group size of the resident life cycle. The unicellular invader cannot spread in a resident population, when ${\mathbf{A}}_{1,1}^{(I,M)}<0$.
Parameters of calculation used in figures
Figure 2
In Figure 2A, we used ${b}_{1}=1,{d}_{1}=0,{K}_{11}=1$. The initial population size was drawn from the uniform distribution $U(0.1,2)$.
In Figure 2B, we used $b=(1,1)$ and $d=(0,0)$. The competition matrix was
The initial number of groups of each size was randomly drawn from the uniform distribution $U(0.1,2)$.
Figure 3
In Figure 3, we considered the life cycles 1 + 1 and 2 + 1. We used $b=(1,0.5)$ and $d=(0,0)$. For each plotted trajectory, the populations were initialized with ${\mathbf{x}}_{1+1}=({s}_{1})$, ${\mathbf{x}}_{2+1}=({s}_{2},0)$, where ${s}_{1},{s}_{2}\in \{0.1,0.2,0.3,\mathrm{\dots},1.0\}$. The dynamics shown in the four panels differ by the competition matrix used:
Panel A
$K=\left(\begin{array}{cc}\hfill 1\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 1\hfill \end{array}\right)$Panel B
$K=\left(\begin{array}{cc}\hfill 3\hfill & \hfill 3\hfill \\ \hfill 1\hfill & \hfill 1\hfill \end{array}\right)$Panel C
$K=\left(\begin{array}{cc}\hfill 1\hfill & \hfill 0.1\hfill \\ \hfill 0.1\hfill & \hfill 1\hfill \end{array}\right)$Panel D
$K=\left(\begin{array}{cc}\hfill 1\hfill & \hfill 1\hfill \\ \hfill 0.6\hfill & \hfill 0.4\hfill \end{array}\right)$
Figure 6
In Figure 6, we considered the life cycles ${\kappa}^{(C1)}=2+2$, ${\kappa}^{(C2)}=4+4$, ${\kappa}^{(M)}=4+2+2$. Panels differ in birth, death, and competition rates. Trajectories in each panel have different initial states. For each initial state, the composite population contains all three life cycles with different fractions ${s}^{(C1)},{s}^{(C2)},{s}^{(M)}$, such that ${s}^{(C1)}+{s}^{(C2)}+{s}^{(M)}=1$. The initial group sizes distribution is proportional to the equilibrium population of that life cycle alone:
where the vectors ${\mathbf{x}}^{\mathbf{\ast}\mathbf{(}\mathbf{C}\mathbf{1}\mathbf{)}}}^{T$, ${\mathbf{x}}^{\mathbf{\ast}\mathbf{(}\mathbf{C}\mathbf{2}\mathbf{)}}}^{T$, ${\mathbf{x}}^{\mathbf{\ast}\mathbf{(}\mathbf{M}\mathbf{)}}}^{T$ are equilibrium population states of corresponding life cycles grown in isolation, computed according to Equation (14).
In panel A (hierarchic dominance), we used ${b}_{i}=1.0$, ${d}_{i}=0$ and the competition matrix ${K}_{ij}=0.1$:
or equivalently
In panel B (bistability), we used ${b}_{i}=1$ and ${d}_{i}=0$, while the competition matrix was
or equivalently
In panel C (coexistence), we used ${b}_{i}=1$ and ${d}_{i}=0$, while the competition matrix was
or equivalently
In panel D (nonhierarchical dominance), we used ${d}_{i}=0$,
and
Data availability
The code and data are available from https://github.com/yuriypichugin/Ecoevolutionarydynamicslifecycles, (copy archived at swh:1:rev:b7122aeec7b11953867e6b5e588701e5a602276d).
References

Ecological communities with LotkaVolterra dynamicsPhysical Review. E 95:042414.https://doi.org/10.1103/PhysRevE.95.042414

The replicator dynamics for multilevel selection in evolutionary gamesJournal of Mathematical Biology 79:101–154.https://doi.org/10.1007/s00285019013525

Analysis of multilevel replicator dynamics for general twostrategy social dilemmaBulletin of Mathematical Biology 82:66.https://doi.org/10.1007/s1153802000742x

Nascent multicellular life and the emergence of individualityJournal of Biosciences 39:237–248.https://doi.org/10.1007/s1203801494205

Community structure follows simple assembly rules in microbial microcosmsNature Ecology & Evolution 1:109.https://doi.org/10.1038/s415590170109

Interacting cells driving the evolution of multicellular life cyclesPLOS Computational Biology 15:e1006987.https://doi.org/10.1371/journal.pcbi.1006987

Group formation and the evolution of socialityEvolution; International Journal of Organic Evolution 67:131–141.https://doi.org/10.1111/j.15585646.2012.01739.x

Games of multicellularityJournal of Theoretical Biology 403:143–158.https://doi.org/10.1016/j.jtbi.2016.04.037

Emergence of evolutionarily stable communities through ecoevolutionary tunnellingNature Ecology & Evolution 2:1644–1653.https://doi.org/10.1038/s4155901806557

Frequencydependent selection in bacterial populationsPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 319:459–472.https://doi.org/10.1098/rstb.1988.0059

Geometry shapes evolution of early multicellularityPLOS Computational Biology 10:e1003803.https://doi.org/10.1371/journal.pcbi.1003803

A unifying framework reveals key properties of multilevel selectionJournal of Theoretical Biology 341:41–52.https://doi.org/10.1016/j.jtbi.2013.09.024

Why have aggregative multicellular organisms stayed simple?Current Genetics 67:871–876.https://doi.org/10.1007/s00294021011930

Aggregative cycles evolve as a solution to conflicts in social investmentPLOS Computational Biology 17:e1008617.https://doi.org/10.1371/journal.pcbi.1008617

Selection for synchronized cell division in simple multicellular organismsJournal of Theoretical Biology 457:170–179.https://doi.org/10.1016/j.jtbi.2018.08.038

Fragmentation modes and the evolution of life cyclesPLOS Computational Biology 13:e1005860.https://doi.org/10.1371/journal.pcbi.1005860

Evolution of simple multicellular life cycles in dynamic environmentsJournal of the Royal Society, Interface 16:20190054.https://doi.org/10.1098/rsif.2019.0054

Evolution of multicellular life cycles under costly fragmentationPLOS Computational Biology 16:e1008406.https://doi.org/10.1371/journal.pcbi.1008406

Modeling and quantifying frequencydependent fitness in microbial populations with crossfeeding interactionsEvolution; International Journal of Organic Evolution 69:1313–1320.https://doi.org/10.1111/evo.12645

Emergent multicellular life cycles in filamentous bacteria owing to densitydependent population dynamicsJournal of the Royal Society, Interface 8:1772–1784.https://doi.org/10.1098/rsif.2011.0102

Mutation, multilevel selection, and the evolution of propagule size during the origin of multicellularityThe American Naturalist 158:638–654.https://doi.org/10.1086/323590

Towards a general theory of group selectionEvolution; International Journal of Organic Evolution 67:1561–1572.https://doi.org/10.1111/j.15585646.2012.01835.x

Emergence of diverse life cycles and life histories at the origin of multicellularityNature Ecology & Evolution 3:1197–1205.https://doi.org/10.1038/s4155901909400

Evolutionary construction by staying together and coming togetherJournal of Theoretical Biology 320:10–22.https://doi.org/10.1016/j.jtbi.2012.11.022

Fast evolution unlocks forbidden communitiesNature Ecology & Evolution 2:1525–1526.https://doi.org/10.1038/s415590180688y

The evolution of germsoma specialization under different genetic and environmental effectsJournal of Theoretical Biology 534:110964.https://doi.org/10.1016/j.jtbi.2021.110964

Phenotypic heterogeneity and the evolution of bacterial life cyclesPLOS Computational Biology 12:04764.https://doi.org/10.1371/journal.pcbi.1004764

Lotkavolterra population modelsAnnual Review of Ecology and Systematics 9:189–218.https://doi.org/10.1146/annurev.es.09.110178.001201
Decision letter

Babak MomeniReviewing Editor; Boston College, United States

Aleksandra M WalczakSenior Editor; CNRS LPENS, France

Denis TverskoiReviewer
Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Decision letter after peer review:
Thank you for submitting your article "Ecoevolutionary dynamics of clonal multicellular life cycles" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Aleksandra Walczak as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Denis Tverskoi (Reviewer #3).
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission. Overall, the reviewers found the manuscript informative and interesting, but have the following suggestions to improve the presentation and clarity of the manuscript. In addition to the consolidated essential revisions listed below, we would highly recommend that you address all the reviewers' comments in your revision.
Essential revisions:
1) Please comment on the impact of the mechanism of interaction and discuss the forms of interactions that might be more relevant in the modeling. In the discussions, it would also be helpful to mention real biological systems that the models represent.
2) Please revise the terminology (e,g, "projection matrix", "linear dynamics", and stationary state") to avoid confusion with technical math definitions.
3) In the discussion "Similarities between models with and without group competition," please include results from models without competition as well to make the comparisons clearer for the reader.
4) Please include additional quantitative discussions and elaborate on the effect of kernel on the outcomes.
5) Please highlight the primary takeaways of the paper among many results that are derived to underscore the main findings. The reviewers suggest that you frame these questions at the beginning of your paper to help your reader follow the modeling and results.
6) Please reorganize the Discussion section to include a recap and improve its flow for your readers.
7) Please doublecheck the equations for accuracy and address the concerns raised by reviewer #3.
Reviewer #1 (Recommendations for the authors):
I have few comments and suggestions to improve the quality of the work:
1. The model is very clearly stated and a connection between current model and previous works has been established. The generalization of interaction between different groups are quite interesting and significant improvement in the modeling approach to the evolution of multicellular life cycles. However, it would be more instructive if the authors identify what form of interactions might be more relevant in the modeling, beside victim and killer. Also, it would be pedagogical if the above is compared, in a couple of sentences or a paragraph, with other frameworks in evolutionary game theory. To me, the model in the absence of interactions is a quasispecies equation where each group is a trait and every growth terms is similar to mutation to a new trait. A given life cycle is a mutational pathway in this analogy. The interaction term is interaction between the different types/traits. It would be instructive if authors, by making such connections with theoretical evolution literature, discuss what qualitative or qualitative changes they already expect from introducing the interactions kernel.
2. The model has a complex form, inevitably. There are several mechanisms interplaying with each other. Growth rates of groups, fragmentation patterns, interaction between groups of the same life cycle, and interaction between groups belonging to different life cycles and most often the latter two are indistinguishable. Even though the model is complex, the paper is well written and the results are clearly enough stated. However, it would be very important to improve the significance of the work by highlighting a few main messages of the paper among many results that are derived to highlight the main findings. This is the part I found somewhat unclear in the current manuscript.
3. Regarding the above comment, I want to suggest that the authors highlight, from the beginning, what questions they are trying to answer and what quantities can be calculated from the model. In principle, I expect the model says the condition for selection of a lifecycle against other ones (invasion of rare) and how it relates to model parameters. Similarly conditions for coexistence between different life cycles. Finally, the steady state values of each group abundances. These are all discussed in the paper to some extent. I just suggest highlighting the main quantities you are calculating and questions you are aiming to answer a little clearer.
4. In the case of two competing life cycles, it is not clear, to me, what is the qualitative result. The kernel is introduced in four forms suggested by the authors (Figure 2). It is concluded that depending on whether one has a killer kernel or victim kernel one life cycle wins and in other cases there might be coexistence. Are there further results? For example can one in principle states the condition for ESS or selection advantage of one life cycle?
5. (Similar to the previous comment) How does different kernels (victim or killer) improve or modify condition for selection of one life cycle while competing with another? The answer is detailed but it is useful to sum up the overall behavior cycles.
6. I want to comment that the three life cycle competition results was interesting and the analogy with rockpaperscissor game was clarifying the findings.
7. Does the introduction of kernel term overall promotes coexistence between two life cycles?
8. While the paper has a theory approach and the generality of results justifies it, it would be very useful to connect the general finding with the observations in experimental evolution of multicellularity. For example, discuss cases where there are unique life cycles and when there are more than one.
The evolution of multicellular life cycles is a central process in the course of the emergence of multicellularity. The model suggested by authors connects evolution of multicellular life cycles to evolutionary game theory. The introduction of the interaction terms seems to be a great modeling way to discuss under what circumstances different life cycles can coexist or when one life cycle is chosen among other potential ones due to a natural selection among life cycles. The results are discussed in some details but due to complexity in some cases examples are used. I recommend this work for publication in Journal eLife after revision. It is a new model with high impact on the field of evolution of multicellular structures.
Reviewer #2 (Recommendations for the authors):
1) The first paragraph does not contain any references to back up the claims presented by the authors. It would be good to address this and add appropriate references.
2) In Figure 1 it would be good to explicitly state in the caption what the different colors for the groups mean.
3) In the discussion "Similarities between models with and without group competition "it would be good to explicitly refer to results from results from models without competition (which I think are only seen in Figure 2A) to make the comparisons clearer for the reader.
Reviewer #3 (Recommendations for the authors):
Please find my detailed recommendations for the authors below
Page 7, line 142: “The projection matrix $A$…” It is not clear to me what the authors mean by “projection matrix”. From a mathematical point of view, a square matrix $P$ is called a projection matrix if it is equal to its square, i.e. if $P^2=P$. However, it does not seem that in general $A^2 = A$.
Page 10, line 219: “Unlike the resident dynamics, the dynamics of the invader population is linear…” This dynamics is linear only if $x^{(R)}$ is constant (e.g., at equilibrium). Therefore, it would be helpful to use different notations for the function $x^{(R)}(t)$ and its value at the equilibrium $x^{(R,*)}$ in formulas 1112 describing these dynamics.
Page 11, lines 245246: “Numerical simulations show that an isolated life cycle always comes to the same stationary state $x^*$ from any initial distribution of group sizes.” Could you please add more details about these numerical simulations? What kind of life cycles were considered? What birth and death rates have been chosen? How was the competition matrix generated?
Throughout the text, the authors discuss the results on various special types of a competition matrix (constant, killer kernel, victim kernel). These special cases are used to obtain analytical conclusions that cannot be drawn in a general case. This helps readers deepen their intuition about this model. However, I wonder if there are any real biological systems that can be described by these special kernels? If so, it would be very helpful to include the corresponding discussion in the text.
Page 25, line 515: “In the linear model, the stationary state is an exponentially growing population…” From a mathematical point of view, a stationary state of a system is a state with all observables independent of time. Therefore, I am not sure if it is appropriate to use the term “stationary state” here.
Page 25, lines 526529. It is not clear what does the index $s$ mean in the corresponding formulas.
Appendix A1. As was shown in the previous work, the solution to Equation (14) depends on the leading eigenvalue $\λ^*$ of the matrix $A$ and the corresponding eigenvector $w^*$. After refreshing that, the authors showed that a stationary state of the dynamics governed by Formula (21) (i.e., by Equation (14) in the case of the killer kernel) depends on $\λ^*$ and $w^*$ (see Formula (22)). This is correct. However, for each eigenvalue $\λ$ and each corresponding eigenvector $w$, Formula (22) produces a stationary state of the dynamics governed by Formula (21). It is not clear to me why, in the case of the killer kernel, the authors consider only the leading eigenvector? If such a conclusion is made on the basis of the above results for the linear case, please explain why these results can be generalized to the killer kernel.
And the last question. Is it possible that the leading eigenvalue $\λ$ has an eigenspace of dimension higher than 1 so that stationary states (22) form a “line” of equilibria?
Page 25, line 528. I guess, it should be $x_{s,j}^*$ instead of $x_{s,j}$.
Page 27, line 551552. What is $N_0$? How is it related to $N^*$?
Section 3.2. The authors state that the dynamics are very complex in the case of competition between multiple life cycles, and therefore consider these dynamics only in some special cases. I agree that the dynamics are complex. However, as a reader, I have no intuition about how the model works in the general case, which is the most interesting question for me. Therefore, it would be helpful to add some numerical simulations exploring the above dynamics in the general case. It would also be useful to present some statistics illustrating the average number of different life cycles presented in a stationary state as well as the number of extinct ones. Maybe some life cycles are rarely observed in a stationary state, while others are widespread under a broad range of parameters and initial conditions?
Page 16, Figure 3D. If possible, could you please mark in different colours trajectories approaching different stationary states?
Lines 450451. “Yet, it is possible to introduce a linear transformation $x \rightarrow Cy$, where $C$ is a matrix, which will make the linear term in our model diagonal.” Where is a proof that the matrix $A$ can be diagonalized?
I think the structure of the Discussion section could be improved. For example, I do not understand why you put paragraph 2 (lines 431442) at the beginning of the Discussion? This paragraph is about a generalization of your model and its results, and it is strange to discuss the generalization of the results before discussing the results themselves. Second, it might be helpful to add a brief recap of the problem under study, and a brief overview of the model at the beginning of the Discussion section before discussing the position of your study in various contexts, and related results (lines 443505).
https://doi.org/10.7554/eLife.78822.sa1Author response
Essential revisions:
1) Please comment on the impact of the mechanism of interaction and discuss the forms of interactions that might be more relevant in the modeling. In the discussions, it would also be helpful to mention real biological systems that the models represent.
We have now described the interactions in more detail and extended the discussion. Concerning the biology, we draw inspiration from several microbial systems, but we do not make detailed predictions for any particular system – in part, because so far there are no experiments considering the emergence and competition of multiple life cycles at the same time. However, these experiments have now become feasible and some experimentalists have already expressed interest.
2) Please revise the terminology (e,g, "projection matrix", "linear dynamics", and stationary state") to avoid confusion with technical math definitions.
We apologize for having used a terminology that is potentially misleading. These terms are commonly used in the field of mathematical demography, but we now switched to a more precise notation.
3) In the discussion "Similarities between models with and without group competition," please include results from models without competition as well to make the comparisons clearer for the reader.
We agree that our results are hard to understand without the background provided by that basic model. We now included a paragraph describing the results established for that model without group competition.
4) Please include additional quantitative discussions and elaborate on the effect of kernel on the outcomes.
We have now included more details on the interaction kernels.
5) Please highlight the primary takeaways of the paper among many results that are derived to underscore the main findings. The reviewers suggest that you frame these questions at the beginning of your paper to help your reader follow the modeling and results.
We agree that it is a good idea to mention the main takeaway already in the introduction and have now done so.
6) Please reorganize the Discussion section to include a recap and improve its flow for your readers.
We have rewritten the discussion along these lines.
7) Please doublecheck the equations for accuracy and address the concerns raised by reviewer #3.
These are concerns that we took very seriously, but we show in the revised manuscript that our previous approach was justified. For example, we can focus on the leading eigenvector as it is the only one that is biologically relevant, as all its elements are positive.
Reviewer #1 (Recommendations for the authors):
I have few comments and suggestions to improve the quality of the work:
1. The model is very clearly stated and a connection between current model and previous works has been established. The generalization of interaction between different groups are quite interesting and significant improvement in the modeling approach to the evolution of multicellular life cycles. However, it would be more instructive if the authors identify what form of interactions might be more relevant in the modeling, beside victim and killer. Also, it would be pedagogical if the above is compared, in a couple of sentences or a paragraph, with other frameworks in evolutionary game theory. To me, the model in the absence of interactions is a quasispecies equation where each group is a trait and every growth terms is similar to mutation to a new trait. A given life cycle is a mutational pathway in this analogy. The interaction term is interaction between the different types/traits. It would be instructive if authors, by making such connections with theoretical evolution literature, discuss what qualitative or qualitative changes they already expect from introducing the interactions kernel.
We agree that the general dynamics of a population subdivided into several classes with transitions between them by a set of differential equations does resemble the quasispecies framework. However, the interpretation of our model is very different: For example, in our model each organism passes through different stages from birth to death. In the quasispecies model, an individual maintains its identity and switching between classes occurs due to mutations. In addition to the main connection to theoretical biology in terms of demographic models, we now have added connections to other related topics.
2. The model has a complex form, inevitably. There are several mechanisms interplaying with each other. Growth rates of groups, fragmentation patterns, interaction between groups of the same life cycle, and interaction between groups belonging to different life cycles and most often the latter two are indistinguishable. Even though the model is complex, the paper is well written and the results are clearly enough stated. However, it would be very important to improve the significance of the work by highlighting a few main messages of the paper among many results that are derived to highlight the main findings. This is the part I found somewhat unclear in the current manuscript.
We have now concluded the introduction with our main question and with the key results.
3. Regarding the above comment, I want to suggest that the authors highlight, from the beginning, what questions they are trying to answer and what quantities can be calculated from the model. In principle, I expect the model says the condition for selection of a lifecycle against other ones (invasion of rare) and how it relates to model parameters. Similarly conditions for coexistence between different life cycles. Finally, the steady state values of each group abundances. These are all discussed in the paper to some extent. I just suggest highlighting the main quantities you are calculating and questions you are aiming to answer a little clearer.
Thank you. Agreed, we now spell out more clearly what these conditions mean and how they arise.
4. In the case of two competing life cycles, it is not clear, to me, what is the qualitative result. The kernel is introduced in four forms suggested by the authors (Figure 2). It is concluded that depending on whether one has a killer kernel or victim kernel one life cycle wins and in other cases there might be coexistence. Are there further results? For example can one in principle states the condition for ESS or selection advantage of one life cycle?
Unfortunately, we have no results for the fully general case. However, even the cases where the interaction kernel does not alter the competition are potentially of a lot of interest – it implies that we can derive results from exponential growth competition and they carry over to restricted systems. We now describe this more clearly.
5. (Similar to the previous comment) How does different kernels (victim or killer) improve or modify condition for selection of one life cycle while competing with another? The answer is detailed but it is useful to sum up the overall behavior cycles.
Thank you for this question. We reworked the results structure and put the presentation of killer and victim kernels to the dedicated subsection 3.5. There, we show that killer and victim kernels are the special cases, where the result of evolution is a single life cycle dominating others. For a killer kernel, the winning life cycle is the one with the largest growth rate. For a victim kernel, the winning life cycle is the one with the largest carrying capacity.
6. I want to comment that the three life cycle competition results was interesting and the analogy with rockpaperscissor game was clarifying the findings.
We agree that it is interesting, but it is not very frequent. Now we also refer to the ecological literature on nontransitive interactions, which is another way to phrase this.
7. Does the introduction of kernel term overall promotes coexistence between two life cycles?
No, not necessarily, as for many kernels there is no coexistence. But without such a kernel, there is no generic scenario for coexistences in the first place.
8. While the paper has a theory approach and the generality of results justifies it, it would be very useful to connect the general finding with the observations in experimental evolution of multicellularity. For example, discuss cases where there are unique life cycles and when there are more than one.
Unfortunately, the experimental literature on the evolution of life cycles is only about to take off. So far, we only know that some experimental scientists start to work on experiments in this direction, but preliminary results from Will Ratcliff lab in Georgia Tech indicate that evolution of colonial yeast may lead to coexistence of two life cycles with different colony sizes  a scenario predicted by our model.
The evolution of multicellular life cycles is a central process in the course of the emergence of multicellularity. The model suggested by authors connects evolution of multicellular life cycles to evolutionary game theory. The introduction of the interaction terms seems to be a great modeling way to discuss under what circumstances different life cycles can coexist or when one life cycle is chosen among other potential ones due to a natural selection among life cycles. The results are discussed in some details but due to complexity in some cases examples are used. I recommend this work for publication in Journal eLife after revision. It is a new model with high impact on the field of evolution of multicellular structures.
Thank you.
Reviewer #2 (Recommendations for the authors):
1) The first paragraph does not contain any references to back up the claims presented by the authors. It would be good to address this and add appropriate references.
Thanks, now we added appropriate references.
2) In Figure 1 it would be good to explicitly state in the caption what the different colors for the groups mean.
Done. Here, colors code for number of cells in an organism. Thanks.
3) In the discussion "Similarities between models with and without group competition "it would be good to explicitly refer to results from results from models without competition (which I think are only seen in Figure 2A) to make the comparisons clearer for the reader.
We understand that the paper was very hard to read without an explicit discussion of these cases, which we have now added.
Reviewer #3 (Recommendations for the authors):
Please find my detailed recommendations for the authors below
Page 7, line 142: “The projection matrix $A$…” It is not clear to me what the authors mean by “projection matrix”. From a mathematical point of view, a square matrix $P$ is called a projection matrix if it is equal to its square, i.e. if $P^2=P$. However, it does not seem that in general $A^2 = A$.
True. Thank you for pointing out this source of possible misunderstanding. We took this term from the literature of formal demography, where the word “projection” uses the meaning of extrapolation into future. Now, we use more precise term “population projection matrix” and added additional explanation where it is introduced.
Page 10, line 219: “Unlike the resident dynamics, the dynamics of the invader population is linear…” This dynamics is linear only if $x^{(R)}$ is constant (e.g., at equilibrium). Therefore, it would be helpful to use different notations for the function $x^{(R)}(t)$ and its value at the equilibrium $x^{(R,*)}$ in formulas 1112 describing these dynamics.
Agreed, it also makes the notation in Equations 11 and 12 consistent with the rest of the paper.
Page 11, lines 245246: “Numerical simulations show that an isolated life cycle always comes to the same stationary state $x^*$ from any initial distribution of group sizes.” Could you please add more details about these numerical simulations? What kind of life cycles were considered? What birth and death rates have been chosen? How was the competition matrix generated?
We added an illustration of how different initial states converge to the same steady state (Figure 2B,C) for a single life cycle and a chosen set of birth/death/interactions rates. We also additionally checked the convergence to the same steady state with the dataset obtained for Figure 4B – with 7 competing life cycles, 13000 randomly sampled sets of rates (all drawn independently from exponential distribution with unit rate parameter), and 100 replicates of evolution dynamics for each set different by the initial conditions. 75% of 13000 sets featured a dominance of a single life cycle. In all these cases, the variance in the final number of groups across replicates of a set was less than 10^(4). Thus, which a high accuracy, the dynamics of a single life cycle comes to the same stationary state.
Throughout the text, the authors discuss the results on various special types of a competition matrix (constant, killer kernel, victim kernel). These special cases are used to obtain analytical conclusions that cannot be drawn in a general case. This helps readers deepen their intuition about this model. However, I wonder if there are any real biological systems that can be described by these special kernels? If so, it would be very helpful to include the corresponding discussion in the text.
We speculate about possible forms of empirical interaction matrices in the discussion, but so far the experimental field does not look into such interaction kernels. But we hope that this kind of data is available in a few years!
Page 25, line 515: “In the linear model, the stationary state is an exponentially growing population…” From a mathematical point of view, a stationary state of a system is a state with all observables independent of time. Therefore, I am not sure if it is appropriate to use the term “stationary state” here.
Agreed, we should use steady state here instead (note that demographers speak of stable demographies when populations are growing exponentially in size, so it depends on the field which words seem to be appropriate to characterize such dynamics).
Page 25, lines 526529. It is not clear what does the index $s$ mean in the corresponding formulas.
Thanks, this was a typo left from the earlier draft. Index $s$ does not represent anything and is removed.
Appendix A1. As was shown in the previous work, the solution to Equation (14) depends on the leading eigenvalue $\λ^*$ of the matrix $A$ and the corresponding eigenvector $w^*$. After refreshing that, the authors showed that a stationary state of the dynamics governed by Formula (21) (i.e., by Equation (14) in the case of the killer kernel) depends on $\λ^*$ and $w^*$ (see Formula (22)). This is correct. However, for each eigenvalue $\λ$ and each corresponding eigenvector $w$, Formula (22) produces a stationary state of the dynamics governed by Formula (21). It is not clear to me why, in the case of the killer kernel, the authors consider only the leading eigenvector? If such a conclusion is made on the basis of the above results for the linear case, please explain why these results can be generalized to the killer kernel.
Right, any eigenvector of the population projection matrix can work in this equation. However, only the eigenvector associated with the leading eigenvalue has positiveonly elements. This follows from PerronFrobenius theorem (we added a paragraph explaining why this theorem applies to our matrices). All other eigenvectors contain negative elements, which would correspond to negative number of groups at the stationary state, so these would not be a biologically meaningful solutions.
And the last question. Is it possible that the leading eigenvalue $\λ$ has an eigenspace of dimension higher than 1 so that stationary states (22) form a “line” of equilibria?
No, such a situation is forbidden by the PerronFrobenius theorem as well.
Page 25, line 528. I guess, it should be $x_{s,j}^*$ instead of $x_{s,j}$.
Yes, thanks!
Page 27, line 551552. What is $N_0$? How is it related to $N^*$?
Thank you, this was a typo – it refers to $N^*$!
Section 3.2. The authors state that the dynamics are very complex in the case of competition between multiple life cycles, and therefore consider these dynamics only in some special cases. I agree that the dynamics are complex. However, as a reader, I have no intuition about how the model works in the general case, which is the most interesting question for me. Therefore, it would be helpful to add some numerical simulations exploring the above dynamics in the general case. It would also be useful to present some statistics illustrating the average number of different life cycles presented in a stationary state as well as the number of extinct ones. Maybe some life cycles are rarely observed in a stationary state, while others are widespread under a broad range of parameters and initial conditions?
Right, we did not have this intuition either and the question you raise is an interesting one! We now added new section 3.3. considering the competition between all life cycles with groups smaller than 4 cells (there are 7 such life cycles). We classified the stationary states: most often we find a dominance of a single life cycle, like in a linear model. However, in rare cases, we observe quite complex outcomes – e.g. bistability between different combinations of coexisting life cycles.
Page 16, Figure 3D. If possible, could you please mark in different colours trajectories approaching different stationary states?
Thank you, that is a very good idea!
Lines 450451. “Yet, it is possible to introduce a linear transformation $x \rightarrow Cy$, where $C$ is a matrix, which will make the linear term in our model diagonal.” Where is a proof that the matrix $A$ can be diagonalized?
We have no such proof but it is not needed. Here we say that even if $A$ can be diagonalized, then our system is still not equivalent to generalized LotkaVolterra due to the complications occurring with competition term. If $A$ cannot be diagonalized, as you suggested, then our system is not equivalent to GLV for sure, and our message still holds. In the revision, we clarified this sentence.
I think the structure of the Discussion section could be improved. For example, I do not understand why you put paragraph 2 (lines 431442) at the beginning of the Discussion? This paragraph is about a generalization of your model and its results, and it is strange to discuss the generalization of the results before discussing the results themselves. Second, it might be helpful to add a brief recap of the problem under study, and a brief overview of the model at the beginning of the Discussion section before discussing the position of your study in various contexts, and related results (lines 443505).
Thank you, we have now rewritten the discussion along the suggestions of you and the other reviewers.
https://doi.org/10.7554/eLife.78822.sa2Article and author information
Author details
Funding
No external funding was received for this work.
Senior Editor
 Aleksandra M Walczak, CNRS LPENS, France
Reviewing Editor
 Babak Momeni, Boston College, United States
Reviewer
 Denis Tverskoi
Version history
 Preprint posted: March 16, 2022 (view preprint)
 Received: March 21, 2022
 Accepted: August 9, 2022
 Version of Record published: September 13, 2022 (version 1)
Copyright
© 2022, Ress et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 499
 Page views

 148
 Downloads

 1
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
The microbial community composition in the human gut has a profound effect on human health. This observation has lead to extensive use of microbiome therapies, including overthecounter 'probiotic' treatments intended to alter the composition of the microbiome. Despite so much promise and commercial interest, the factors that contribute to the success or failure of microbiometargeted treatments remain unclear. We investigate the biotic interactions that lead to successful engraftment of a novel bacterial strain introduced to the microbiome as in probiotic treatments. We use pairwise genomescale metabolic modeling with a generalized resource allocation constraint to build a network of interactions between taxa that appear in an experimental engraftment study. We create induced subgraphs using the taxa present in individual samples and assess the likelihood of invader engraftment based on network structure. To do so, we use a generalized LotkaVolterra model, which we show has strong ability to predict if a particular invader or probiotic will successfully engraft into an individual's microbiome. Furthermore, we show that the mechanistic nature of the model is useful for revealing which microbemicrobe interactions potentially drive engraftment.

 Computational and Systems Biology
Angiogenesis is a morphogenic process resulting in the formation of new blood vessels from preexisting ones, usually in hypoxic microenvironments. The initial steps of angiogenesis depend on robust differentiation of oligopotent endothelial cells into the Tip and Stalk phenotypic cell fates, controlled by NOTCHdependent cell–cell communication. The dynamics of spatial patterning of this cell fate specification are only partially understood. Here, by combining a controlled experimental angiogenesis model with mathematical and computational analyses, we find that the regular spatial Tip–Stalk cell patterning can undergo an order–disorder transition at a relatively high input level of a proangiogenic factor VEGF. The resulting differentiation is robust but temporally unstable for most cells, with only a subset of presumptive Tip cells leading sprout extensions. We further find that sprouts form in a manner maximizing their mutual distance, consistent with a Turinglike model that may depend on local enrichment and depletion of fibronectin. Together, our data suggest that NOTCH signaling mediates a robust way of cell differentiation enabling but not instructing subsequent steps in angiogenic morphogenesis, which may require additional cues and selforganization mechanisms. This analysis can assist in further understanding of cell plasticity underlying angiogenesis and other complex morphogenic processes.