Modeling gene content across a phylogeny to determine when genes become associated

Abstract

We consider a model for inferring functional links between genes. We begin with the simple case of two genes whose presence or absence evolves stochastically along a phylogenetic tree. We develop a hidden Markov model where the hidden states of the model correspond to whether or not the genes perform a joint function. In the case that two genes do perform a joint function, the rates of gain or loss of each gene depend on the presence or absence of the other gene. Otherwise, those two genes are assumed to be gained and lost independently. Using simulation, we investigate the conditions under which the package corHMM can infer the hidden state correctly, and we also investigate when the Akaike information criterion (AIC) has the power to reject the simpler model when it is incorrect. We find that we can more accurately determine the dependent and independent rate class regimes when the trees have more tips and when the differences in the two rate classes are larger. We find the accuracy of corHMM is not overly affected by whether or not there are multiple transitions between the rate classes or just a single transition. We show how the two-gene case can be extended to a more general n-gene model with a level-dependent quasi-birth-and-death (LD-QBD) framework. We assume that the level n of the QBD corresponds to the number of genes that are required for some beneficial function, and the phases within each level record the presence/absence of particular genes.

Keywords:

• Hidden Markov models
• level-dependent quasi-birth-and-death model
• neofunctionalization
• stochastic mapping

1. Introduction

In this work, we develop a stochastic model of gene gain and loss with the aim of inferring when (if at all), in evolutionary history, an association between two (or more) genes arises. The data we consider is a species tree along with information on the presence or absence of genes in each species. One biological motivation for our model is that if genes are involved in the same biochemical pathway and therefore they are both required for some function, then the rate of gain or loss of one gene in the pathway should depend upon the presence or absence of the other gene in the pathway[12]. However, if genes are not functionally linked, then the rate of gain or loss of one gene should be independent of the state of the other genes[24, 25]. Functional linkage could also occur if multiple genes code for subunits that form oligomeric proteins[17].

We do not explicitly model the mechanism by which genes are gained, but the assumptions of our model are consistent with some of the ways new genes are thought to arise in nature. Genes can be gained via de novo gene birth, or by gene duplication followed by subsequent subfunctionalization or neofunctionalization (sometimes known as the patchwork model[30]), or by horizontal gene transfer (HGT). Of these mechanisms, de novo gene gain seems like it would be highly unlikely to occur independently in different evolutionary lineages, but there is evidence for convergent gene gain via both HGT[2, 15] and gene duplication followed by subsequent subfunctionalization or neofunctionalization[11].

We are not the first to consider mathematical models of how genes gain and maintain interactions. Francis and Tanaka[12] developed a model for the evolution of bacterial genetic pathways (where genes were gained via HGT) under the assumption that the complete pathway was required to confer any selective benefit. They considered the dynamics of such a system in a single species rather than looking at patterns of gene presence and absence over a phylogeny. Lynch proposed models[17, 18] for the evolution of multimeric protein assemblages - his models focused on population genetic aspects such as mutation and drift rather than looking at a phylogenetic scale.

The model proposed in this article takes a more phylogenetic approach compared to the models mentioned above but could be applied in similar biological settings. It builds upon some existing models for the evolution of discrete traits. Barker et al.[3] and Barker and Pagel[4] were the first to note that the presence and absence of two genes can be modeled as a four-state Markov process. We assume that states $\{00,01,10,11\}$ record each gene’s presence (1) or absence (0). Each transition can only have one gene loss or gain (e.g. 00→01 or 00→10 are permitted but not 00→11). To test if two genes were associated, Barker et al.[3] and Barker and Pagel[4] used likelihood ratio tests to compare a model where substitution rates were consistent with independent evolution (e.g., $q_{00\to 01}=q_{10\to 11}$) versus a model where rates were dependent. Marazzi et al.[19] proposed a precursor model for the evolution of a binary trait. Instead of transiting directly from 0→1 or 1→0, they developed a model with an intermediate unobservable state. Boyko and Beaulieu[7] further generalized this idea to models with hidden rate classes. For example, two different rate classes may have the same observable states but different transition rates among the states.

The following are our contributions in this article. First, in Section 2, we construct a two-gene model based on the dependent and independent models considered by Barker et al.[3] as two hidden rate classes in the framework proposed by Boyko and Beaulieu[7]. We simulate data under this two-gene model and determine under what conditions a shift from the independent rates class to the dependent rates class can be detected using a hidden Markov model as implemented in the R package corHMM from Boyko and Beaulieu[7]. For example, how large a tree is required and how large a shift in the rates is needed before the Akaike information criterion (AIC) supports a model with two rate classes over a simpler model with just one rate class? If a model with two rate classes is preferred, can it correctly detect where on the evolutionary tree the shift (neofunctionaliszation) occurred?

Further, in Section 3, building on the theory of the quasi-birth-and-death models (QBDs) in matrix-analytic methods[16, 21], we present an extension of our two-gene model to the case where an association arises between more than two genes. We construct a level-dependent QBD (LD-QBD) to model the processes by which the association between N−1 genes can be extended to N genes for some N > 2. In our model, we assume that the rates of gene gain or loss may depend on the number of genes that have become associated. We then present numerical examples to illustrate biological insights that may be gained from the analysis of the stationary distribution.

Finally, in Section 4, we discuss a wide range of possible applications and possible extensions of our model to the analysis of biological gene presence-absence data.

2. Neofunctionalization two-gene model I

First, we propose the following continuous-time Markov chain $\{X(t):t\ge 0\}$ with state space (1) $$\begin{array}{ll}\mathcal{S}& =\{00,01,10,11,\widetilde{00},\widetilde{01},\widetilde{10},\widetilde{11}\},\end{array}$$(1) where 1 and 0 correspond to the presence or absence of a gene, states in the subset ${\mathcal{S}}_{ind}=\{00,01,10,11\}\subset \mathcal{S}$ correspond to the case where the gain/loss of each gene does not depend on the presence/absence of the other gene, and states in the subset ${\mathcal{S}}_d=\{\widetilde{00},\widetilde{01},\widetilde{10},\widetilde{11}\}\subset \mathcal{S}$ correspond to genes that do depend on one another. For example, state $X(t)=\widetilde{11}$ means that both genes are present in the species, and are also dependent on one another, at time t.

The transition rates of this process are collected in the generator matrix $\mathbf{A}=[A_{i,j}{]}_{i,j\in \mathcal{S}}$ given by, (2) $$\begin{array}{ll}\mathbf{A}& =\begin{array}{cc}& \begin{array}{cccccccc}00& 01& 10& 11& \widetilde{00}& \widetilde{01}& \widetilde{10}& \widetilde{11}\end{array}\\ \begin{array}{c}00\\ 01\\ 10\\ 11\\ \widetilde{00}\\ \widetilde{01}\\ \widetilde{10}\\ \widetilde{11}\end{array}& \left[\begin{array}{cccccccc}\ast & q_{12}& q_{13}& 0& 0& 0& 0& 0\\ q_{21}& \ast & 0& q_{24}& 0& 0& 0& 0\\ q_{31}& 0& \ast & q_{34}& 0& 0& 0& 0\\ 0& q_{42}& q_{43}& \ast & 0& 0& 0& \lambda \\ {\text{hlineq}}_1& 0& 0& 0& \ast & \widetilde{q_{12}}& \widetilde{q_{13}}& 0\\ 0& q_2& 0& 0& \widetilde{q_{21}}& \ast & 0& \widetilde{q_{24}}\\ 0& 0& q_3& 0& \widetilde{q_{31}}& 0& \ast & \widetilde{q_{12}}\\ 0& 0& 0& q_4& 0& \widetilde{q_{42}}& \widetilde{q_{43}}& \ast \end{array}\right]\end{array},\end{array}$$(2) where 0 entries correspond to events that may not occur; the rates of transitions from being dependent to being independent are set to $q_1,q_2,q_3,q_4$; the transition from independence to dependence may only occur from state 11 to state $\widetilde{11}$, at some rate denoted λ; and the ondiagonals, marked with *, are given by $A_{i,i}=-\sum _{j\ne i}{A}_{i,j}$, for all $i,j\in \mathcal{S}$. An alternative assumption that could be made about the transition from an independent to a dependent state is that a change in selective pressure could occur regardless of whether both genes are present, in this case, the diagonal entries in the upper right block of matrix A would all be set to λ. This latter assumption might be more appropriate in the case where genes are gained via horizontal transfer, whereas the assumption we make here perhaps fits better with gene gain being the result of duplication followed by subsequent subfunctionalization or neofunctionalization[30].

We assume that genes are gained or lost at some Poisson rates α_k and β_k, per gene k. Therefore, when in one of the ‘independent-genes’ states in the set 𝒮_ind, we have ${\alpha }_1=q_{13}=q_{24}$ and ${\alpha }_2=q_{12}=q_{34}$, which are some Poisson rates at which species may gain gene 1 and gene 2, respectively; while ${\beta }_1=q_{31}=q_{42}$ and ${\beta }_2=q_{21}=q_{43}$, which are some Poisson rates at which species may lose gene 1 and gene 2, respectively.

For our simulation study, we make some further simplifying assumptions about the model in 2(2) $$\begin{array}{ll}\mathbf{A}& =\begin{array}{cc}& \begin{array}{cccccccc}00& 01& 10& 11& \widetilde{00}& \widetilde{01}& \widetilde{10}& \widetilde{11}\end{array}\\ \begin{array}{c}00\\ 01\\ 10\\ 11\\ \widetilde{00}\\ \widetilde{01}\\ \widetilde{10}\\ \widetilde{11}\end{array}& \left[\begin{array}{cccccccc}\ast & q_{12}& q_{13}& 0& 0& 0& 0& 0\\ q_{21}& \ast & 0& q_{24}& 0& 0& 0& 0\\ q_{31}& 0& \ast & q_{34}& 0& 0& 0& 0\\ 0& q_{42}& q_{43}& \ast & 0& 0& 0& \lambda \\ {\text{hlineq}}_1& 0& 0& 0& \ast & \widetilde{q_{12}}& \widetilde{q_{13}}& 0\\ 0& q_2& 0& 0& \widetilde{q_{21}}& \ast & 0& \widetilde{q_{24}}\\ 0& 0& q_3& 0& \widetilde{q_{31}}& 0& \ast & \widetilde{q_{12}}\\ 0& 0& 0& q_4& 0& \widetilde{q_{42}}& \widetilde{q_{43}}& \ast \end{array}\right]\end{array},\end{array}$$(2) . We assume that both genes are gained/lost at the same rate. We further assume that the transition from dependence to independence occurs at a constant rate, denoted by $q_1=q_2=q_3=q_4=\gamma$.

Finally, we denote by $r_{00},r_{01},r_{10},r_{11},{\overline{r}}_{01},{\overline{r}}_{10}>0$ the parameters which modify the base rates of transition between states within the set 𝒮_d, so that given some state i∈𝒮_d, gene $k=1,2$ is gained at rate α_kr_i or lost at rate ${\beta }_k{\overline{r}}_i$.

By the above assumptions, the generator matrix 𝐀 of the resulting process, referred to as Model I, is given by (3) $$\begin{array}{l}\mathbf{A}=\begin{array}{cc}& \begin{array}{cccccccc}00& 01& 10& 11& \widetilde{00}& \widetilde{01}& \widetilde{10}& \widetilde{11}\end{array}\\ \begin{array}{c}00\\ 01\\ 10\\ 11\\ \widetilde{00}\\ \widetilde{01}\\ \widetilde{10}\\ \widetilde{11}\end{array}& \left[\begin{array}{cccccccc}\ast & {\alpha }_2& {\alpha }_1& 0& 0& 0& 0& 0\\ {\beta }_2& \ast & 0& {\alpha }_1& 0& 0& 0& 0\\ {\beta }_1& 0& \ast & {\alpha }_2& 0& 0& 0& 0\\ 0& {\beta }_1& {\beta }_2& \ast & 0& 0& 0& \lambda \\ \gamma & 0& 0& 0& \ast & {\alpha }_2{r}_{00}& {\alpha }_1{r}_{00}& 0\\ 0& \gamma & 0& 0& {\beta }_2{r}_{01}& \ast & 0& {\alpha }_1{\overline{r}}_{01}\\ 0& 0& \gamma & 0& {\beta }_1{r}_{10}& 0& \ast & {\alpha }_2{\overline{r}}_{10}\\ 0& 0& 0& \gamma & 0& {\beta }_1{r}_{11}& {\beta }_2{r}_{11}& \ast \end{array}\right].\end{array}\end{array}$$(3)

In the mutation-selection framework, beneficial mutations are more likely to fix in a population, so the rate of substitution of a beneficial mutation should be higher[13, 20, 22]. Similarly, deleterious mutations are less likely to be fixed in a population, so the rate of losing one of those genes should be lower[8, 23, 32]. We consider that when two genes become associated, i.e., both are needed to perform some useful function, then a transition from state 11 to state $\widetilde{11}$ occurs. Since the function is useful, we assume that the rate of losing a gene, when both genes are present ($\widetilde{11}$ moves to $\widetilde{10}$ or $\widetilde{01}$), would be substantially less than the rate of losing a gene when one of the genes is already absent ($\widetilde{10}$ moves to $\widetilde{00}$, or $\widetilde{01}$ moves to $\widetilde{00}$). Therefore, r₁₁≪1 and for similar reasons ${\overline{r}}_{01},{\overline{r}}_{10}>1$. In contrast, we expect r₀₀, r₀₁, and r₁₀ to be ≈1. The reason is that when one of the genes is absent, the useful function that the genes perform together has already been lost, so a further loss will not be protected by selection compared to the independent case. Here we enforce that these rates (r₀₀, r₀₁, and r₁₀) are exactly equal to 1 as we assume that any other shifts in selective pressure between the dependent and independent regimes will be minor, furthermore, it is advantageous to reduce the number of parameters that corHMM needs to estimate.

2.1. Simulation study

For all simulations and example figures in this section, species trees were simulated under a birth-and-death process with constant rates using the sim.bdtree function from the R package geiger[26]. The speciation rate was set to 1, and the extinction rate was 0. 7. Given a species tree, we then randomly generated the presence or absence of genes for each species on the tree, using two different transition rate classes (𝐐_Ind and 𝐐_Dep). In our first two sets of simulations, we picked a single node in the tree where a transition from the independent to the dependent regime occurred. The node was selected randomly, but with the constraint that the descendants of the node were no less than 25% of the total taxa and no more than 75%.

We considered two main measures of performance. The first was whether a Hidden Markov Model (HMM) following the structure given in Equation (3) gave a better fit than a standard Markov model with only one transition rate class. We compared it to a standard Markov model with one transition rate class that either had the same structure as 𝐐_Ind or 𝐐_Dep. We compared the AICc values between a two-rate classes HMM model and a standard Markov model following Anderson and Burnham’s[1] interpretation of the difference in AICc between models, where a difference of 2 or more indicates that the model with the hidden rate classes is considered significantly better compared to a simpler model.

The second measure of performance is related to whether the transition points between the independent and dependent regimes can be identified. To infer the ancestral states, we first applied the HMM model in corHMM to estimate the transition rates in 𝐀. We then applied the function makeSimmap from the corHMM package to simulate the ancestral states along the tree, given the inferred transition rates. For our model, the function makeSimmap assigns one of 8 different states to each segment of the branches. However, we only take note of whether or not the state is in the independent/dependent rate class. Although the results of makeSimmap are stochastic, we found that multiple runs produced similar results, so we report only one “simmap” per simulation. We computed the percentage of matching branch length between the true generating model and the “simmap” based on the inferred model.

To illustrate the structure of the data we generate, in Figure 1, we present an example of a 50-taxon species tree and the presence or absence of two genes. In the blue portion of the tree shown in Figure 1, the two genes are associated and having them both is beneficial. This means that the rate of loss of either gene from a state with both genes is very low, in this example, all tip species in this portion of the tree have both genes. In the rest of the tree, the genes are not associated and we see all possible combinations of presence and absence.

Figure 1. A phylogeny with the presence or absence of two genes shown at the tips. Gene 1 (circles) and gene 2 (squares), are either absent (red) or present (black) in each of the species at the tips of the tree. Black edges correspond to the parts of the tree where the genes are independent of each other, while blue edges indicate the two genes are dependent.

In their analysis of the corHMM package, Boyko and Beaulieu[7] found that datasets simulated with a mean rate of 0. 1 for transitions between hidden rate classes and between observable states gave more accurate results than simulations using higher mean rates. In keeping with this, our first simulation used the generator matrices 𝐐_Ind, 𝐐_Dep given by (4) $$\begin{array}{ll}{\mathbf{Q}}_{Ind}& =\left[\begin{array}{cccc}\ast & 0.04& 0.04& 0\\ 0.03& \ast & 0& 0.04\\ 0.03& 0& \ast & 0.04\\ 0& 0.03& 0.03& \ast \end{array}\right],\\ {\mathbf{Q}}_{Dep}& =\left[\begin{array}{cccc}\ast & 0.04& 0.04& 0\\ 0.03& \ast & 0& 0.4\\ 0.03& 0& \ast & 0.4\\ 0& 0.00375& 0.00375& \ast \end{array}\right].\end{array}$$(4)

These matrices correspond to values of ${\overline{r}}_{01},{\overline{r}}_{10}=10$ and r₀₁ = 0. 125 in Equation (3(3) $$\begin{array}{l}\mathbf{A}=\begin{array}{cc}& \begin{array}{cccccccc}00& 01& 10& 11& \widetilde{00}& \widetilde{01}& \widetilde{10}& \widetilde{11}\end{array}\\ \begin{array}{c}00\\ 01\\ 10\\ 11\\ \widetilde{00}\\ \widetilde{01}\\ \widetilde{10}\\ \widetilde{11}\end{array}& \left[\begin{array}{cccccccc}\ast & {\alpha }_2& {\alpha }_1& 0& 0& 0& 0& 0\\ {\beta }_2& \ast & 0& {\alpha }_1& 0& 0& 0& 0\\ {\beta }_1& 0& \ast & {\alpha }_2& 0& 0& 0& 0\\ 0& {\beta }_1& {\beta }_2& \ast & 0& 0& 0& \lambda \\ \gamma & 0& 0& 0& \ast & {\alpha }_2{r}_{00}& {\alpha }_1{r}_{00}& 0\\ 0& \gamma & 0& 0& {\beta }_2{r}_{01}& \ast & 0& {\alpha }_1{\overline{r}}_{01}\\ 0& 0& \gamma & 0& {\beta }_1{r}_{10}& 0& \ast & {\alpha }_2{\overline{r}}_{10}\\ 0& 0& 0& \gamma & 0& {\beta }_1{r}_{11}& {\beta }_2{r}_{11}& \ast \end{array}\right].\end{array}\end{array}$$(3) ). We generated 100 species trees, each with 150, 200, or 250 tips. The results of the simulation are shown in Figure 2.

Figure 2. Panel (a) shows the AICc difference between the HMM with two hidden rate classes and a standard Markov with the same structure as 𝐐_Ind. Panel (b) shows the AICc difference between the HMM and a standard Markov model with the same structure as 𝐐_Dep. Panel (c) shows the percentage of the edge length where the state of dependence or independence was correctly assigned when comparing between the true generating model and a “simmap” using the inferred transition rates in the HMM. Results are shown for 100 trees each of size of 150, 200 and 250. The true generating model used rates from Equation (4(4) $$\begin{array}{ll}{\mathbf{Q}}_{Ind}& =\left[\begin{array}{cccc}\ast & 0.04& 0.04& 0\\ 0.03& \ast & 0& 0.04\\ 0.03& 0& \ast & 0.04\\ 0& 0.03& 0.03& \ast \end{array}\right],\\ {\mathbf{Q}}_{Dep}& =\left[\begin{array}{cccc}\ast & 0.04& 0.04& 0\\ 0.03& \ast & 0& 0.4\\ 0.03& 0& \ast & 0.4\\ 0& 0.00375& 0.00375& \ast \end{array}\right].\end{array}$$(4) ) and had a single transition from independence to dependence.

We found that, while it was possible to correctly identify that a model with two hidden rate classes was superior to a model with a single rate class (Figure 2, panels a and b), it was not possible to accurately identify where in the tree the transition occurred (Figure 2, panel c).

In their study, Barker et al.[3] investigated rates of gene gain or loss in the range [0. 1, 6] (per unit of branch length). They found that the best values for being able to identify an association between genes were in the range [0. 8, 1. 5]. In our second set of simulations, we used higher rates of gene gain and loss to investigate if it would lead to better performance in detecting where transitions in the hidden rate classes occur. We chose the rates of gene gain or loss as follows: (5) $$\begin{array}{l}{\mathbf{Q}}_{Ind}=\left[\begin{array}{cccc}\ast & 3& 3& 0\\ 2& \ast & 0& 3\\ 2& 0& \ast & 3\\ 0& 2& 2& \ast \end{array}\right],{\mathbf{Q}}_{Dep}=\left[\begin{array}{cccc}\ast & 3& 3& 0\\ 2& \ast & 0& 8h\\ 2& 0& \ast & 8h\\ 0& \frac{1}{2h}& \frac{1}{2h}& \ast \end{array}\right],\end{array}$$(5) where $h=1,2,3$. In terms of Equation (3(3) $$\begin{array}{l}\mathbf{A}=\begin{array}{cc}& \begin{array}{cccccccc}00& 01& 10& 11& \widetilde{00}& \widetilde{01}& \widetilde{10}& \widetilde{11}\end{array}\\ \begin{array}{c}00\\ 01\\ 10\\ 11\\ \widetilde{00}\\ \widetilde{01}\\ \widetilde{10}\\ \widetilde{11}\end{array}& \left[\begin{array}{cccccccc}\ast & {\alpha }_2& {\alpha }_1& 0& 0& 0& 0& 0\\ {\beta }_2& \ast & 0& {\alpha }_1& 0& 0& 0& 0\\ {\beta }_1& 0& \ast & {\alpha }_2& 0& 0& 0& 0\\ 0& {\beta }_1& {\beta }_2& \ast & 0& 0& 0& \lambda \\ \gamma & 0& 0& 0& \ast & {\alpha }_2{r}_{00}& {\alpha }_1{r}_{00}& 0\\ 0& \gamma & 0& 0& {\beta }_2{r}_{01}& \ast & 0& {\alpha }_1{\overline{r}}_{01}\\ 0& 0& \gamma & 0& {\beta }_1{r}_{10}& 0& \ast & {\alpha }_2{\overline{r}}_{10}\\ 0& 0& 0& \gamma & 0& {\beta }_1{r}_{11}& {\beta }_2{r}_{11}& \ast \end{array}\right].\end{array}\end{array}$$(3) ) this corresponds to ${\overline{r}}_{01},{\overline{r}}_{10}=\frac{8}{3}h$ and $r_{11}=\frac{1}{4h}$. Results relating to model selection are shown in Figure 3. Similarly to our first set of simulations, the AICc differences in Figure 3 show that the two rate-class HMMs are usually significantly better than a standard Markov model.

Figure 3. The AICc difference between the HMM and the standard Markov models under different values of h. (a) shows the AICc difference between the HMM and a standard Markov, which has the same structure as 𝐐_Ind, and (b) shows the AICc difference between the HMM and a standard Markov model, which has the same structure as 𝐐_Dep The dashed horizontal line is the value equal to 2.

We were particularly interested to see if larger rates of gene gain and loss would lead to better performance in identifying where on the tree transitions between the independent and dependent rate classes occurred. Figure 4 shows two contrasting results, in (a) and (b), there is a 94. 4% match between the true and inferred states with respect to dependence and independence, however, in (c) and (d) there is only a 42. 3% match. These two examples were both generated with the same choice of 𝐐_Ind and 𝐐_Dep and illustrate that there is a lot of variability in how well we are able to recover the hidden state.

Figure 4. Results of two simulation instances using the values of 𝐐_Ind, 𝐐_Dep in Equation (5(5) $$\begin{array}{l}{\mathbf{Q}}_{Ind}=\left[\begin{array}{cccc}\ast & 3& 3& 0\\ 2& \ast & 0& 3\\ 2& 0& \ast & 3\\ 0& 2& 2& \ast \end{array}\right],{\mathbf{Q}}_{Dep}=\left[\begin{array}{cccc}\ast & 3& 3& 0\\ 2& \ast & 0& 8h\\ 2& 0& \ast & 8h\\ 0& \frac{1}{2h}& \frac{1}{2h}& \ast \end{array}\right],\end{array}$$(5) ) with h = 3. Panels (a) and (c) are examples of the true generating model. Panels (b) and (d) are examples of a character history simulated under the HMM using “simmap”. Black edges correspond to the times when the genes are independent of each other (transition rates in 𝐐_Ind), and blue edges indicate the dependence of the two genes (transition rates in 𝐐_Dep). In panel (b) and (d) black edges indicate the inferred independence of the two genes, and the red lines correspond to the inferred dependence.

Figure 5 indicates that bigger differences between Q_Ind and Q_Dep make it easier to accurately assign ancestral states. Figure 6 shows that when h = 3 it is possible to accurately infer which edges are in the dependent versus independent state, provided the number of taxa is more than 80. We note that more taxa might be required for different values of h.

Figure 5. The percentage of the edge-length where the state of dependence or independence was correctly assigned when comparing between the true generating model and a “simmap” using the inferred transition rates in the HMM. Results are shown for 100 trees of size 50. The true generating model used rates from Equation (5(5) $$\begin{array}{l}{\mathbf{Q}}_{Ind}=\left[\begin{array}{cccc}\ast & 3& 3& 0\\ 2& \ast & 0& 3\\ 2& 0& \ast & 3\\ 0& 2& 2& \ast \end{array}\right],{\mathbf{Q}}_{Dep}=\left[\begin{array}{cccc}\ast & 3& 3& 0\\ 2& \ast & 0& 8h\\ 2& 0& \ast & 8h\\ 0& \frac{1}{2h}& \frac{1}{2h}& \ast \end{array}\right],\end{array}$$(5) ) with $h\in \{1,2,3\}$. The red dashed line at 0. 5 is the expected match if we randomly assign the ancestral state (i.e., guessing).

Figure 6. The percentage of the edge length where the state of dependence or independence was correctly assigned when comparing between the true generating model and a “simmap” using the inferred transition rates in the HMM. Results are shown for 100 trees each of size of $20,50,80,110,140$, and 170. The true generating model used rates from Equation (5(5) $$\begin{array}{l}{\mathbf{Q}}_{Ind}=\left[\begin{array}{cccc}\ast & 3& 3& 0\\ 2& \ast & 0& 3\\ 2& 0& \ast & 3\\ 0& 2& 2& \ast \end{array}\right],{\mathbf{Q}}_{Dep}=\left[\begin{array}{cccc}\ast & 3& 3& 0\\ 2& \ast & 0& 8h\\ 2& 0& \ast & 8h\\ 0& \frac{1}{2h}& \frac{1}{2h}& \ast \end{array}\right],\end{array}$$(5) ) with h = 3.

In our final simulation, we also considered cases where there were multiple transitions from the independent rate class to the dependent rate class. Here we simulated 100 trees and chose a random number of nodes in the range of [3, 8] that transitioned from the independent rate class to the dependent rate class.

Figures 7 and 8 show that in the case of multiple transitions in the hidden rate class, it is still possible to identify where in the tree these transitions occur with reasonable accuracy. The results have a larger variance and a lower median accuracy compared to the same rates and the same size of tree with only a one-node transition from the independent rate class to the dependent rate class, see Figure 6.

Figure 7. One example of comparing the correct matching edges between the true generating model (a) and the fitting model under the HMM (b). The proportion of matching edges is 82. 9%.

Figure 8. The percentage of the edge-length where the state of dependence or independence was correctly assigned when comparing between the true generating model and a “simmap” using the inferred transition rates in the HMM. Results are shown for 100 trees of a size of 140. The true generating model used rates from Equation (4(4) $$\begin{array}{ll}{\mathbf{Q}}_{Ind}& =\left[\begin{array}{cccc}\ast & 0.04& 0.04& 0\\ 0.03& \ast & 0& 0.04\\ 0.03& 0& \ast & 0.04\\ 0& 0.03& 0.03& \ast \end{array}\right],\\ {\mathbf{Q}}_{Dep}& =\left[\begin{array}{cccc}\ast & 0.04& 0.04& 0\\ 0.03& \ast & 0& 0.4\\ 0.03& 0& \ast & 0.4\\ 0& 0.00375& 0.00375& \ast \end{array}\right].\end{array}$$(4) ). Each true generating tree has between [3, 8] internodes of interest that transition from the independent rate class to the dependent rate class.

Overall, the simulation study showed that it is possible to detect that a model with two rate classes is more appropriate based on AICc scores. However, to detect exactly where these transitions occur the rate of gene gain and loss has to be quite large—perhaps unrealistically so.

3. Neofunctionalization model II

We extend the idea of Model I in Section 2 to the following more general model of neofunctionalization, referred to as Model II. Given n potential genes that may interact with one another, let $$\begin{array}{ll}{\mathcal{S}}^{(n)}& ={\mathcal{S}}_{ind}^{(n)}\cup {\mathcal{S}}_d^{(n)},\\ {\mathcal{S}}_{ind}^{(n)}& =\{(0;i_1,\dots ,i_n):i_k\in \{0,1\},k=1,\dots ,n\},\\ {\mathcal{S}}_d^{(n)}& =\{(1;i_1,\dots ,i_n):i_k\in \{0,1\},k=1,\dots ,n\},\end{array}$$ where $i_k=1,0$ corresponds to the presence or absence of gene k, respectively; while $(1;i_1,\dots ,i_n)$ and $(0;i_1,\dots ,i_n)$ denotes the dependence and independence of the genes.

Assume as before that, given state $s=(0;i_1,\dots ,i_n)$, gene $k=1,\dots ,n$ may be gained or lost at a Poisson rate α_k and β_k, respectively. On the other hand, given a state $s=(1;i_1,\dots ,i_n)$, gene $k=1,\dots ,n$ may be gained or lost at a Poisson rate ${\alpha }_k{\overline{r}}_s$ and β_kr_s, respectively.

Suppose that an additional gene may be added to the functional pathway at some Poisson rate δ > 0. Then, given the current state $(m;i_1,\dots ,i_n)$, a transition to some state $(m;i_1,\dots ,i_n,1)$ with $i_{n+1}=1$ occurs at rate δ, for $m=0,1$.

Then, to model the evolution of genes, we construct a level-dependent Quasi-Birth-and-Death process (LD-QBD) $\{(X(t),\phi (t)):t\ge 0\}$ with state space (6) $$\begin{array}{l}\mathcal{S}=\{(n,\phi ):n=0,1,\dots ;\phi \in {\mathcal{S}}^{(n)}\},\end{array}$$(6) in which the level n is the number of genes that may potentially interact in some functional pathway, and phase φ records the information about the presence or absence of the genes, as well as dependence or independence.

Furthermore, assuming that N is the maximum number of potential genes involved in an interaction, the generator of the process $\{(X(t),\phi (t)):t\ge 0\}$ is matrix $$\begin{array}{ll}\mathbf{Q}& =[{\mathbf{Q}}^{[n,n^{{}^{\prime }}]}{]}_{n,n^{{}^{\prime }}\in \{2,3,\dots ,N\}}\\ & =\left[\begin{array}{cccccc}{\mathbf{Q}}^{[2,2]}& {\mathbf{Q}}^{[2,3]}& \mathbf{0}& \cdots & \cdots & \mathbf{0}\\ \mathbf{0}& {\mathbf{Q}}^{[3,3]}& {\mathbf{Q}}^{[3,4]}& \cdots & \cdots & \mathbf{0}\\ ⋮& ⋮& ⋮& \cdots & \cdots & ⋮\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \cdots & \mathbf{0}& {\mathbf{Q}}^{[N,N]}\end{array}\right],\end{array}$$ with the following expressions for the block matrices in 𝐐.

For $n=2,\dots ,N-1$, the block matrix $$\begin{array}{ll}{\mathbf{Q}}^{[n,n+1]}& =[q_{(m;i_1,\dots ,i_n)(m^{{}^{\prime }};i_1^{{}^{\prime }},\dots ,i_{n+1}^{{}^{\prime }})}]\end{array}$$ is such that when $(m;i_1,\dots ,i_n)=(1;1,\dots ,1)$ then $$\begin{array}{ll}{q}_{(m;i_1,\dots ,i_n)(m;i_1,\dots ,i_n,1)}& =\delta \end{array}$$ and $$\begin{array}{ll}{q}_{(m;i_1,\dots ,i_n)(m^{{}^{\prime }};i_1^{{}^{\prime }},\dots ,i_{n+1}^{{}^{\prime }})}& =0\end{array}$$ otherwise.

For $2,\dots ,N$, the block matrix $$\begin{array}{ll}{\mathbf{Q}}^{[n,n]}& =[q_{(m;i_1,\dots ,i_n)(m^{{}^{\prime }};i_1^{{}^{\prime }},\dots ,i_n^{{}^{\prime }})}]\end{array}$$ is such that $$\begin{array}{ll}{\mathbf{Q}}^{[n,n]}& =[{\mathbf{Q}}_{[\ell ,{\ell }^{\prime }]}^{[n,n]}{]}_{\ell ,{\ell }^{{}^{\prime }}\in \{0,1,\dots ,n\}}\\ & =\left[\begin{array}{cccccc}{\mathbf{Q}}_{[0,0]}^{[n,n]}& {\mathbf{Q}}_{[0,1]}^{[n,n]}& \mathbf{0}& \cdots & \cdots & \mathbf{0}\\ {\mathbf{Q}}_{[1,0]}^{[n,n]}& {\mathbf{Q}}_{[1,1]}^{[n,n]}& {\mathbf{Q}}_{[1,2]}^{[n,n]}& \cdots & \cdots & \mathbf{0}\\ ⋮& ⋮& ⋮& \cdots & \cdots & ⋮\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \cdots & {\mathbf{Q}}_{[n-1,n]}^{[n,n]}& {\mathbf{Q}}_{[n,n]}^{[n,n]}\end{array}\right],\end{array}$$ where the block matrix $$\begin{array}{ll}{\mathbf{Q}}_{[\ell ,{\ell }^{{}^{\prime }}]}^{[n,n]}& =[q_{(m;i_1,\dots ,i_n)(m^{{}^{\prime }};i_1^{{}^{\prime }},\dots ,i_n^{{}^{\prime }})}{]}_{i_1+\dots +i_n=\ell ;i_1^{{}^{\prime }}+\dots +i_n^{{}^{\prime }}={\ell }^{\prime }}\end{array}$$ collects all transition rates from states with ℓ present genes ($i_1+\dots +i_n=\ell$) to states with with ℓ^′ present genes ($i_1^{{}^{\prime }}+\dots +i_n^{{}^{\prime }}={\ell }^{\prime }$).

As a simple example, matrix 𝐀 in (3(3) $$\begin{array}{l}\mathbf{A}=\begin{array}{cc}& \begin{array}{cccccccc}00& 01& 10& 11& \widetilde{00}& \widetilde{01}& \widetilde{10}& \widetilde{11}\end{array}\\ \begin{array}{c}00\\ 01\\ 10\\ 11\\ \widetilde{00}\\ \widetilde{01}\\ \widetilde{10}\\ \widetilde{11}\end{array}& \left[\begin{array}{cccccccc}\ast & {\alpha }_2& {\alpha }_1& 0& 0& 0& 0& 0\\ {\beta }_2& \ast & 0& {\alpha }_1& 0& 0& 0& 0\\ {\beta }_1& 0& \ast & {\alpha }_2& 0& 0& 0& 0\\ 0& {\beta }_1& {\beta }_2& \ast & 0& 0& 0& \lambda \\ \gamma & 0& 0& 0& \ast & {\alpha }_2{r}_{00}& {\alpha }_1{r}_{00}& 0\\ 0& \gamma & 0& 0& {\beta }_2{r}_{01}& \ast & 0& {\alpha }_1{\overline{r}}_{01}\\ 0& 0& \gamma & 0& {\beta }_1{r}_{10}& 0& \ast & {\alpha }_2{\overline{r}}_{10}\\ 0& 0& 0& \gamma & 0& {\beta }_1{r}_{11}& {\beta }_2{r}_{11}& \ast \end{array}\right].\end{array}\end{array}$$(3) ) corresponds to n = 2, and so, when rearranged, is equivalent to matrix 𝐐^{[2, 2]}, given by $$\begin{array}{ll}{\mathbf{Q}}^{[2,2]}& =[{\mathbf{Q}}_{[\ell ,{\ell }^{\prime }]}^{[2,2]}{]}_{\ell ,{\ell }^{{}^{\prime }}\in \{0,1,2\}}\end{array}$$ with $$\begin{array}{ll}{\mathbf{Q}}_{[0,0]}^{[2,2]}& =\left[\begin{array}{cc}\ast & 0\\ \gamma & \ast \end{array}\right],\\ {\mathbf{Q}}_{[0,1]}^{[2,2]}& =\left[\begin{array}{cccc}{\alpha }_2& {\alpha }_1& 0& 0\\ 0& 0& {\alpha }_2{r}_{00}& {\alpha }_1{r}_{00}\end{array}\right],\end{array}$$ and $$\begin{array}{ll}{\mathbf{Q}}_{[1,0]}^{[2,2]}& =\left[\begin{array}{cc}{\beta }_2& 0\\ {\beta }_1& 0\\ 0& {\beta }_2{r}_{01}\\ 0& {\beta }_1{r}_{10}\end{array}\right],\\ {\mathbf{Q}}_{[1,1]}^{[2,2]}& =\left[\begin{array}{cccc}\ast & 0& 0& 0\\ 0& \ast & 0& 0\\ \gamma & 0& \ast & 0\\ 0& \gamma & 0& \ast \end{array}\right],\\ {\mathbf{Q}}_{[1,2]}^{[2,2]}& =\left[\begin{array}{cc}{\alpha }_1& 0\\ {\alpha }_2& 0\\ 0& {\alpha }_1{\overline{r}}_{01}\\ 0& {\alpha }_2{\overline{r}}_{10}\end{array}\right],\end{array}$$ and $$\begin{array}{ll}{\mathbf{Q}}_{[2,1]}^{[2,2]}& =\left[\begin{array}{cccc}{\beta }_1& {\beta }_2& 0& 0\\ 0& 0& {\beta }_1{r}_{11}& {\beta }_2{r}_{11}\end{array}\right],\\ {\mathbf{Q}}_{[2,2]}^{[2,2]}& =\left[\begin{array}{cc}\ast & \lambda \\ \gamma & \ast \end{array}\right].\end{array}$$

That is, $$\begin{array}{ll}{\mathbf{Q}}^{[2,2]}& =\left[\begin{array}{ccc}{\mathbf{Q}}_{[0,0]}^{[2,2]}& {\mathbf{Q}}_{[0,1]}^{[2,2]}& \mathbf{0}\\ {\mathbf{Q}}_{[1,0]}^{[2,2]}& {\mathbf{Q}}_{[1,1]}^{[2,2]}& {\mathbf{Q}}_{[1,2]}^{[2,2]}\\ \mathbf{0}& {\mathbf{Q}}_{[2,1]}^{[2,2]}& {\mathbf{Q}}_{[2,2]}^{[2,2]}\end{array}\right]\\ & =\left[\begin{array}{cccccccc}\ast & 0& {\alpha }_2& {\alpha }_1& 0& 0& 0& 0\\ \gamma & \ast & 0& 0& {\alpha }_2{r}_{00}& {\alpha }_1{r}_{00}& 0& 0\\ {\beta }_2& 0& \ast & 0& 0& 0& {\alpha }_1& 0\\ {\beta }_1& 0& 0& \ast & 0& 0& {\alpha }_2& 0\\ 0& {\beta }_2{r}_{01}& \gamma & 0& \ast & 0& 0& {\alpha }_1{\overline{r}}_{01}\\ 0& {\beta }_1{r}_{10}& 0& \gamma & 0& \ast & 0& {\alpha }_2{\overline{r}}_{10}\\ 0& 0& {\beta }_1& {\beta }_2& 0& 0& \ast & \lambda \\ 0& 0& 0& 0& {\beta }_1{r}_{11}& {\beta }_2{r}_{11}& \gamma & \ast \end{array}\right].\end{array}$$

If the maximum number of potential genes involved in an interaction is N = 2, then we simply have (9) $$\begin{array}{ll}\mathbf{Q}& ={\mathbf{Q}}^{[2,2]},\end{array}$$(9) however, if N > 2, then 𝐐 has the more general form, as in (3), and then 𝐐^{[2, 3]} is an 8×16 matrix, with (10) $$\begin{array}{ll}{\mathbf{Q}}^{[2,3]}& =\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& \delta \end{array}\right],\end{array}$$(10) since the only possible transition in such a matrix is from state $(m;i_1,\dots ,i_n)=(1;1,\dots ,1)$ to state $(m^{{}^{\prime }};i_1^{{}^{\prime }},\dots ,i_{n+1}^{{}^{\prime }})=(1;1,\dots ,1,1)$, which occurs at rate δ.

The expressions for various metrics of interest for such defined LD-QBD can then be obtained by applying results from the literature of matrix-analytic methods, including, e.g., Ramaswami[28], Joyner and Fralix[14], Phung-Duc et al.[27], and Samuelson et al.[31]. Therefore, the analysis of the evolution of gene content and their association across a phylogeny is possible with our model proposed here.

3.1. Model where N genes become associated

We note that the process will eventually be absorbed into the set of states with N associated genes. Therefore, we now consider the generator matrix 𝐐^{[N, N]}, which by 3, is given by (11) $$\begin{array}{ll}{\mathbf{Q}}^{[N,N]}& =[{\mathbf{Q}}_{\ell ,{\ell }^{{}^{\prime }}\in \{0,1,\cdots ,N\}}^{[N,N]}]\\ & =\left[\begin{array}{ccccc}{\mathbf{Q}}_{[0,0]}^{[N,N]}& {\mathbf{Q}}_{[0,1]}^{[N,N]}& \mathbf{0}& \mathbf{0}& \cdots \mathbf{0}\\ {\mathbf{Q}}_{[1,0]}^{[N,N]}& {\mathbf{Q}}_{[1,1]}^{[N,N]}& {\mathbf{Q}}_{[1,2]}^{[N,N]}& \mathbf{0}& \cdots \mathbf{0}\\ \mathbf{0}& {\mathbf{Q}}_{[2,1]}^{[N,N]}& {\mathbf{Q}}_{[2,2]}^{[N,N]}& {\mathbf{Q}}_{[2,3]}^{[N,N]}& \cdots \mathbf{0}\\ ⋮& ⋮⋮& ⋮\cdots & ⋮\\ \mathbf{0}& \mathbf{0}& \cdots \mathbf{0}& {\mathbf{Q}}_{[N,N-1]}^{[N,N]}& {\mathbf{Q}}_{[N,N]}^{[N,N]}\end{array}\right].\end{array}$$(11)

Note that while the structure of Equation (11(11) $$\begin{array}{ll}{\mathbf{Q}}^{[N,N]}& =[{\mathbf{Q}}_{\ell ,{\ell }^{{}^{\prime }}\in \{0,1,\cdots ,N\}}^{[N,N]}]\\ & =\left[\begin{array}{ccccc}{\mathbf{Q}}_{[0,0]}^{[N,N]}& {\mathbf{Q}}_{[0,1]}^{[N,N]}& \mathbf{0}& \mathbf{0}& \cdots \mathbf{0}\\ {\mathbf{Q}}_{[1,0]}^{[N,N]}& {\mathbf{Q}}_{[1,1]}^{[N,N]}& {\mathbf{Q}}_{[1,2]}^{[N,N]}& \mathbf{0}& \cdots \mathbf{0}\\ \mathbf{0}& {\mathbf{Q}}_{[2,1]}^{[N,N]}& {\mathbf{Q}}_{[2,2]}^{[N,N]}& {\mathbf{Q}}_{[2,3]}^{[N,N]}& \cdots \mathbf{0}\\ ⋮& ⋮⋮& ⋮\cdots & ⋮\\ \mathbf{0}& \mathbf{0}& \cdots \mathbf{0}& {\mathbf{Q}}_{[N,N-1]}^{[N,N]}& {\mathbf{Q}}_{[N,N]}^{[N,N]}\end{array}\right].\end{array}$$(11) ) looks almost identical to Equation (3), the matrices 𝐐^{[n, n]} for n < N in Equation (3) correspond to non-conservative generators, with some row sums being strictly negative due to the possibility of transitions from level n to level (n+1). On the other hand, 𝐐^{[N, N]} in Equation (11(11) $$\begin{array}{ll}{\mathbf{Q}}^{[N,N]}& =[{\mathbf{Q}}_{\ell ,{\ell }^{{}^{\prime }}\in \{0,1,\cdots ,N\}}^{[N,N]}]\\ & =\left[\begin{array}{ccccc}{\mathbf{Q}}_{[0,0]}^{[N,N]}& {\mathbf{Q}}_{[0,1]}^{[N,N]}& \mathbf{0}& \mathbf{0}& \cdots \mathbf{0}\\ {\mathbf{Q}}_{[1,0]}^{[N,N]}& {\mathbf{Q}}_{[1,1]}^{[N,N]}& {\mathbf{Q}}_{[1,2]}^{[N,N]}& \mathbf{0}& \cdots \mathbf{0}\\ \mathbf{0}& {\mathbf{Q}}_{[2,1]}^{[N,N]}& {\mathbf{Q}}_{[2,2]}^{[N,N]}& {\mathbf{Q}}_{[2,3]}^{[N,N]}& \cdots \mathbf{0}\\ ⋮& ⋮⋮& ⋮\cdots & ⋮\\ \mathbf{0}& \mathbf{0}& \cdots \mathbf{0}& {\mathbf{Q}}_{[N,N-1]}^{[N,N]}& {\mathbf{Q}}_{[N,N]}^{[N,N]}\end{array}\right].\end{array}$$(11) ) corresponds to a conservative generator with all row sums equal to zero. The process with the generator 𝐐^{[N, N]} can be analyzed using the results from the literature on the LD-QBDs[14, 27, 29]. The stationary distribution of this process, relevant to the analysis here, is considered in the next section.

Denote by ${\alpha }_i,{\beta }_i,i\in \{1,2,\cdots ,N\}$ the rates at which the i-th gene is gained and lost, respectively. As discussed in Section 2, we also consider the parameters that affect the transition rates within the set 𝒮_d, which we denote by $r_{\ell },{\bar{r}}_{\ell },r_N,{\bar{r}}_N>0$. Parameters $r_{\ell },{\bar{r}}_{\ell }$ affect the rate of gene gain and loss when ℓ < N, and we expect that the values of ${\bar{r}}_{\ell },r_{\ell }$ are close to 1. Parameters ${\bar{r}}_N,r_N$ affect the rate of gene gain and loss when ℓ = N, when all the genes are present and perform some useful function, and we expect r_N < 1 and ${\bar{r}}_N>1$.

We now describe the blocks ${\mathbf{Q}}_{[\ell ,{\ell }^{{}^{\prime }}]}^{[N,N]}$ in the generator 𝐐^{[N, N]}. For ℓ = 0, we have, $$\begin{array}{ll}{\mathbf{Q}}_{[0,0]}^{[N,N]}& =\left[\begin{array}{cc}\ast & 0\\ \gamma & \ast \end{array}\right]\end{array}$$ and (12) $$\begin{array}{l}{\mathbf{Q}}_{[0,1]}^{[N,N]}=\left[\begin{array}{cccccc}{\alpha }_N& {\alpha }_{N-1}& \cdots {\alpha }_1& 0& 0& \cdots 0\\ 0& 0& \cdots 0& {\alpha }_N{\bar{r}}_0& {\alpha }_{N-1}{\bar{r}}_0& \cdots {\alpha }_1{\bar{r}}_0\end{array}\right].\end{array}$$(12)

For ℓ = 1, (13) $$\begin{array}{l}{\mathbf{Q}}_{[1,0]}^{[N,N]}=\left[\begin{array}{cc}{\beta }_N& 0\\ ⋮& ⋮\\ {\beta }_1& 0\\ 0& {\beta }_N{r}_1\\ ⋮& ⋮\\ 0& {\beta }_1{r}_1\end{array}\right],\end{array}$$(13) (14) $$\begin{array}{l}{\mathbf{Q}}_{[1,1]}^{[N,N]}=\left[\begin{array}{cc}\ast & \mathbf{O}\\ \mathit{\gamma }& \ast \end{array}\right],\end{array}$$(14) where 𝜸 is a N×N diagonal matrix with diagonal values γ. Next, assuming the order of the rows in ${\mathbf{Q}}_{[1,2]}^{[N,N]}$ according to $(0,\dots ,0,1),(0,\dots ,0,1,0),\dots ,(1,0,\dots ,0)$ and columns in ${\mathbf{Q}}_{[1,2]}^{[N,N]}$ according to $(1,0,\dots ,0,1),(0,1,0\dots ,0,1,0),\dots ,(1,0,\dots ,0,1)$, we can write $$\begin{array}{ll}& {\mathbf{Q}}_{[1,2]}^{[N,N]}\\ & =\left[\begin{array}{cccccccccccccc}{\alpha }_1& {\alpha }_2& \cdots & {\alpha }_{N-1}& 0& 0& 0& 0& 0& \cdots & \cdots & 0& 0& 0\\ 0& 0& \cdots & {\alpha }_N& {\alpha }_{N-2}& \cdots & {\alpha }_1& 0& 0& \cdots & \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \cdots & \cdots & ⋮& ⋮& ⋮& ⋮\\ {\alpha }_N& 0& \cdots & 0& 0& 0& {\alpha }_{n-1}& 0& \cdots & \cdots & \cdots & 0& 0& 0\\ 0& 0& \cdots & 0& 0& 0& 0& {\alpha }_1{\bar{r}}_1& {\alpha }_2{\bar{r}}_1& \cdots & {\alpha }_{N-1}{\bar{r}}_1& 0& \cdots & 0\\ 0& 0& \cdots & 0& 0& 0& 0& 0& 0& \cdots & {\alpha }_N{\bar{r}}_1& {\alpha }_{N-2}{\bar{r}}_1& \cdots & {\alpha }_1{\bar{r}}_1\\ ⋮& ⋮& \cdots & ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \cdots & 0& 0& 0& 0& {\alpha }_N{\bar{r}}_1& 0& \cdots & 0& 0& \cdots & {\alpha }_{N-1}{\bar{r}}_1\end{array}\right].\end{array}$$

For $\ell =1,2,\cdots ,N-1$, we notice that the number of possible states of $i_1+\cdots +i_n=\ell$ is that total number of combinations equal to $(\genfrac{}{}{0ex}{}{N}{\ell })$, and so the number of possible states to transition to such that $i_1^{\prime }+\cdots +i_n^{\prime }={\ell }^{\prime }$ where ${\ell }^{\prime }=\ell \pm 1$, is equal to $(\genfrac{}{}{0ex}{}{N}{{\ell }^{\prime }})$.

We apply similar methods to write blocks ${\mathbf{Q}}_{[\ell ,{\ell }^{{}^{\prime }}]}^{[N,N]}$, ${\ell }^{{}^{\prime }}\in \{\ell -1,\ell ,\ell +1\}$ for $\ell =2,\dots ,N-1$. For ℓ = N, we have (16) $$\begin{array}{l}{\mathbf{Q}}_{[N,N-1]}^{[N,N]}=\left[\begin{array}{ccccc}{\beta }_1& {\beta }_2& \cdots {\beta }_N& 0& \cdots 0\\ 0& 0& \cdots 0& {\beta }_1{r}_N& \cdots {\beta }_N{r}_N\end{array}\right],\end{array}$$(16) (17) $$\begin{array}{l}{\mathbf{Q}}_{[N-1,N]}^{[N,N]}=\left[\begin{array}{cc}{\alpha }_1& 0\\ ⋮& ⋮\\ {\alpha }_N& 0\\ 0& {\alpha }_1{\bar{r}}_N\\ ⋮& ⋮\\ 0& {\alpha }_N{\bar{r}}_N\end{array}\right],\end{array}$$(17) and (18) $$\begin{array}{l}{\mathbf{Q}}_{[N,N]}^{[N,N]}=\left[\begin{array}{cc}\ast & \lambda \\ \gamma & \ast \end{array}\right].\end{array}$$(18)

Assume that the gene gain or loss for each gene is equally likely. The rates of gaining or losing a gene are equal to α_i = α and β_i = β for all $i\in \{1,2,\cdots ,N\}$.

3.2. Numerical results

Define the stationary distribution vector $\mathit{\pi }={\mathit{\pi }}_{\ell =0,1,\cdots ,N}$, ${\mathit{\pi }}_{\ell }=[{\pi }_{(\ell ;i_1,i_2\cdots ,i_N)}{]}_{i_1+i_2+\cdots +i_N=\ell }$, given by (19) $$\begin{array}{l}{\pi }_{({\ell }^{{}^{\prime }};i_1^{{}^{\prime }},i_2^{{}^{\prime }}\cdots ,i_N^{{}^{\prime }})}=\underset{t\to \mathrm{\infty }}{lim}\mathcal{P}((X(t),\phi (t))=(\ell {}^{\prime };i_1^{{}^{\prime }},i_2^{{}^{\prime }}\cdots ,i_N^{{}^{\prime }})),\end{array}$$(19) where the limits exist and do not depend on the initial state because of the stability of the process.

In order to compute 𝝅_ℓ, for $\ell =0,1,\cdots ,N$, we apply an iterative scheme in Baumann and Sandmann[5].

1. Compute ${\mathbf{R}}_{N-1}=-{\mathbf{Q}}_{[N-1,N]}{\mathbf{Q}}_{[N,N]}^{-1}$;

2. Iterate ${\mathbf{R}}_{\ell -1}={\mathbf{Q}}_{[\ell -1,\ell ]}({\mathbf{Q}}_{[\ell ,\ell ]}+{\mathbf{R}}_{\ell }{\mathbf{Q}}_{[\ell +1,\ell ]}{)}^{-1}$ for $\ell =N-1,N-2,\cdots ,1$;

3. Compute a solution 𝝅₀≠𝟎 of ${\mathit{\pi }}_0({\mathbf{Q}}_{[0,0]}+{\mathbf{R}}_0{\mathbf{Q}}_{[1,0]})=\mathbf{0}$;

4. Iterate ${\mathit{\pi }}_{\ell +1}={\mathit{\pi }}_{\ell }{\mathbf{R}}_{\ell }$ for $\ell =0,1,\cdots ,N-1$;

5. Renormailse 𝝅_ℓ using ${\mathit{\pi }}_{\ell }=\frac{{\mathit{\pi }}_{\ell }}{\sum _{\ell =0}^N{\mathit{\pi }}_{\ell }}$.

We present examples of stationary distribution when N = 4. We apply the value in Equation (5(5) $$\begin{array}{l}{\mathbf{Q}}_{Ind}=\left[\begin{array}{cccc}\ast & 3& 3& 0\\ 2& \ast & 0& 3\\ 2& 0& \ast & 3\\ 0& 2& 2& \ast \end{array}\right],{\mathbf{Q}}_{Dep}=\left[\begin{array}{cccc}\ast & 3& 3& 0\\ 2& \ast & 0& 8h\\ 2& 0& \ast & 8h\\ 0& \frac{1}{2h}& \frac{1}{2h}& \ast \end{array}\right],\end{array}$$(5) ) when h = 1 in Figure 9.

Figure 9. The stationary distribution π_ℓ with $\ell =0,1,\cdots ,4$. Here the γ = 1, is the transition rate from the dependent rate class to the independent rate class. In (a), λ = 4 is the transition rate from the independent rate class to the dependent rate class, and in (b), λ = 2.

Now we calculate the stationary distribution using the value in Equation (4(4) $$\begin{array}{ll}{\mathbf{Q}}_{Ind}& =\left[\begin{array}{cccc}\ast & 0.04& 0.04& 0\\ 0.03& \ast & 0& 0.04\\ 0.03& 0& \ast & 0.04\\ 0& 0.03& 0.03& \ast \end{array}\right],\\ {\mathbf{Q}}_{Dep}& =\left[\begin{array}{cccc}\ast & 0.04& 0.04& 0\\ 0.03& \ast & 0& 0.4\\ 0.03& 0& \ast & 0.4\\ 0& 0.00375& 0.00375& \ast \end{array}\right].\end{array}$$(4) ) in Figure 10. Since the processes have smaller rates of gene gain and loss, the transition rates from the independent rate class to the dependent rate class and from the dependent rate class to the independent rate class are both smaller.

Figure 10. The stationary distribution π_ℓ with $\ell =0,1,\cdots ,4$. Here the γ = 0. 025, is the transition rate from the dependent rate class to the independent rate class. In (a), λ = 0. 1 is the transition rate from the independent rate class to the dependent rate class, and in (b), λ = 0. 2.

4. Discussion

In this article, we have developed a new approach for modeling neofunctionalization across an evolutionary tree. In previous work by Diao et al.[9, 10] neofunctionalization has been modeled in terms of a single gene acquiring a new beneficial function, here, we instead think of it in terms of groups of two or more genes gaining an association, i.e., becoming linked in a biochemical pathway that performs some new function. We applied the package corHMM from Boyko and Beaulieu[7] to generate gene presence-absence data under the model and then assessed how well we were able to detect if a HMM was the best model, and if so, where on the tree transitions between the dependent and independent rate classes occurred. While we were often able to reject simpler one-rate-class models in favor of the true HMM, we were not, in general, able to get accurate knowledge of where in the tree the transitions between rate classes occurred unless rates of gene gain and loss were high. We observed a large variance in the percentage of edges where we could correctly infer the hidden rate class.

The corHMM package allows users to specify which rates should be equal, but it does not allow us to impose additional constraints on whether some rates should be larger or smaller than others. In the biological scenario, we are modeling it makes sense that rates of loss from the 11 state should be lower in 𝐐_Dep than 𝐐_Ind. It is possible that an inference algorithm that was able to enforce this constraint would achieve higher levels of accuracy in determining where transitions between dependence and independence occur.

In our simulations with a single transition from independence to dependence, we found that using larger trees and larger mean rates of gain and loss had better accuracy at identifying the hidden rate classes along the edges, in the best cases, the median accuracy was as high as 95%. However, using a lower mean rate, which Boyko and Beaulieu[7] mentioned as leading to better performance of corHMM, gave us poorer results. Additionally, we conducted some simulations where there were multiple transitions from the independent rate class to the dependent rate class. Those simulations led to a larger variance and less accuracy, but still had a median accuracy of around 90%.

To extend to the general case where N genes become associated, we developed a level-dependent quasi-birth-and-death (LD-QBD) model. The level, n, records the current number of genes that are potentially required to perform some joint function, and the phase records which genes are present or absent. In our model, all n genes must be present for selection acting on some useful function to offer protection (reduce the rates of gene loss). The genes act like components of a series circuit, where if any one component is missing, the beneficial function is lost. One possible extension of the work is to consider more flexible configurations where there is some redundancy.

Pan-genome studies are becoming increasingly prevalent, so in the future, there will be many more data sets where gene presence-absence information is available for numerous species, e.g., for crops and their wild relatives[6]. For this reason, it will be important to continue to develop stochastic models that can advance our understanding of biological processes. The models presented here are a useful step in this direction.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

We would like to thank the Australian Research Council for funding this research through Discovery Project DP180100352. This work was partly funded by The Australian Research Council Centre of Excellence for Plant Success in Nature and Agriculture (CE200100015).