Electronic Supplementary Material from “ Competing metabolic strategies in a multilevel selection model ”

This electronic supplementary material provides detailed analytical calculations for the equilibria and stability analysis for well-mixed populations. Additionally, it provides the parameters summary with their description and units.


The Jacobian matrix
Here we calculate the eigenvalues of the Jacobian matrix of the discrete model. If the absolute values of the eigenvalues are smaller than one then the solution in Section 3.2.1 is stable. If we write the set of equations (10) of our two-variable model as n D (t+1) = f D (n D , n C ) and n C (t+1) = f C (n D , n C ), the Jacobian matrix is defined as The Jacobian matrix must be evaluated at the equilibrium of interest, i.e. n C = 0 andn D = − α AT P D S ln(1−ν/a D ) . After some algebra, the coefficients of the matrix at the equilibrium of interest are determined where ∆ AT P = α AT P D α AT P C and = A D A C . The above matrix is upper diagonal which means that the eigenvalues of the Jacobian are simply the diagonal coefficients. Because ν < a D for the referred equilibrium, the element J 11 is always smaller than one and positive regardless the values of ν and a D . Note that the logarithmic term is always negative. Therefore, the stability of the solution is uniquely settled by the element J 22 . Of particular interest is the point at which the equilibrium becomes unstable which occurs when J 22 equals one. Making J 22 equal to one a relation between Γ AT P and ∆ AT P is obtained.
The same analysis is done in Section 3.2.2, but instead the solutionn D = 0

Coexistence solution
The system of equations (10) also admits a solution where the two strains coexist. By rearranging the set of equations we verify that the coexistence solution obeys the following set of equations: showing that there is an indeterminacy, and the coexistence solution only exists when ∆ AT P = 1 ln(1−ν/a D ) ln(1−ν/a C ) . Actually Eq. (SI 1) describes a family of solutions. Of course, the range of values at which n D and n C are meaningful are n C ∈ [0, N ] and n D ∈ [0, N ], where N = − α AT P C S ln(1−ν/a C ) with the constraint that n D + n C = − α AT P C S ln(1−ν/a C ) . Note that when n C = 0 we recover the solution n D = − α AT P D S ln(1−ν/a D ) , whereas if n D = 0 we recover the solution n C = − α AT P C S ln(1−ν/a C ) . By determining the Jacobian of the system at equilibrium and hence calculating its corresponding eigenvalues we check whether the coexistence solution is stable or not. With the help of mathematica, we performed a numerical computation of the eigenvalues for different set of parameters by varying all the range of acceptable values of n D (since n C is written as a function of n D ). As expected in this situation, one of the eigenvalues is equal to one, which means that any disturbance in the same direction as the line drives the system to the new values of (n D , n C ). This is called marginal stability as found in a simple predator-prey model studied by Lotka and Volterra (Parker and Kamenev 2010). The stability of the family of solutions will be dictated by the second eigenvalue. For the set of parameter values used in our simulations the second eigenvalue is slightly smaller than one, meaning that any transverse disturbance is damped by the system, thus warranting stability of the family of solution in Eq. (SI 1). It is important to highlight that the coexistence solution only holds for ∆ AT P = 1 ln(1−ν/a D ) ln(1−ν/a C ) , which describes a line in the diagram ∆ AT P versus Γ AT P .
In addition, from Eq. (10) we also notice thatn D = 0 andn C = 0 are also equilibrium solutions, where both types go extinct. The solution is only stable when ν > a D and ν > a C .
The coexistence of the two strains can not be observed in the simulations as the family of solutions embodies the absorbing states n C = 0 and n D = 0. As stochasticity plays a key role in finite populations, in a finite time the system will eventually reach one of these boundaries. Since we always start from n C = 1, the system will probably evolve to n C = 0 in a short time.

Population size at equilibrium
Here we show how the stationary population size of individuals type D (before introducing a single strain C), N st , depends on the set of parameters. From this we gain some insight about the strength of stochasticity. As intuitively expected, a linear relation between population size and the influx of resource into the system is observed (data not shown). The linearity is not followed when studying the dependence of the population size on the other parameters.
In Figure S1 we , respectively. The range of values of Γ AT P and ∆ AT P are, as argued before, those that are meaningful to the problem here addressed. Please note that a logarithmic scale is adopted for the color gradation in Fig. S1. The range of ∆ AT P = α AT P D /α AT P C in which population size is shown for the well-mixed populations is much smaller than the range displayed for structured population. Outside of this range the population already falls in the regime where the efficient strain is counter-selected -as can be seen in Fig. 2 (Main Text). The right panel of Fig. S1 shows the equilibrium solution for the discrete-time model of a wellmixed population, given byn = − , as derived in Section 3.1.1. The grey region denotes the extinction of the population, where the above solution becomes unstable, while the solutionn = 0 becomes stable. The plot evinces a very good agreement between the theoretical prediction and simulations. The simulation results for a well-mixed population display a larger grey area in comparison to the theoretical prediction, which owes to stochastic effects, which are not captured by the time-discrete model. In the right panel, the population size is too small -around five or less individuals -being prone to extinction.
Especially, in structured populations the stationary population size N st of a pure population of strain D can vary from less than 100 individuals (small Γ AT P and ∆ AT P ) to nearly 5, 000 individuals (large Γ AT P and ∆ AT P ). The wide range of the population size at equilibrium explains why we should adopt a relative measurement of fixation probability since the strength of drift, namely 1/N st , is variable under the change of the parameter values. . The other parameter values are influx rate of resource S = 25, death rate ν = 0.01, group carrying capacity P max = 10, internal energy for one individual to split to two individuals E max = 10, and A D = 10. The simulation data points correspond to 40 distinct populations and for each population 10, 000 independent runs were performed. The black lines denote isoclines, along which the population size is a constant (corresponding values indicated). The grey region denotes a region where the defector population is not sustainable and goes to extinction (before inserting the cooperator).

Parameters summary
For the sake of clarity we summarise the parameters of the model, including a short description and units. We use as basic units the time step, resource and energy.

Parameter Interpretation
Units ν death rate time-step −1 S total available resource resource A C rate of resource capture of C resource/time-step A D rate of resource capture of D resource/time-step A AT P C maximum rate of resource processing of C energy/time-step A AT P D maximum rate of resource processing of D energy/time-step α AT P D efficiency of D in resource processing resource −1 α AT P C efficiency of C in resource processing resource −1 E max energy threshold for cell split energy P max threshold size for group split dimensionless m migration rate (probability of each cell to migrate to a random group per time-step) time-step −1 Also we include information about some derived quantities that are important to analyse the results, despite not being parameters of the model.

Quantity Definition Interpretation
Units S G S/N groups resource available per group resource J S C A C resource uptake of C resource/time-step J S D A D resource uptake of D resource/time-step J AT P C (see eq. 4) rate of resource processing of C energy/time-step J AT P D (see eq. 4) rate of resource processing of D energy/time-step A D /A C fraction of the resource acquired by D in a pairwise interaction dimensionless ∆ AT P α D /α C relative efficiency of D and C dimensionless Γ AT P A AT P D /A AT P C ratio between the maximum rates of resource processing of D and C dimensionless r J AT P D /J AT P C relative advantage of D over C in a pairwise interaction dimensionless a D A AT P D /E max maximum growth rate of the D time-step −1 a C A AT P C /E max maximum growth rate of the C time-step −1 ∆E ij J AT P ∆t internal energy increase energy