Joining

Joint runs

Levels can also be simulated jointly, as the figure shows. Two pairs reach the species tree: a trait or the gene content drives speciation, and the tree is then a result of the run. Three pairs join two different levels along a tree used as input: a trait and a genome, a trait and a sequence, a genome and a sequence. In the remaining three cases the two joined parts belong to the same level: a gene family changes the rates of the genome that carries it, one trait changes how another trait evolves, and each of two genes changes how fast the other one evolves. A joint run holds exactly two parts, and asking for a third is refused. The species tree and a sequence cannot be joined: a sequence evolves along a gene tree, which evolves along the species tree, so joining the two would require simulating the genome as well.

Naming the driver

The connections of a joint run are written with the same verbs and mappings as everywhere else. The difference is the driver. It cannot be passed as a finished result, because it does not exist yet, so it is named instead: "trait", "genomes:toxin", "sequences:rpoB". The name refers to another level of the same run, and the run grows both together. A run with a single trait can call it plainly "trait"; names exist so that a run can hold two. Appendix C lists the named drivers and what each one offers.

A family whose presence drives must be declared in the genome spec, with families=[family("toxin")]. On a tree passed with tree=, a declared family can also have its own duplication, transfer and loss, and these rates can read the other level of the run: with family("A", loss=PerCopy(0.3).scaled_by("trait", {"on": 6.0, "off": 1.0})), A is lost six times faster where the trait is on. A run that simulates the tree reads only the family's name. A gene whose composition drives a rate must declare it with offers=composition("KR", absent=0.02): which letters are counted, and the value used on lineages that carry no copy. A gene that drives a rate usually also declares start=, a second substitution model whose stationary frequencies are used only to found the gene. The gene starts with that composition and then evolves under its own model, so its composition moves toward its own equilibrium. Without start= a gene is founded at its equilibrium, its composition barely changes, and the rate that depends on it stays constant (Jo2).

A joint run is written from the specs of the levels it simulates, for example joint.simulate(species.birth_death(...), traits.discrete(...)). If the species spec is present, the tree is simulated, and the stop condition (n_extant=, total_time=) is written on that spec; if it is not, the tree is passed with tree= and the run ends where the tree ends. The result is a JointResult that carries both simulated levels, the same result objects the ordinary functions return, and both levels write their usual files. When the two joined parts belong to the same level, the run stays on that level's own function, with joint=True: simulate_traits, simulate_genomes_family or simulate_sequences. The flag is checked both ways: asking for it without a live driver is an error, and depending on a live driver without it is an error. One more restriction: a transfer needs the set of lineages alive at an instant, which is unknown while the tree is still being simulated, so a genome run that grows the tree refuses transfer. The gallery presents every joint model in full, with its code (Jo1 to Jo14).

Exact or sliced

A joint run advances by racing every possible event against every other, with the same Gillespie algorithm as everywhere in ZOMBI2. A discrete driver changes only at events, so between two events every rate is constant and the race is exact. When a speciation fires, it uses the trait state of the lineage it happens on; when a trait switches, the tree is untouched; both daughters of a split start from the parent's state.

Two kinds of driver change between events: a continuous trait diffuses at every instant, and a composition moves with every substitution. There is no interval over which the driven rate holds still, so there is nothing exact to race. These models slice time instead. Time is cut into steps of step, and the driver is held fixed within each slice, so the rates are constant inside it. At the slice boundary the driver moves: a diffusing trait by its exact transition law, a composition by the substitutions that accumulated. The approximation is only in the timing: the target uses the value the driver had at the start of the slice.

Model	Tree	Advances by
a trait drives speciation	simulated	events, exact
gene content drives speciation	simulated	events, exact
a diffusing trait drives speciation	simulated	slices of step
a trait and a genome	input	events, exact
two discrete traits, each driving the other	input	events, exact
a continuous trait and another trait, each driving the other	input	slices; the driver is taken at the start of the slice
a trait and a gene's sequence	input	slices; the composition is taken at the start of the slice
a genome and a gene's sequence	input	slices; the composition is held within a slice
a gene family drives its own genome	input	events, exact
two genes, each depending on the other	input	slices; both compositions held within a slice

step is written on the link, because the right size depends on how the value is taken: a steep curve needs a finer step than a flat one. In a joint run it has no default, and leaving it out is an error, because the timescale belongs to a model that does not exist yet.

A sliced run is first order in step: halve step and the error halves. Halving step also changes the numbers for a given seed, because it changes how many draws the run takes. The two runs are samples from two nearby laws, not one answer computed twice, so comparing them tells you nothing. Check convergence across seeds instead (Jo14):

1. pick a summary the model is about — the mean value at the tips, the variance at the tips, the correlation between two traits at the tips, the fraction of tips in a given state;
2. run a set of seeds, say fifty, at step, then the same seeds at step / 2, then at step / 4;
3. for each step, take the mean of that summary over the seeds, and its standard error over the seeds;
4. compare the shift between consecutive steps against that standard error. When the shift is smaller, the step is fine enough to report; when it is larger and still falling by about half each time, halve again.

The third step is what makes this a test rather than an impression. A single seed at two steps differs by the Monte Carlo spread whatever the step is, and that spread is usually far larger than the error you are trying to see.

On the command line

zombi2 joint covers the two models that drive speciation from a discrete driver. The rate is written exactly as in Python, and the remaining flags are the ones zombi2 traits or zombi2 genomes would take. Giving flags from both levels at once is an error.

zombi2 joint out/ --death 0.2 --states small,large --switch 0.3 \
    --n-extant 100 --seed 1 \
    --birth "PerLineage(1.0).scaled_by('trait', {'small': 1.0, 'large': 3.0})"

zombi2 joint out/ --origination 0.2 --loss 0.1 --family-names toxin \
    --n-extant 60 --seed 1 \
    --birth "PerLineage(1.0).scaled_by('genomes:toxin', {'present': 3.0, 'absent': 1.0})"

The other models are Python only, deliberately: each needs something no flag can carry, such as a curve, a rate for a single named family, a list of specs, or a start model built as a matrix.

Literature

The state-dependent diversification models are usually known by their acronyms. This table gives the correspondence.

What it does	From the literature	Gallery
a binary trait drives speciation (and extinction)	BiSSE (Maddison et al. 2007)	Jo10
a multi-state trait drives speciation	MuSSE (FitzJohn 2012)	Jo12
a trait drives speciation and jumps at the split	ClaSSE (Goldberg & Igić 2012)	Jo13
a continuous trait drives speciation	QuaSSE (FitzJohn 2010)	Jo14