Eco-evolutionary dynamics of clonal multicellular life cycles

  1. Vanessa Ress
  2. Arne Traulsen
  3. Yuriy Pichugin  Is a corresponding author
  1. Max Planck Institute for Evolutionary Biology, Germany
  2. Hamburg Center for Health Economics, University of Hamburg, Germany
  3. Department of Ecology and Evolutionary Biology, Princeton University, United States

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 short-term dynamics, they fall short in predicting the population-scale 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.sa0

Introduction

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árquez-Zacarí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 long-term 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 density-dependent 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 density-dependent 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 density-dependent growth. Therefore, even in models with density-dependent growth, the evolutionary success of the life cycle is still fully determined by the population growth rates.

Nevertheless, density-dependent growth is also a simplification as different groups may differ in their competitiveness. For instance, large-cell 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 eco-evolutionary dynamics.

From a broader empirical perspective, frequency-dependent 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, frequency-dependent 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 frequency-dependent dynamics, implemented in the form of frequency-dependent colony death rate. We focus on three questions:

  • What is the population dynamics of a single life cycle?

  • What kind of evolutionary outcomes does frequency-dependent 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 frequency-dependent 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 eco-evolutionary 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 bi, thus a group of size i grows at rate ibi. Groups also die due to both external environmental factors and within-population competition for resources or space. The death rate of groups of size i due to external factors is di. Frequency-dependent competition is modeled as the death of groups of size i upon encounter with groups of size j at rate Ki,j (see Figure 1A).

Model of clonal life cycles.

(A) There are four processes occurring in the model: groups grow by cell division, which occurs with rates bi, groups spontaneously die with rates di, groups fragment immediately upon exceeding maturity size m (3 in this example), according to predefined fragmentation pattern κ (1 + 1 + 1 + 1 here), and groups of size i die due to the competition with groups of size j at rates Kij. (B) Together, group growth and fragmentation constitute a life cycle, in which initially small groups grow in size and eventually fragment, giving rise to more small groups. Colors represent group sizes: red, solitary cells; green, bicellular groups; blue, tricellular groups.

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 κ – 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 κ=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 κ=2+1, or fission into three solitary cells, κ=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 κ is equal to m+1.

If we denote the abundance of cell groups containing i cells as xi, then the dynamics of population is described by a system of differential equations

(1) dx1dt=b1x1d1x1+mbmπ1(κ)xmx1j=1mk1,jxj,dxidt|i>1=(i1)bi1xi1ibixidixi+mbmπi(κ)xmxij=1mk1,jxj,

where the first two terms (i1)bi1xi1ibixi describe the growth of groups – the positive term represents growth from size i-1 to i and the negative term represents growth from i to i+1. The third term -dixi is the environmentally caused death. The term mbmπi(κ)xm describes the birth of new groups of size i via fragmentation of larger groups, where πi(κ) is the number of groups of size i produced in the result of that fragmentation. Finally, -xij=1mKi,jxj is the death of groups due to the competition between groups.

Summarizing the dynamics into matrix notation, the system (1) can be written as

(2) dxdt=Axdiag(Kx)x.

Here, x is the column vector of group abundances

(3) x=(x1,x2,x3,,xm)T.

The linear term Ax represents the processes of group growth, fragmentation, and frequency-independent death. The matrix A of size m×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 A is given by

(4) A=(-b1-d1000mbmπ1(κ)b1-2b2-d200mbmπ2(κ)02b2-3b3-d30mbmπ3(κ)003b30mbmπ4(κ)000(m-1)bm-1mbmπm(κ)-mbm-dm).

The elements of the population projection matrix Ai,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 K are given by Ki,j for i,j=1,,m. The operation diag() takes an input vector of length m and transforms it into a diagonal matrix of size m×m with the entries of the input vector on the main diagonal.

Population dynamics of multiple life cycles

To investigate the eco-evolutionary dynamics of clonal life cycles, we consider a composition of subpopulations executing various life cycles: κ(1),κ(2),κ(r). In this composite population, the cell growth (bi), environmentally caused (constant) group death (di), and group fragmentation (πi(κ)) occur independently in each subpopulation. However, frequency-dependent 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 xi(j) the number of groups containing i cells in a subpopulation executing the life cycle κ(j), the dynamics of the composite population is described by the differential equations

(5) dx1(j)dt=b1x1(j)d1x1(j)+m(j)bm(j)π1(κ(j))xm(j)(j)x1(j)p=1rk=1nK1,kxk(p),dxi(j)dt|i>1=(i1)bi1xi1(j)ibixi(j)dixi(j)+m(j)bm(j)πi(κ(j))xm(j)(j)xi(j)p=1rk=1nKi,kxk(p),

where m(j) is the maturity size of the life cycle κ(j), and n=max(m(1),m(2),,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 frequency-dependent death.

In vector form, the state of the composite population is described by a concatenation of vector states of each subpopulation:

(6) x~=(x(1)T,x(2)T,,x(r)T)T=(x1(1),,xn(1),x1(2),,xn(2),,x1(r),,xn(r))T,

where x~ is the column vector describing the state of the composite population, x(j) is the column vector describing jth subpopulation in a form (3). Note that the last entries of any 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):

(7) dx~dt=A~x~diag(K~x~)x~.

Here, the composite population projection matrix representing the cell growth, environmentally caused (constant) group death, and group fragmentation is a diagonal block matrix

(8) A~=(A(1)0000A(2)0000A(3)0000A(r)),

where A(i) is the population projection matrix of the life cycle κ(i) extended to size n×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 top-left m(i)×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 K~ is constructed as

(9) K~=(KKKKKKKKKKKKKKKK),

where each block 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 κ(R) and rare ‘invader’ executing different life cycle κ(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 frequency-dependent competition is negligible. The members of the resident population compete predominantly between themselves, so the resident dynamics is effectively a single-species scenario,

(10) dx(R)dt[A(R)diag(Kx(R))]x(R)=0=A(R,R)x(R),

where the vector x(R) represents the composition of the resident population, A(R) is the population projection matrix of the resident, x(R)* is the equilibrium composition, and we introduced the self-invasion population projection matrix A(R,R)=A(R)-diag(Kx(R)*). Since the resident is assumed to be at a stable equilibrium in the absence of invaders, the self-invasion matrix A(R,R) has an eigenvalue that is zero, and the equilibrium population composition 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 self-competition is negligible. Thus, its dynamics is given by

(11) dx(I)dt[A(I)diag(Kx(R))]x(I)=A(I,R)x(I),

where vector x(I) represents the composition of the invader population, and we introduced the invasion matrix

(12) A(I,R)=A(I)diag(Kx(R)).

Unlike the resident dynamics, the dynamics of the invader population is linear – the invasion population projection matrix A(I,R) is independent from the invader population state 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 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 (Kij=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 (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 (bi) and group death rates (di). 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 (κ=1+1), where population is composed only of solitary cells, the dynamics of our model given by Equation (2) reduces to logistic growth,

(13) dx1dt=(b1d1)x1K1,1x12,

where the net growth rate is equal to b1-d1 and the carrying capacity is (b1-d1)/K1,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 x*, where

(14) dxdt=Axdiag(Kx)x=0.
A single life cycle comes to a stationary state.

(A) For a unicellular life cycle 1+1, our model is equivalent to logistic growth. There, the number of solitary cells (red) approaches the carrying capacity from any initial state. (B, C) For multicellular life cycles (2 + 1 depicted), the population also approaches such a stationary state from any initial state. (B) shows the dynamics of solitary cells (red) and bicellular groups (green) with time. (C) shows the trajectories of the population composition, the dot marks the stationary state. (B) and (C) show the same dataset.

Numerical simulations show that a population executing a single life cycle always comes to the same stationary state 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 (κ(1),κ(2)) with the special initial conditions, where one of these life cycles is abundant, while the other one is rare. The life cycle κ(1) can invade from rare into the abundant κ(2) if the largest eigenvalue of the invasion matrix A(1,2) is positive; see Equation (12). Otherwise, the rare life cycle κ(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).

Table 1
Competitive interactions in a pair of life cycles can lead to the dominance of either life cycle, their coexistence, or bistability.
λ(A(1,2))<0
κ(1) cannot invade into κ(2)
λ(A(1,2))>0
κ(1) can invade into κ(2)
λ(A(2,1))<0
κ(2) cannot invade into κ(1)
Bistability between κ(1) and κ(2)Dominance of κ(1)
λ(A(2,1))>0
κ(2) cannot invade into κ(1)
Dominance of Coexistence of κ(1) and κ(2)
Competitive interactions can lead to the bistability between life cycles (A), a dominance of either of them (B,C), or their coexistence (D).

Each panel shows the evolutionary trajectories of populations with various initial abundances of the life cycles 1 + 1 and 2 + 1 (lines) and all final states (dots). In (A), colors indicate to which of two stationary states the trajectory converges. Birth and death rates are b=(1,0.5) and d=(0,0); they favor the unicellular life cycle 1 + 1 in the absence of competition. Simulations performed on each panel only differ in the competition matrix K shown in panels.

Life cycle κ(1) dominates life cycle κ(2) if κ(1) spreads from rare, but κ(2) does not. This is equivalent to the leading eigenvalue of the invasion matrix A(1,2) being positive, while the leading eigenvalue of A(2,1) is negative. Then, independently of the initial conditions, only the life cycle κ(1) survives in the long run (see Figure 3B). The opposite signs result in the dominance of κ(2) over κ(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 A(1,2) and 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).

Evolution in a population with multiple competing life cycles.

(A) The schematics of the life cycles taking part in evolution. Groups colored according to the number of cells. (B) The frequencies of various classes of outcomes. A single life cycle (yellow) is the most common. Coexistence of life cycles (red) and bistability between them (blue) are less common. More complex composite configurations (gray) are also possible, but found rarely.

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).

Table 2
Observation frequencies of different kinds of evolutionary outcomes in a population with multiple competing life cycles.
OutcomeFrequency
Single life cycle0.751
Coexistence ofTwo life cycles0.121
Three life cycles0.044
Four or more life cycles0.003
Multistability betweenTwo life cycles0.065
Three life cycles0.002
Four or more life cycles0
Composite configurationsBistability with a coexisting pair0.007
Other configurations0.007

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 κ(1), κ(2), κ(3). If the life cycle κ(1) is a resident, there are four variants of its stability against invasions: (i) either it is stable against invasions from both κ(2) and κ(3) (then κ(1) is an ESS), or (ii) stable against only κ(2), or (iii) stable against only κ(3), or (iv) both κ(2) and κ(3) can invade. Similar four variants exist for the two other life cycles. As a result, for the whole triplet, there are 43=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 higher-order life cycle dominates (always invade) lower-order 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).

Constrained triplets demonstrate fewer patterns of pairwise invasion.

(A) For a combination of three life cycles, there are 223=64 pairwise invasion patterns possible. For the triplet κ(1)=5+3,κ(2)=4+4,κ(3)=4+2+2, the rates of all processes (bi,di,Ki,j) were randomly sampled from an exponential distribution with unit rate parameter. Then the pairwise invasion pattern was identified. All 64 possible patterns were observed in these simulations. (B) In a similar investigation for the triplet κ(C1)=2+2,κ(C2)=4+4,κ(M)=4+2+2, in which κ(C1) and κ(C2) constrain κ(M), only eight patterns were found; see the main text for a discussion.

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 κ(M) can be represented as a combination of fragmentation rules of two other (constraining) life cycles κ(C1) and κ(C2). The simplest example is the triplet κ(C1)=2+1, κ(C2)=1+1, κ(M)=1+1+1, where the fragmentation of a three-celled group into three solitary cells (31+1+1) can be presented as a combination of the detachment of a single cell (32+1) immediately followed by the dissolution of the two-cell group (21+1). A lot of constrained triplets can be constructed, for example, in our illustrations we use the triplet κ(C1)=2+2,κ(C2)=4+4,κ(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 κ(M) was found to be always between growth rates of the constraining life cycles κ(C1) and κ(C2). It follows that in the present model with competition (Kij0), each of two constraining life cycles (κ(C1) and κ(C2)) must be either stable against two others or unstable against both. The constrained life cycle κ(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 222=8 (see Figure 5B). Among these eight patterns, two feature hierarchical dominance, where the life cycles can be ordered in a way that a higher-order life cycle dominates the lower-order life cycles (with the constrained life cycle κ(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 κ(C1) and κ(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 non-hierarchical 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.

There are only eight patterns of pairwise invasion in a constrained triplet of life cycles.

Four patterns are shown here and four more are symmetric to them. Blue lines are population dynamics trajectories with different initial composition of the population. Red points are final states. Black arrows show directions of invasion from rare. (A) In a hierarchical dominance, one of the constraining life cycles outcompetes the two others from any initial condition. (B) In a bistability situation, each of the two constraining life cycles is able to outcompete the third life cycle. Which of them will survive depends on the initial conditions. Trajectories converging to different stationary states are highlighted with different colors. (C) In a coexistence, all three life cycles survive to a stationary state. (D) In a non-hierarchical dominance, one of the constraining life cycles outcompetes another; however, in the abundance of the constrained life cycle, the dominance is reversed. The parameters used to produce these figures are presented in Appendix 5.

The fundamental feature of constrained triplets is that any constrained life cycle (κ(M)) can always be invaded by exactly one constraining life cycles (κ(C1) and κ(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 (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

(15) Kij=kj.

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 no-interactions model (Kij=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 no-interaction 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 (A~) and competition (K~) matrices. If the competition matrix K is a killer kernel, the composite competition matrix 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 A~ corresponding to its leading eigenvalue. The composite population projection matrix 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 A~ is the largest among the eigenvalues of all population projection matrices comprising A~. Additionally, the corresponding eigenvector has nonzero components only at the positions in 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 (Kij=kj) 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

(16) Kij=ki.

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*=ixi* is equal to the growth rate of this life cycle in the no-interaction model (Kij=0) with modified cell division rates bi=bi/ki and group death rates di=di/ki (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 κ(M), its invasion matrix (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,

(17) di=di+b1d1K1,1Ki,1.

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 (A(1,M)) has size 1×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

(18) b1d1i=1m(M)K1,ixi(M)<0,

where m(M) is the maximal group size in the life cycle κ(M), and xi(M)* is the composition of the population when κ(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 eco-evolutionary 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, frequency-dependent 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, κ(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 κ(C1)=2+2 and κ(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 eco-evolutionary 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 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 A is not diagonal. Of course, if the population projection matrix A can be diagonalized, we can perform a linear transformation xCy (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 (-xijKijxj, see Equation (1)) and will become a general second-order term instead (-jkKijxjxk). 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 fine-tuned to mathematically beautiful outcomes. What might empirical competition matrices in a stage-structured 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 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 (Ki,11) and are the most susceptible to oxygen denial (K1,i1). In opposite limit, an established mat of millions of cells just drives everything else to extinction (Ki,j11). For arbitrary competitors size, the terms Kij should increase with the size of an opponent group (Kij<Ki,j+1) but decrease with the size of the focal group (Kij>Ki+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 κ(R), in which more than two offspring groups are produced as a result of fragmentation. The initial dynamics of any life cycle κ(I) invading from rare can be described by a linear model with death rates modified as

(19) didi+jKi,jxj(R).

The invasion is successful if the leading eigenvalue of the corresponding population projection matrix 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 (bi) and group death rates (di), 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 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 self-invasion 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 (Ki,j=0) and under the killer kernel (Ki,j=kj).

First, we consider a linear model, where the population dynamics is governed by the population projection matrix A:

(20) dxdt=Ax.

A number of useful properties of this dynamics comes from the Perron–Frobenius theorem for irreducible non-negative matrices. Here, non-negative matrix means a matrix with all elements greater or equal to zero. Matrix M is irreducible if its representation by a directed graph, where node i has an edge to node j only if Mij>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 A describes the dynamics of the system, the irreducibility of A means that group of any given size i can reach any other size j through growth and fragmentation.

However, the population projection matrix A itself is not non-negative as its main diagonal is negative (Aii=-ibi-di). It is not always irreducible as well: if the fragmentation mode produces only multicellular offspring (e.g., κ=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 A=A+Imax(ibi+di) is non-negative. The matrix A has the same eigenvectors as A, and its eigenvalues (λ) relate to eigenvalues of the population projection matrix (λ) as λ=λ+max(ibi+di). 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 long-term 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 non-negative and irreducible, hence, the Perron–Frobenius theorem applies. From this theorem, it immediately follows that for arbitrary birth rates (bi), death rates (di), and the life cycle executed (κ), the population projection matrix has a simple leading eigenvalue λ, its eigenspace is one-dimensional, 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)

(21) x(t)=weλt,

where λ is the leading eigenvalue of the matrix A, and W* is the corresponding eigenvector,

(22) Aw=λw.

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

(23) dxdt=Axdiag(Kx)x=Axdiag(jkjxj1)x=(AjkjxjI)x,

where 1 is the vector of ones, and I is the identity matrix (diag(1)=I). It can be shown by a direct calculation that the vector

(24) x=λikiwiw

is the stationary state of the dynamics (23):

(25) dxdt=Axjkjxjx=λikiwiAw(λikiwi)2jkjwjw=λ2ikiwiwλ2ikiwiw=0.

Note that while any other eigenvector W of the population projection matrix A can be used in Equation (25), only the leading eigenvector 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 x* under the killer kernel is proportional to the vector describing the stationary distribution 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*=ixi*) under the victim kernel (Ki,j=ki) is mathematically equivalent to the task of finding the population growth rate in the linear model (Ki,j=0) with modified cell birth rates (bibi/ki) and group death rates (didi/ki).

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

(26) 0=det(A-λI)=|-b1-d1-λ00mbmπ1(κ)b1-2b2-d2-λ0mbmπ2(κ)02b2-3b3-d3-λmbmπ3(κ)000mbmπm(κ)-mbm-dm-λ|.

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

(27) dxdt=Axdiag(Kx)x=Axdiag(kjxj)x=(ANdiag(k))x=0,

where k is a vector constructed from any column of the competition matrix K (they are all identical), and N*=ixi* 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):

(28) x=0,ordet(ANdiag(k))=0

Limiting ourselves to scenarios where the stationary state x* is not an empty population, we can conclude that the population at the equilibrium satisfies

(29) 0=det(A-N*diag(k))=|-b1-d1-k1N*00mbmπ1(κ)b1-2b2-d2-k2N*0mbmπ2(κ)02b2-3b3-d3-k3N*mbmπ3(κ)000mbmπm(κ)-mbm-dm-kmN*|=k1k2k3km|-b1k1-d1k1-N*00mbmkmπ1(κ)b1k1-2b2k2-d2k2-N*0mbmkmπ2(κ)02b2k2-3b3k3-d3k3-N*mbmkmπ3(κ)000mbmkmπm(κ)-mbmkm-dmkm-N*|,

where in the last step, we divided the ith column by ki for all i. Comparing the determinants in Equations (26) and (29), we find that they are identical after the substitution

(30) λN,bibiki,didiki.

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 (Ki,j=ki), 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 (Ki,j=0) with modified cell division (bibi/ki) and group death (didi/ki) rates.

If the competition matrix K constitutes a victim kernel, then the composite competition matrix 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 (N~*=p=1ri=1nxi(p)) is determined by

(31) 0=det(A~N~diag(k~))=p=1rdet(A(p)N~diag(k)),

where k~ is the vector made from a single column of the composite competition matrix K~ (it is a concatenation of r vectors k). In the second step of Equation (31), we used the property of block-diagonal matrices: the determinant of a block-diagonal matrix is the product of determinants of its blocks. Equation (31) is satisfied if one of the multipliers is equal to zero. The pth multiplier becomes zero when the composite population reaches a size equal to the carrying capacity of the pth life cycle; cf. Equation (28). At that moment, the pth 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 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 (κ(R)=1+1), its steady state is given by the solution of

(32) (b1d1)x1K1,1x12=0,

equal to

(33) x1=b1d1K1,1.

Then, for an arbitrary invader multicellular life cycle κ(M), the invasion matrix is given by

(34) A(M,R)=A(M)diag(Kx(R))=A(M)diag(Ki,1x1)=A(M)b1d1K1,1diag(Ki,1).

This is equivalent to the linear growth of the invader life cycle in an environment with modified death rates

(35) di=di+b1d1K1,1Ki,1.

If the invader is unicellular (κ(I)=1+1), then the invasion matrix is effectively a 1×1 matrix since the invader contains only isolated cells. Even if A(I,M) formally has a larger size, it is a block matrix, with the element A1,1(I,M) being a block, with a value that is equal to the growth rate of the unicellular life cycle. This value is

(36) A1,1(I,M)=A1,1(I)diag(Kx(M))1,1=b1d1i=1m(M)K1,ixi(M),

where xi(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 A1,1(I,M)<0.

Parameters of calculation used in figures

Figure 2

In Figure 2A, we used b1=1,d1=0,K11=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

K=(10.20.20.5).

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 x1+1=(s1), x2+1=(s2,0), where s1,s2{0.1,0.2,0.3,,1.0}. The dynamics shown in the four panels differ by the competition matrix used:

  • Panel A

    K=(1111)
  • Panel B

    K=(3311)
  • Panel C

    K=(10.10.11)
  • Panel D

    K=(110.60.4)

Figure 6

In Figure 6, we considered the life cycles κ(C1)=2+2, κ(C2)=4+4, κ(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:

(37) x~t=0=(s(C1)x(C1)T,s(C2)x(C2)T,s(M)x(M)T)T,

where the vectors x(C1)T, x(C2)T, x(M)T are equilibrium population states of corresponding life cycles grown in isolation, computed according to Equation (14).

  • In panel A (hierarchic dominance), we used bi=1.0, di=0 and the competition matrix Kij=0.1:

­

(38) KHier. dom.=(0.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.11.11.11.11.11.11.11.1),

or equivalently

(39) Kij={1.1,i=70.1,otherwise
  • In panel B (bistability), we used bi=1 and di=0, while the competition matrix was

­

(40) KBi-stability=(0.10.10.10.10.10.10.10.10.10.10.10.10.1100.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.10.1100.10.10.10.10.1),

or equivalently

(41) Kij={10,(i,j)=(2,7) or (7,2)0.1,otherwise
  • In panel C (coexistence), we used bi=1 and di=0, while the competition matrix was

­

(42) KCoex.=(0.10.10.10.10.10.100.10.10.10.10.100.10.10.10.10.100.10.10.10.10.100.10.10.10.10.100.10.10.10.10.100.10.10.10.10.100.10.10.10.10.10.1),

or equivalently

(43) Kij={0,i+j=80.1,otherwise
  • In panel D (non-hierarchical dominance), we used di=0,

­

(44) bNonh. dom.=(6.48076,2.79693,2.3088,5.34057,1.0,1.62478,1.32615)

and

(45) KNonh. dom.=(0.124910.164530.409720.139810.034960.793640.060970.364770.108590.090990.703910.017140.533540.497780.164320.5828510.019180.012680.080710.112080.135190.567710.118790.026010.089050.11720.146610.012510.313530.026390.074330.053120.228770.148410.062580.438330.306790.33230.010140.096370.247510.199450.0220.000870.24690.097330.082470.37168).

Data availability

The code and data are available from https://github.com/yuriypichugin/Eco-evolutionary-dynamics-life-cycles, (copy archived at swh:1:rev:b7122aeec7b11953867e6b5e588701e5a602276d).

References

  1. Book
    1. Hoffman K.
    (1971)
    Linear algebra
    Englewood Cliffs, NJ: Prentice-Hall.
    1. Levin BR
    (1988) Frequency-dependent selection in bacterial populations
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 319:459–472.
    https://doi.org/10.1098/rstb.1988.0059
  2. Book
    1. Maynard Smith J
    2. Szathmáry E
    (1995)
    The Major Transitions in Evolution
    Oxford: W. H. Freeman.

Decision letter

  1. Babak Momeni
    Reviewing Editor; Boston College, United States
  2. Aleksandra M Walczak
    Senior Editor; CNRS LPENS, France
  3. Denis Tverskoi
    Reviewer

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 "Eco-evolutionary 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 double-check 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 quasi-species 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 life-cycle 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 rock-paper-scissor 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 11-12 describing these dynamics.

Page 11, lines 245-246: “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 526-529. 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 551-552. 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 450-451. “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 431-442) 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 443-505).

https://doi.org/10.7554/eLife.78822.sa1

Author 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 double-check 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 quasi-species 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 quasi-species 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 quasi-species 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 life-cycle 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 rock-paper-scissor 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 non-transitive 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 11-12 describing these dynamics.

Agreed, it also makes the notation in Equations 11 and 12 consistent with the rest of the paper.

Page 11, lines 245-246: “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 526-529. 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 positive-only elements. This follows from Perron-Frobenius 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 Perron-Frobenius 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 551-552. 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. bi-stability 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 450-451. “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 Lotka-Volterra 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 431-442) 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 443-505).

Thank you, we have now rewritten the discussion along the suggestions of you and the other reviewers.

https://doi.org/10.7554/eLife.78822.sa2

Article and author information

Author details

  1. Vanessa Ress

    1. Max Planck Institute for Evolutionary Biology, Plön, Germany
    2. Hamburg Center for Health Economics, University of Hamburg, Hamburg, Germany
    Contribution
    Software, Formal analysis, Investigation, Methodology, Writing - original draft
    Competing interests
    No competing interests declared
  2. Arne Traulsen

    Max Planck Institute for Evolutionary Biology, Plön, Germany
    Contribution
    Formal analysis, Methodology, Writing - original draft, Project administration
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0669-5267
  3. Yuriy Pichugin

    1. Max Planck Institute for Evolutionary Biology, Plön, Germany
    2. Department of Ecology and Evolutionary Biology, Princeton University, Princeton, United States
    Contribution
    Conceptualization, Software, Formal analysis, Supervision, Visualization, Methodology, Writing - original draft
    For correspondence
    pichugin@princeton.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3078-2499

Funding

No external funding was received for this work.

Senior Editor

  1. Aleksandra M Walczak, CNRS LPENS, France

Reviewing Editor

  1. Babak Momeni, Boston College, United States

Reviewer

  1. Denis Tverskoi

Publication history

  1. Preprint posted: March 16, 2022 (view preprint)
  2. Received: March 21, 2022
  3. Accepted: August 9, 2022
  4. 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

  • 248
    Page views
  • 78
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

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)

  1. Vanessa Ress
  2. Arne Traulsen
  3. Yuriy Pichugin
(2022)
Eco-evolutionary dynamics of clonal multicellular life cycles
eLife 11:e78822.
https://doi.org/10.7554/eLife.78822

Further reading

    1. Cancer Biology
    2. Computational and Systems Biology
    Erika K Ramos, Chia-Feng Tsai ... Huiping Liu
    Research Article

    Tumor-initiating cells with reprogramming plasticity or stem-progenitor cell properties (stemness) are thought to be essential for cancer development and metastatic regeneration in many cancers; however, elucidation of the underlying molecular network and pathways remains demanding. Combining machine learning and experimental investigation, here we report CD81, a tetraspanin transmembrane protein known to be enriched in extracellular vesicles (EVs), as a newly identified driver of breast cancer stemness and metastasis. Using protein structure modeling and interface prediction-guided mutagenesis, we demonstrate that membrane CD81 interacts with CD44 through their extracellular regions in promoting tumor cell cluster formation and lung metastasis of triple negative breast cancer (TNBC) in human and mouse models. In-depth global and phosphoproteomic analyses of tumor cells deficient with CD81 or CD44 unveils endocytosis-related pathway alterations, leading to further identification of a quality-keeping role of CD44 and CD81 in EV secretion as well as in EV-associated stemness-promoting function. CD81 is co-expressed along with CD44 in human circulating tumor cells (CTCs) and enriched in clustered CTCs that promote cancer stemness and metastasis, supporting the clinical significance of CD81 in association with patient outcomes. Our study highlights machine learning as a powerful tool in facilitating the molecular understanding of new molecular targets in regulating stemness and metastasis of TNBC.

    1. Cell Biology
    2. Computational and Systems Biology
    Julie Paxman, Zhen Zhou ... Nan Hao
    Research Article

    Chromatin instability and protein homeostasis (proteostasis) stress are two well-established hallmarks of aging, which have been considered largely independent of each other. Using microfluidics and single-cell imaging approaches, we observed that, during the replicative aging of S. cerevisiae, a challenge to proteostasis occurs specifically in the fraction of cells with decreased stability within the ribosomal DNA (rDNA). A screen of 170 yeast RNA-binding proteins identified ribosomal RNA (rRNA)-binding proteins as the most enriched group that aggregate upon a decrease in rDNA stability induced by inhibition of a conserved lysine deacetylase Sir2. Further, loss of rDNA stability induces age-dependent aggregation of rRNA-binding proteins through aberrant overproduction of rRNAs. These aggregates contribute to age-induced proteostasis decline and limit cellular lifespan. Our findings reveal a mechanism underlying the interconnection between chromatin instability and proteostasis stress and highlight the importance of cell-to-cell variability in aging processes.