READ ME This file illustrates the R code for numerical simulations in ‘Evolution of fixed demographic heterogeneity from a game of stable coexistence’ by Giaimo, Li, Traulsen and Baudisch in Demographic Research. In what follows, ppm stands for population projection matrix. payoff.m: stores as a 2 × 2 matrix the values of the payoff matrix, Eq. (1) in the manuscript, assumed for the game. equilibrium.m: stores as a scalar the value of equilibrium frequency of A in the coexistence game in absence of demographic structure for the given payoff matrix det.payoffmatrix.m: stores as a scalar the value of the determinant of the given payoff matrix expected.payoffs: function that takes as argument a scalar x indicating the frequency of A in the population and returns a two components column vector with first component and the second component are the expected payoffs to A and B, respectively, when the frequency of A is x. make.Leslie: function that returns a N × N Leslie matrix from two arguments: survival, a (N − 1)-components vector of age-specific survival probabilities, i.e. subdiagonal of the Leslie matrix; fertility, a N-components vector of age-specific fertilities, i.e. first row of the Leslie matrix. P.x: stores the time-independent survival probabilities of A for the simulation in Fig. 2 of the manuscript. M.x: stores the time-independent fertilities of A for the simulation in Fig. 2 of the manuscript. L.base: stores the values of the (Leslie) ppm shared by A and B in Fig. 2 of the manuscript as a 6 × 6 matrix . lambda: function that takes a ppm as argument and returns a scalar corresponding to the leading eigenvalue of the matrix. u.x: function that takes a ppm as argument and returns the right leading eigenvector scaled so that the sum of the components equals 1, i.e. stable stage distribution. v.x: function that takes a ppm as argument and returns the left leading eigenvector scaled so that its inner product with the appropriately scaled right leading eigenvector equals 1. birth.rate: function that takes a ppm as argument and returns the birth rate of the population at the demographic equilibrium implied by the ppm. The computation is based on Eq. (37) in the manuscript. matw.payoff: function that takes a ppm and a scalar (i.e. payoff) as arguments and returns the same ppm in which payoff has been added to each first row entry. transfer: function that takes as arguments a ppm and a scalar (i.e. x). The scalar x corresponds to the leading eigenvalue one wants to the ppm to have by adding entrywise to the first row a scalar δ of unknown magnitude. The function returns the value of the transfer function at x for the given ppm, i.e. G(x), see Section 4.1 of the manuscript. As the transfer function satisfies the relation δG(x)=1 for the required x, the output of the function is 1/δ. inverse.hat.x.A: function that takes as arguments a scalar hat.z and the time-independent ppm of type A. It returns a scalar corresponding to the relative frequency x of A required so that the (time-dependent) ppm of type A as a function of x (i.e. with fertilities including the expected payoff from the game) has leading eigenvalue equal to hat.z. equilibrium.equation: function that takes as arguments a scalar z, the time-independent ppm of type A and the time independent ppm of type B. It returns 1 minus the value of the right hand side of Eq. (23) in the manuscript. Roots of this function correspond to the growth rate shared by the A and B population at the interior equilibrium, i.e. coexistence, for the assumed payoff matrix. outofboundpayoff: stores the expected payoffs to both types for values of the relative frequency of A slightly above 1 and slightly below 0. fixed.point.age.str: function that takes as arguments the time-independent ppm of type A and the time independent ppm of type B. It returns, if it exists, the (assumed unique) root of the function equilibrium.equation. singlevtperturbed.eq: function that takes as arguments a ppm and a vector perturbation.range. The argument ppm is assumed to be the shared time-independent ppm of A and B and having the structure of a Leslie matrix. The vector contains the values of the perturbations to apply to each nonzero entry of ppm separately. The function returns a list with a number of elements equal to the number of nonzero entries of ppm, i.e.2N − 1 for a N × N Leslie matrix. Each element of the list is a vector with the same number of components as the vector perturbation.range. Element i of the list refers to the i nonzero entry of ppm, where such entries are numbered by reading the matrix row-wise starting from the first row. The component j of the vector in element i of the list corresponds to the equilibrium frequency of A at coexistence with B under the assumed payoff matrix when B has time-independent ppm equal to ppm and A has time-independent ppm equal to ppm with the i nonzero element being perturbed by an amount equal to component j of perturbation.range. If the i element of the list refers to a survival probability (i.e. i > N), the perturbation is multiplicative. If the i element of the list refers to a fertility (i.e. i ≤ N), the perturbation is additive. tangent.lines: function that takes as argument a N × N ppm with the structure of a Leslie matrix and returns a vector where component i is the leading eigenvalue sensitivity to an additive perturbation in entry (1, i) if N ≥ i or to a multiplicative perturbation in entry (i + 1, i) if i > N. par.ppm: function that takes as arguments 4 separate scalars: max.age, a, b and h. It returns a N × N Leslie matrix where N =max.age, age-specific survival is computed from the a and b Gompertz parameters, see Eq. (35) in the manuscript, and age-specific fertility is computed from max.age and h (= c in the manuscript), see Eq. (36). eq.parmat: function that takes as arguments 8 separate scalars: max.age.A, a.A, b.A, h.A, max.age.B, a.B, b.B and h.B. Calling par.ppm it constructs the ppms of A and B from the given parameters. Then it calls fixed.point.age.str and returns a scalar corresponding to the equilibrium frequency of A at coexistence. gompertz.mortality: function that takes as arguments a, b and x. From the Gompertz parameters a and b, it returns a scalar corresponding to mortality at age x. gompertz.survivorship: function that takes as arguments a, b and x. After calling gompertz.mortality, it returns a scalar corresponding to survivorship at age x for the Gompertz parameters a and b. average.cohort.mort: function that takes as arguments 8 separate scalars:max.age.A, a.A, b.A, h.A, max.age.B, a.B, b.B, h.B and y. Calling eq.parmat it constructs the ppms of A and B from the given parameters and retrieves the equilibrium frequency xˆ of A at coexistence. The function returns a 2 elements list with xˆ as first element and as a second element a 3-components vector: mortality at age y in a cohort of 1-aged A individuals, mortality at age y in a cohort of 1-aged B individuals and mortality at age y in a mixed cohort of 1-aged A and B individuals where the initial proportions of types reflect the equilibrium birth rates of A and B at coexistence, see Eq. (38). total.pop: stores the initial total population size for the simulation in Fig. 2 of the manuscript. max.time.steps: stores the maximum number of iteration of the simulation in Fig. 2 of the manuscript. mat.A: stores the initial time-independent ppm of A, set equal to L.base for the simulation in Fig. 2 of the manuscript. mat.B: stores the initial time-independent ppm of B, set equal to L.base for the simulation in Fig. 2 of the manuscript. S.x.Aprime: stores the mutated time-independent survival probabilities of A for the simulation in Fig. 2 of the manuscript. M.x.Aprime: stores the mutated time-independent fertilities of A for the simulation in Fig. 2 of the manuscript. mat.Aprime: stores the mutated time-independent ppm of A for the simulation in Fig. 2 of the manuscript. A.x: function that takes as arguments the scalars num.A and num.B which are the abundances of A and B, respectively, and returns a scalar representing the relative frequency of A. A.x.current: stores a vector with last element equal to the current frequency of A for the simulation in Fig. 2 of the manuscript. n.A: stores the initial size of A population for the simulation in Fig. 2 of the manuscript. n.B: stores initial size of B population for the simulation in Fig. 2 of the manuscript. # loop: loop that iterates for max.time.steps times the dynamics in Eq. (12) of the manuscript up to t = 500, then A mutates and acquires ppm equal to mat.Aprime. The dynamics are then those in Eq. (16). At every iteration, the loop updates the vector A.x.current by appending at the end of it the current frequency of A.