(A) Schematic of the five different model features (F1-F5) used in the four models considered (see Table 1). F1: assumes the burst frequency oscillates over 24 hours; F2: assumes a fixed burst size for each gene. The burst frequency (f) and burst size (b) determine the shape of the NB distribution used to model the mRNA counts (see Methods), where the mean μ = bf and the variance σ2= μ + μ2/f. F3: the burst size in each cell scales with measured cell area using a gene-specific exponent βg. F4: to model incomplete synchronisation, each cell has a cell-specific phase offset compared to the average phase. F5: the burst size in each cell is proportional to the cell area multiplied by a cell-specific burst parameter εi, g. These cell-specific εi, g are further modelled using a bivariate log-normal distribution. (B) PSIS-LOO calculated for each of the four models using the whole dataset of both gene-pairs across all time points. 