The evolution of biogeochemical recycling by persistence-based selection
Model descriptionThe model involves a discrete time, discrete valued stochastic Markov process. Model variables and parameters are given in Tables 1 and 2 respectively. Both time and the number of individuals of each type are constrained to be integer valued. Death and reproductive mutation are stochastic processes derived from sampling from binomial distributions given by the relevant probabilities. All ensemble results give the 100-replicate average for the parameter choices in question.Growth of individuals from species ({S}_{1}) and ({S}_{2}) is proportional to the bio-available level of environmental substances ({R}_{1}) and ({R}_{2}) respectively. At time (t) (where time is in units of biological generations) the change in the number ({N}_{q,j}) of individuals of genotype (j) (non-producer, producer, plastic) within species (q) (({S}_{1}) or ({S}_{2})) can be written as a function of the state of the variables at the previous time-step:$${N}_{q,j}left(t+1right)=left(left({N}_{q,j}left(tright)-{rho }_{q,j}left(tright)right)cdot {G}_{q,j}left(tright)-{{{{{{rm{{Upsilon }}}}}}}}_{q,j}left(tright)-{delta }_{q,j}left(tright)+{{{{{{rm{{Upsilon }}}}}}}}_{q,xne j}left(tright)right)cdot left(1-frac{{S}_{q}left(tright)}{K}right)$$
(1)
The leftmost bracket on the right-hand side represents the number of individuals escaping starvation (death due to insufficient environmental substance) at the previous time-step and ({G}_{q,j}left(tright)) is the per capita reproductive growth rate. ({{{{{{rm{{Upsilon }}}}}}}}_{q,j}left(tright)) gives the number of mutant offspring individuals produced during reproduction from parent individuals of genotype (j). ({{{{{{rm{{Upsilon }}}}}}}}_{q,xne j}left(tright)) represents the number of (j) genotype individuals derived from mutation in parent individuals of other genotypes. ({delta }_{q,j}left(tright)) is the number of individuals of genotype (j) lost to random death, and the rightmost bracket relates the total number ({S}_{q}left(tright)) of individuals of species (q) to carrying capacity (K), which represents limitation of growth by any factor other than the relevant environmental substance, e.g. space. (The steady state population size in all simulations shown is below (K) and limited by the environmental substance influx. The carrying capacity is included in the model for computational reasons and as a “crash preventer” but has no qualitative effect on the results).The total number of individuals ({S}_{q}left(tright)) in species (q) is the sum of the number of individuals of each genotype (producer, non-producer and plastic, as discussed in the main text):$${S}_{q}left(tright)=mathop{sum }limits_{j=1}^{{j}_{{total}}}{N}_{q,j}left(tright)={N}_{q,{prod}}left(tright)+{N}_{q,{non}-{prod}}left(tright)+{N}_{q,{plast}}left(tright)$$
(2)
The genotype-specific reproductive growth rate ({G}_{q,j}left(tright)) (again for genotype (j) within species (q), time (t)), gives the number of offspring individuals produced per parent individual, per time-step. Growth rate is an increasing function of the bio-available level of environmental substance ({R}_{q,{BIOAVAIL}{ABLE}}) (the subscript (q) being identical because species ({S}_{1}) and ({S}_{2}) assimilate substances ({R}_{1}) and ({R}_{2}) respectively). Growth rate also includes a substance-to-biomass conversion efficiency parameter ({f}_{{conv}}) and a genotype-specific per capita term ({G}_{q,j,{PR}}) (number of offspring per parent, per unit environmental substance assimilated, per unit time). In the absence of growth-limitation by environmental substance levels, growth rate is capped at a genotype-specific maximum ({G}_{q,{jMAX}}):$${G}_{q,j}left(tright)={MIN}[{G}_{q,j,{PR}}cdot {R}_{q,{BIOAVAILABLE}}(t)cdot {f}_{{conv}},{G}_{q,{jMAX}}]$$
(3)
$${G}_{q,{non}-{prod},{PR}}={G}_{0}$$
(4)
$${G}_{q,{non}-{prod},{MAX}}={G}_{0}cdot {R}_{{assimMAX}}$$
(5)
({G}_{0}) is the baseline number of offspring, per parent, per unit substance assimilated. ({R}_{{assimMAX}}) is a universal maximum potential number of units of environmental substance that can be assimilated by a single individual per time-step (i.e. representing basic physiological constraints on growth). The producer genotype incurs a per capita reproductive growth rate cost ({kappa }_{{prod}}) relative to the non-producer:$${G}_{q,{prod},{PR}}={G}_{0}cdot (1-{kappa }_{{prod}})$$
(6)
$${G}_{q,{prod},{MAX}}={G}_{0}cdot (1-{kappa }_{{prod}})cdot {R}_{{assimMAX}}$$
(7)
This growth rate formulation is therefore a highly simplified linearization of the Michaelis-Menten kinetics normally used in models of resource and nutrient assimilation.The plastic genotype switches phenotype depending upon the level of environmental substance relative to a fixed threshold ({{R}_{q,{BIOAVAILABLE}}}_{{crit}}), in effect becoming a second non-producer genotype below this threshold and a second producer genotype above it:$${IF}[{R}_{q,{BIOAVAILABLE}}(t)ge {{R}_{q,{BIOAVAILABLE}}}_{{crit}}],{G}_{q,{plast}}left(tright)={G}_{q,{prod}}left(tright)$$$${ELSEIF}[{R}_{q,{BIOAVAILA}{BLE}}left(tright) , < , {{R}_{q,{BIOAVAILABLE}}}_{{crit}}],{G}_{q,{plast}}left(tright)={G}_{q,{non}-{prod}}left(tright)$$
(8)
There is no spatial structure whatsoever, thus access to environmental substance is uniform across individuals. The bioavailable quantity of each environmental substance is simply the total amount ({R}_{q,{NET}}(t)) divided by the total number of individuals assimilating it:$${R}_{q,{BIOAVAILABLE}}left(tright)=frac{{R}_{q,{NET}}(t)}{{S}_{q}left(tright)}$$
(9)
We allow the per capita reproductive growth rate to fall below ({G}_{q,j}left(tright)=1), which, if interpreted deterministically at the individual level would correspond to an individual failing to sustain its biomass to the next time-step and thus dying. However, a population-level average ({G}_{q,j}left(tright) , < , 1) is interpretable in terms of a thinning factor that maps between discretized individuals and continuously distributed environmental substance. Thus, a thinning factor of (left(1-{G}_{q,j}left(tright)right)) is used to calculate the total number of individuals dying of starvation ({rho }_{q,j}) (again genotype (j), species (q)). This represents pre-reproduction deaths, corresponding to the difference between the actual population size and the population size that the environmental substance pool is capable of supporting. ({rho }_{q,j}left(tright)) is constrained to be an integer and is zero for ({G}_{q,j}left(tright) , > , 1):$${rho }_{q,j}left(tright)={N}_{q,j}left(tright)cdot {MAX}left[0,left(1-{G}_{q,j}left(tright)right)right]$$
(10)
A subset of offspring are a different genotype from their parent via mutation. For parent genotype (j), the number ({{{{{{rm{{Upsilon }}}}}}}}_{q,j}left(tright)) of mutant offspring with genotype (ne j) is calculated using baseline mutation probability per reproductive event ({mu }_{0}), with the total number of new individuals produced by the parent individuals surviving starvation ({G}_{q,j}left(tright)cdot left({N}_{q,j}left(tright)-{rho }_{q,j}left(tright)right)). The total number of mutant offspring ({{{{{{rm{{Upsilon }}}}}}}}_{q,j}left(tright)) is thus a binomially distributed random variable with success probability ({mu }_{0}) and number of trials ({G}_{q,j}left(tright)cdot left({N}_{q,j}left(tright)-{rho }_{q,j}left(tright)right)). The expected value (Eleft[{{{{{{rm{{Upsilon }}}}}}}}_{q,j}left(tright)right]) is the product of these two numbers:$$ {{{{{{rm{{Upsilon }}}}}}}}_{q,j}left(tright) sim Bleft({G}_{q,j}left(tright)cdot left({N}_{q,j}left(tright)-{rho }_{q,j}left(tright)right),{mu }_{0}right), \ Eleft[{{{{{{rm{{Upsilon }}}}}}}}_{q,j}left(tright)right]={G}_{q,j}left(tright)cdot left({N}_{q,j}left(tright)-{rho }_{q,j}left(tright)right)cdot {mu }_{0}$$
(11)
Mutation to genotype (j) from the other genotypes is calculated in exactly the same way using the number and reproductive growth rates of the relevant (other) genotypes. Any particular mutant offspring is randomly allocated to one of the other genotypes with equal probability ({p}_{kto j}=frac{1}{{j}_{{total}}-1}=0.5) (where ({j}_{{total}}=3) is the total number of genotypes per species). The expected number of offspring with genotype (j) produced by mutation within parent offspring of other genotypes (kne j) is therefore:$$E[{{{{{{rm{{Upsilon }}}}}}}}_{q,x , ne , j}left(tright)]={left(mathop{sum }limits_{k=1}^{{k}_{{to}{tal}}}{{{{{{rm{{Upsilon }}}}}}}}_{q,k}left(tright)right)}_{kne j}cdot frac{1}{{j}_{{total}}-1}$$
(12)
Independently of reproduction and assimilation of environmental substance, any given individual has a probability ({delta }_{0}) at each time point of death due to stochastic factors. The genotype/species specific number of such deaths is again a random sample from a binomial distribution, with success probability ({delta }_{0}):$${delta }_{q,j}left(tright) sim Bleft({N}_{q,j}left(tright),{delta }_{0}right),Eleft[{delta }_{q,j}left(tright)right]={N}_{q,j}left(tright)cdot {delta }_{0}$$
(13)
The net quantity of growth-limiting environmental substance at each time-step is given by the difference between total biotic assimilation ({A}_{{R}_{q}}) and the production ({P}_{{R}_{q}}) and abiotic input ({varphi }_{{R}_{q}}) fluxes:$${R}_{q,{NET}}(t+1)={varphi }_{{R}_{q}}(t)+{P}_{{R}_{q}}(t)-{A}_{{R}_{q}}(t)$$
(14)
The abiotic net influx is the sum of two fluxes. First, an input term that is the product of a baseline scaling factor ({{varphi }_{0}}_{{R}_{q}}) and a model forcing (frac{partial {t}_{{geo}}}{partial {t}_{{bio}}}) representing the mapping between abiotic-geological and biotic-evolutionary timescales. In practice (frac{partial {t}_{{geo}}}{partial {t}_{{bio}}}(t)) was set to either (1) or (0) or (in fluctuation runs) a time-dependent switching between the two. (More sophisticated implementations of (frac{partial {t}_{{geo}}}{partial {t}_{{bio}}}(t)), e.g. sinusoidal oscillations and stochastic time dependence, were attempted but made little qualitative difference to the results). Second, an abiotic removal term that scales linearly with the quantity of environmental substance:$${varphi }_{{R}_{q}}(t)={{varphi }_{0}}_{{R}_{q}}cdot frac{partial {t}_{{geo}}}{partial {t}_{{bio}}}left(tright)-frac{{R}_{q,{NET}(t)}}{{R}_{q,{NET}0}}$$
(15)
where ({R}_{q,{NET}0}) is a normalization factor representing the sensitivity of the abiotic efflux to the influx. In the absence of any biota and for (frac{partial {t}_{{geo}}}{partial {t}_{{bio}}}=1) the steady state environmental substance level is immediately given by (15) as ({R}_{q,{NET}(t)}={R}_{q,{NET}0}cdot {{varphi }_{0}}_{{R}_{q}}), thus the numerical value of ({R}_{q,{NET}0}) corresponds to the abiotic steady state residence time.Total biotic assimilation ({A}_{R}) of each environmental substance is given by:$${A}_{{R}_{q}}left(tright)=mathop{sum }limits_{k=1}^{{j}_{{total}}}frac{{G}_{q,k}left(tright)cdot left({N}_{q,k}left(tright)-{rho }_{q,k}left(tright)right)}{{G}_{q,{kPR}}}$$
(16)
The numerator gives the total number of individuals produced as a result of biological assimilation of environmental substance ({R}_{q}) and the denominator is the genotype specific number of individuals produced per unit substance assimilated, dividing through by which therefore converts to total units of substance assimilated by the population as a whole.Net biotic ({P}_{{R}_{q}{NET}}) production of substance ({R}_{q}) by the producer genotype in the other species (p) is calculated equivalently, via the product of the per capita production rate ({P}_{q,{prod}}) and the total number of reproducing individuals:$${P}_{q,{prod}}left(tright)=frac{{MIN}[{G}_{q,{prod}}left(tright)cdot {f}_{{convprod}},{G}_{q,{prodMAX}}]}{{G}_{p,{PRODUCER;PR}}}$$
(17)
$${P}_{{R}_{q}{NET}}left(tright)={P}_{q,{prod}}left(tright)cdot left({N}_{p,{PRODUCER}}left(tright)-{rho }_{p,P{RODUCER}}left(tright)right)$$
(18)
Where ({f}_{{conv},{PROD}}) is the per capita efficiency by which producers convert the environmental substance that they assimilate into the by-product they produce (note that the equivalent conversion efficiency for assimilation ({f}_{{conv}}) already appears in the growth functions of each genotype, therefore does not appear in Eq. (17)).The residence time ({T}_{{R}_{q}})of each environmental substance is given by the net quantity of this substance divided by the influx, to give units of the average number of biological generations a unit of environmental substance spends in the relevant pool before being removed. In those simulations in which the abiotic influx ({varphi }_{{R}_{q}}(t)) was set to zero (i.e. during the shut-off intervals) production ({P}_{{R}_{q}{NE}T}left(tright)) was used as an alternative denominator:$${IF}left[{varphi }_{{R}_{q}}left(tright) , > , 0,{T}_{{R}_{q}}left(tright)=frac{{R}_{q,{NET}}left(tright)}{{varphi }_{{R}_{q}}left(tright)}right]!,{IF}left[left(left({varphi }_{{R}_{q}}left(tright)=0right){& }left({P}_{{R}_{q}{NET}}left(tright) , > ,0right)right)!,{T}_{{R}_{q}}left(tright)=frac{{R}_{q,{NET}}left(tright)}{{P}_{{R}_{q}{NET}}left(tright)}right]{ELSE}[{T}_{{R}_{q}}left(tright)=0]$$
(19)
The cycling ratio ({{CR}}_{{R}_{q}}) of each substance is given by the ratio between net biotic assimilation of that substance ({A}_{{R}_{q}}left(tright)) and the abiotic influx of that substance ({varphi }_{{R}_{q}}left(tright)). As with the residence time, when the abiotic influx was zero, the input from biological production was used as an alternative denominator:$${IF}left[{varphi }_{{R}_{q}}left(tright) , > , 0,{{CR}}_{{R}_{q}}left(tright)=frac{{A}_{{R}_{q}}left(tright)}{{varphi }_{{R}_{q}}left(tright)}right]{IF}left[left(left({varphi }_{{R}_{q}}left(tright)=0right){{& }}left({P}_{{R}_{q}{NET}}left(tright) , > ,0right)right),{{CR}}_{{R}_{q}}left(tright)=frac{{A}_{{R}_{q}}left(tright)}{{P}_{{R}_{q}{NET}}left(tright)}right]{ELSE}[{{CR}}_{{R}_{q}}left(tright)=0]$$
(20)
Deterministic approximation to steady stateAssume that at steady state substance assimilation will reach a maximal state such that the level of environmental substance is limiting to population size. Assume that such a state is below the level ({{R}_{q,{BIOAVAILABLE}}}_{{crit}}) at which the plastic genotype effectively becomes a second non-producer genotype and can thus be subsumed into non-producer frequency, such that (2) becomes ({S}_{q}left(tright)=mathop{sum }nolimits_{j=1}^{{j}_{{total}}}{N}_{q,j}left(tright)={N}_{q,{prod}}left(tright)+{N}_{q,{non}-{prod}}left(tright)). Assume that there are non-zero starvations at each time-step for all genotypes, which implies growth rate ({G}_{q,j}left(tright) , < , 1,ll {G}_{q,{jMAX}},forall j,q), which gives by (3)({G}_{q,j}left(tright)={G}_{q,j,{PR}}cdot {R}_{q,{BIOAVAILABLE}}left(tright)cdot {f}_{{conv}}). Substituting this into (10), then the first bracketed term in (1), then labeling the post-starvation number of individuals as ({({N}_{q,j}left(tright))}_{{NET}}):$${({N}_{q,j}left(tright))}_{{NET}}=left({N}_{q,j}left(tright)-{rho }_{q,j}left(tright)right)={N}_{q,j}left(tright)cdot left(1-left(1-{G}_{q,j}left(tright)right)right)={N}_{q,j}left(tright)cdot {G}_{q,j}left(tright)$$
(21)
Approximate (14) deterministically by a fixed fractional parameter corresponding to the baseline random death rate:$${delta }_{q,j}left(tright)approx {N}_{q,j}left(tright)cdot {delta }_{0}$$
(22)
Doing the same for mutation:$${{{{{{rm{{Upsilon }}}}}}}}_{q,j}left(tright)={G}_{q,j}left(tright)cdot {left({N}_{q,j}left(tright)right)}_{{NET}}cdot {mu }_{0}={N}_{q,j}left(tright)cdot {{G}_{q,j}left(tright)}^{2}cdot {mu }_{0}$$
(23)
The term in (12) for mutation to (j) from other genotypes simplifies to$${{{{{{rm{{Upsilon }}}}}}}}_{q,x ,ne , j}left(tright)={sum }_{k,=,1}^{{j}_{{total}}}frac{{G}_{q,k}left(tright)cdot {left({N}_{q,k}left(tright)right)}_{{NET}}cdot {mu }_{0}}{{j}_{{total}}-1}={N}_{q,k}left(tright)cdot {{G}_{q,k}left(tright)}^{2}cdot {mu }_{0}$$Substituting Eqs. (21–23) into (1):$${N}_{q,j}left(tright)cdot {{G}_{q,j}left(tright)}^{2}cdot (1-{mu }_{0})-{N}_{q,j}left(tright)cdot {delta }_{0}+{N}_{q,k}left(tright)cdot {{G}_{q,k}left(tright)}^{2}cdot {mu }_{0}=0$$
(24)
Noting that by Eqs. (3)–(5) combined with the above assumptions, the growth rate of the producer can be written as:$${G}_{q,{prod}}left(tright)={G}_{q,{cheat}}left(tright)cdot (1-{kappa }_{{prod},q})$$
(25)
Writing (24) explicitly for each genotype:$${N}_{q,{non}-{prod}}left(tright)cdot {{G}_{q,{non}-{prod}}left(tright)}^{2}cdot (1-{mu }_{0})-{N}_{q,{non}-{prod}}left(tright)cdot {delta }_{0}+{N}_{q,{prod}}left(tright)cdot {{G}_{q,{prod}}left(tright)}^{2}cdot {mu }_{0}=0$$$${N}_{q,{prod}}left(tright)cdot {{G}_{q,{prod}}left(tright)}^{2}cdot (1-{mu }_{0})-{N}_{q,{prod}}left(tright)cdot {delta }_{0}+{N}_{q,{non}-{prod}}left(tright)cdot {{G}_{q,{non}-{prod}}left(tright)}^{2}cdot {mu }_{0}=0$$Adding:$${N}_{q,{non}-{prod}}left(tright)cdot {{G}_{q,{non}-{prod}}left(tright)}^{2}cdot left(1-{mu }_{0}right)-{N}_{q,{non}-{prod}}left(tright)cdot {delta }_{0}$$$$+{N}_{q,{prod}}left(tright)cdot {left({G}_{q,{non}-{prod}}left(tright)cdot left(1-{kappa }_{{prod},q}right)right)}^{2}cdot {mu }_{0}+{N}_{q,{prod}}left(tright)cdot {left({G}_{q,{non}-{prod}}left(tright)cdot left(1-{kappa }_{{prod},q}right)right)}^{2}cdot left(1-{mu }_{0}right)$$$$-{N}_{q,{prod}}left(tright)cdot {delta }_{0}+{N}_{q,{non}-{prod}}left(tright)cdot {{G}_{q,{non}-{prod}}left(tright)}^{2}cdot {mu }_{0}=0$$Because the mutation terms cancel:$$left({N}_{q,{non}-{prod}}left(tright)+{N}_{q,{prod}}left(tright)cdot {left(1-{kappa }_{{prod},q}right)}^{2}right)cdot left({{G}_{q,{non}-{prod}}left(tright)}^{2}-{delta }_{0}right)=0$$
(26)
Substituting in for the growth rate terms (3–7) gives:$$left({N}_{q,{non}-{prod}}(t)+{N}_{q,{prod}}(t)cdot {left(1-{kappa }_{{prod},q}right)}^{2}right)cdot left({left({G}_{0}cdot {R}_{q,{BIOAVAILABLE}}(t)cdot {f}_{{conv}}right)}^{2}-{delta }_{0}right)=0$$
(27)
By (16), (4), (6) and the above, total steady state assimilation of the growth limiting environmental substance by species (q) is:$${A}_{{R}_{q}}left(tright)=mathop{sum }limits_{k=1}^{{j}_{{total}}}frac{{G}_{q,k}left(tright)cdot left({N}_{q,k}left(tright)-{rho }_{q,k}left(tright)right)}{{G}_{q,{kPR}}}$$$$kern2.4pc=frac{{N}_{q,{non}-{prod}}left(tright)cdot {left({G}_{q,{non}-{prod}}left(tright)right)}^{2}}{{G}_{0}}+frac{left({N}_{q,{prod}}left(tright)cdot {left(1-{kappa }_{{prod},q}right)}^{2}right)cdot {left({G}_{q,{non}-{prod}}left(tright)right)}^{2}}{{G}_{0}cdot left(1-{kappa }_{{prod}}right)}$$$$=left({N}_{q,{non}-{prod}}left(tright)+{N}_{q,{prod}}left(tright)cdot left(1-{kappa }_{{prod},q}right)right)cdot {{G}_{0}cdot ({R}_{q,{BIOAVAILABLE}}left(tright)cdot {f}_{{conv}})}^{2}$$
(28)
By (16–18), production of this substance by the producer allele in the other species (p , ne , q), assuming the various arguments above simultaneously apply to this species, is:$${P}_{{R}_{q}}left(tright)= frac{{G}_{p,{prod}}left(tright)cdot left({N}_{p,{prod}}left(tright)-{rho }_{p,{PRODUCER}}left(tright)right)}{{G}_{p,{prodPR}}}cdot {f}_{{conv},{PROD}}\ = {N}_{p,{prod}}left(tright)cdot left(1-{kappa }_{{prod},p}right) cdot {G}_{0}cdot {({R}_{p,{BIOAVAILABLE}}left(tright)cdot {f}_{{conv}})}^{2}cdot {f}_{{conv},{PROD}}$$
(29)
Balance between input and output fluxes of each environmental substance requires ({varphi }_{{R}_{q}}left(tright)+{P}_{{R}_{q}}left(tright)={A}_{{R}_{q}}left(tright)), meaning that by substituting in ({A}_{{R}_{q}}left(tright)) from (28) it is possible to solve for bioavailable substance level, then substitute in the production flux of this substance derived from the producer allele in the other species (p,ne, q):$${R}_{q,{BIO}{AVAILABLE}}left(tright) =frac{1}{{f}!_{{conv}}}sqrt{frac{{varphi }_{{R}_{q}}left(tright)+{P}_{{R}_{q}}left(tright)}{left({N}_{q,{non}-{prod}}left(tright)+{N}_{q,{prod}}left(tright)cdot left(1-{kappa }_{{prod},q}right)right)cdot {G}_{0}}}\ =frac{1}{{f}!_{{conv}}}sqrt{frac{{varphi }_{{R}_{q}}left(tright)+{N}_{p,{prod}}left(tright)cdot left(1-{kappa }_{{prod},p}right)cdot {G}_{0}cdot {({R}_{p,{BIOAVAILABLE}}left(tright)cdot {f}_{{conv}})}^{2}cdot {f}_{{conv},{PROD}}}{left({N}_{q,{non}-{prod}}left(tright)+{N}_{q,{prod}}left(tright)cdot left(1-{kappa }_{{prod},q}right)right)cdot {G}_{0}}}$$
(30)
Substituting this into (27) gives a symmetrical condition for steady state genotype frequencies and substance levels across the system:$$ left({N}_{q,{non}-{prod}}left(tright)+{N}_{q,{prod}}left(tright)cdot {left(1-{kappa }_{{prod},q}right)}^{2}right)cdot\ left(frac{{varphi }_{{R}_{q}}left(tright)+{N}_{p,{prod}}left(tright)cdot left(1-{kappa }_{{prod},p}right)cdot {G}_{0}cdot {({R}_{p,{BIOAVAILABLE}}left(tright)cdot {f}_{{conv}})}^{2}cdot {f}_{{conv},{PROD}}}{left({N}_{q,{non}-{prod}}left(tright)+{N}_{q,{prod}}left(tright)cdot left(1-{kappa }_{{prod},q}right)right)cdot {G}_{0}}-{delta }_{0}right)=0$$
(31)
This solution illustrates the intuitive ideas that growth and reproduction balance random death at steady state and that the associated producer frequency is lower than that of the non-producer by a factor of the cost. (This factor is of second order because the growth rate is used both directly and (by (10)) in the calculation of starvations). Because our model is a discrete stochastic process, (31) can be viewed as an approximation to a steady state condition, subject to the above assumptions combined with the continuous generation of producers by mutation at a sufficient rate to preclude their extinction. The key point is that over long timescales in the finite populations with which we deal, organism-level selection unavoidably favors the non-producer, with no possibility for multi-level fecundity selection. The producer’s stable presence is thus attributable to the combination of mutation and cycle-level selection. More
