BACKGROUND
a gene family is a group of genes with similar sequences that show evidence of descent from a common ancestor .
this includes orthologs that originate through speciation as well as duplicates  that can be found within a species or shared between species from an older duplication event that predated speciation.
the large number of genes per family suggests that the newly arisen gene duplicates are potentially major contributors to evolutionary novelties .
gene duplication can provide raw genetic material for evolutionary forces to act on.
although a majority of duplicate genes may be silenced by degenerative mutations or lost due to population dynamics, some duplicated genes are able to evolve novel functions permanently preserved in the population  <cit> .
accurately estimating the timing and mode of gene duplications along the evolutionary history of species can provide invaluable information about underlying mechanisms by which the genomes of organisms evolved and the genes with novel functions arose  <cit> .
several biological models have been proposed to depict the mechanisms that lead to different evolutionary fates of a gene duplicate .
nonfunctionalization refers to the process in which mutations occur on one of the gene duplicates and produce a non-functional protein  <cit> .
the neofunctionalization model  <cit>  assumes that duplication itself does not affect fitness.
although a duplicate is most likely to be pseudogenized by degenerative mutation  or lost due to population dynamics  <cit> , the redundant copy may occasionally acquire a new beneficial function through mutation that will be preferentially preserved in the population.
while this function may subsequently be optimized and accommodated within the genome structure  by an evolutionary stokes shift  <cit> , the initial event leading to retention is a single beneficial change.
the waiting time for this single change gives rise to a convexly decaying hazard function when modeled together with non-functionalizing changes and is referred to as the neofunctionalization model .
the duplication-degeneration-complementation model  <cit>  describes a so-called subfunctionalization mechanism in which two gene copies are partially damaged by degenerated mutations.
both copies must be maintained in order to perform the original function of the gene  <cit> .
this model, called subfunctionalization, involves a waiting time for multiple events to occur as deleterious substitutions accumulate in both copies before the retaining mutation can occur.
this waiting time for multiple changes gives rise to a switch from a convex to a concave  hazard function when modeled together with non-functionalizing mutations .
in addition to the processes acting on individual genes, large-scale gene duplication events  may have occurred and produced multiple interacting genes together creating an additional retention mechanism.
dosage balance promotes the retention of duplicated interaction networks, as loss of interaction stoichiometry can lead to declines in fitness.
this gives rise to very different retention dynamics compared to neofunctionalization or subfunctionalization .
the mechanistic models described for nonfunctionalization, subfunctionalization, neofunctionalization, and dosage balance represent one of many conceivable modeling frameworks for duplicate gene retention .
the models here are used within a single population, reflecting a lineage of a phylogenetic tree, but the ultimate aim is to extend their use into an interspecific phylogenetic framework with the population genetic assumptions that accompany this.
simpler models have already been incorporated into a fuller phylogenetic framework of this nature .
accurately reconstructing the evolution of gene families requires informative datasets, powerful mathematical models, and efficient computational algorithms.
advanced biotechnologies provide a vast amount of genetic data for understanding the evolution of gene families  <cit> .
meanwhile, probabilistic models, describing the process of gene family evolution, significantly enhance our ability to extract useful information from genetic data .
the birth-death  model  <cit> , which has been broadly applied in analyzing species phylogenies  <cit> , could also be adopted in phylogenetic analysis of gene families  <cit> .
in  <dig>  thompson  <cit>  introduced a phylogenetic model based on the birth-death process to understand the evolution of human populations.
under the generalized birth-death model, nee et al.
<cit>  derived a reconstructed evolutionary process  <cit>  to estimate birth and death rates in a interspecific phylogenetic framework.
rannala and yang  <cit>  developed a birth-death phylogenetic model for estimating phylogenetic trees from molecular sequence data.
aldous and popovic  <cit>  proposed a continuous-time critical branching process conditioned on the number of species in the present, with the assumption that the birth and death rates are identical in macroevolution, which was later relaxed by gernhard  <cit>  to allow uncorrelated birth and death rates.
with the assumption of constant birth and death rates, stadler  <cit>  derived the probability density function of a phylogenetic tree under the birth-death model.
recently, time-dependent bd processes have attracted more attention as a mode of performing hypothesis-driven research .
rabosky  <cit>  distinguished rate-variable models of diversification from rate-constant models by fitting bd models using likelihood methods.
hohna  <cit>  and hallinan  <cit>  studied the reconstructed process with time-dependent rates in a more general setting by relaxing the assumptions about the number of species and the time of the process.
the bd model was first adopted in  <cit>  and further extended by other researchers to reconcile gene and species trees .
recently, boussau et al.
<cit>  established a bd phylogenetic model for co-estimating gene and species trees without the need of estimation of divergence times in species trees and duplication and loss rates.
the current computational methods for analyzing gene family data  suffer a variety of weakness that need to be addressed.
there is no probabilistic model embedded in a birth-death phylogenetic modeling framework that can take into account the effects of different evolutionary mechanisms of gene retention when analyzing gene family data.
it is desirable to build a stochastic model as a good approximation to the real biological process of gene duplication and loss.
such probabilistic models can both add biological realism to improve the fit of the model to the data as well as enable mechanistic inference that is currently not possible.
in this study, we integrate several evolutionary mechanisms of gene retention into the age-dependent bd model , in which the loss rate is a function of the ages of gene copies.
moreover, we derive the probability density function of gene duplication times for each mechanism.
the conditional density function of a duplication time given the previous duplication time is derived from the reconstructed process under the generalized birth-death model  <cit> .
the conditional density function can be utilized to calculate the joint density of duplication times, and to efficiently simulate duplication times under the generalized bd model.
the simulation results suggest that the maximum likelihood approach can accurately estimate the parameters in the generalized bd model for different mechanisms of gene retention, and the proposed gene-retention model can be used to detect the underlying mechanism that drives the evolutionary process of duplicates within a gene family.
methods
modeling the loss rate
for simplicity, we consider the process of gene duplication/loss in a single population.
for a single population, we assume that a gene copy may duplicate or die at time t. the homogeneous birth-death model assumes that the rate of loss  of a duplicated gene is constant through time  <cit> .
this expectation is consistent with the nonfunctionalization process, but does not take into account any of the processes of neofunctionalization and subfunctionalization, which can affect the loss rate of gene duplicates.
the birth-death model for the fates of gene retention  includes a time-dependent loss rate and a constant duplication rate λ. the time-dependent loss rates will be extended to age-dependent loss rates in the age-dependent birth-death model .
the process starts at time  <dig>  and the number of gene copies at time  <dig> is  <dig>  the process of gene duplication and loss occurs under the following postulates  <cit> :  the probability that a duplication will occur during an infinitesimal interval , while the probability that no duplication will occur is 1- ntλΔt + o, and  the probability that a gene duplicate will be lost during an infinitesimal interval , while the probability that no loss will occur is 1- ntμt Δt + o, in which the loss rate μt is a function of time t.
we introduce three formulas for the loss rate μt based on the processes of nonfunctionalization, neofunctionalization, and subfunctionalization, with assumptions about these processes made in the introduction and also described in  <cit> .
for nonfunctionalization, the loss rate μt is constant over time t, i.e., μt = μ. the neofunctionalization hazard rate  declines with time  <cit> .
averaging across the probability of hitting a neofunctionalizing substitution, the nonfunctionalization probability for duplicate genes declines, leading to the overall decline of duplicate loss over long evolutionary time periods  <cit> .
this convexly declining loss rate has been described with a weibull hazard function to characterize the average process   <cit> .
we use an exponential function to model the loss rate of neofunctionlization, i.e., μt = αe− tα for 0 < α <  <dig>  further, the subfunctionalization loss rate behavior has been characterized to be concavely  declining based upon theoretical expectations of a waiting time for complementary mutations  <cit> .
konrad  <cit>  introduced an extended exponential hazard function to describe the instantaneous rate of loss.
we adopt a generalized logistic function for the loss rate μt of subfunctionalization, i.e., μt=αeγ−t1+eγ−t, in which the scale parameter 0 < α <  <dig> and known location parameter γ >  <dig>
the time-dependent birth-death model
we are interested in the probability distribution of duplication times of the reconstructed lineages , because the phylogeny reconstructed from the sequences of contemporary species does not include the extinct lineages  <cit> .
the pure birth process of the reconstructed lineages can be derived from a generalized birth-death process  <cit> .
we use the following notations which are defined closely to nee et al.
<cit>  throughout this paper.
let t2 =  <dig> be the first duplication time at the root of the tree, and t be the present time .
let nt be the number of lineages at the present time t. let ni be the number of reconstructed lineages alive at ti that survive to the present.
we use {ti | i =  <dig>  …, nt} to denote the duplication times of nt lineages at the tips of a phylogenetic tree, and t2
< t3
< t4
< … < t. let p be the probability that one lineage at time τ leaves multiple descendants at the present time t, i.e., p = p ,  <dig> pτt=1+∫τtμteρτtdt− <dig>
in eq.
, ρ = ∫τtds.
since the integral ∫τtμteρ in eq.
is analytically intractable, it is approximated by a monte carlo method.
we define uij as the probability p that one lineage at time ti leaves multiple descendant reconstructed lineages at a later time tj.
this probability has been derived under the time-dependent bd model, i.e., uij=pnj>1|ni=1=1−ptitjeρtitj  in  <cit> ).
given the number nt of lineages at the present time t and the number n <dig> of lineages at time  <dig>  the probability density function of the duplication times t = {ti | i = n0 +  <dig>  …, nt} is given by  <cit>   <dig> ft|nt,n <dig> t=∏i=n0+1nti−1λptit1−ηti− <dig> tii−1nt−1n0−11−η <dig> tn0η <dig> tnt−n <dig>
in , ηij=1−1−uit1−ujt.
the conditional probability distribution of duplication time ti , given its previous duplication time ti- <dig>  t and nt, is given by  <cit>   <dig> fti|ti− <dig> nt,t=fti|ti−1pnt|nti,tpnt|nti− <dig> t.
in eq.
, fti|ti−1=i−1λptit1−ηti− <dig> tii− <dig>  and  in  <cit> ).
with the conditional densities f of duplication times, the duplication events between times  <dig> and t can be simulated recursively in forward direction.
the conditional density in  differs from the density of duplication times derived by hohna  <cit> , in which the duplication events are treated as a random sample from a common probability distribution.
the age-dependent birth-death model
the time-dependent birth-death model described in the previous section starts with a single population corresponding to the lineage of a phylogenetic tree and assumes a molecular clock that starts ticking for all duplicates at the root.
thus, in the time-dependent birth-death model, the loss rate μt of a gene copy is a function of time t. however, the loss rate μt should be a function of the ages of gene copies.
in this section, the time-dependent birth-death process is extended to the age-dependent process, where the clock for each duplicate starts ticking at its birth.
when the loss rate is constant , the age-dependent model is identical with the time-dependent model.
thus, we only describe the age-dependent model for neofunctionalizaiton and subfunctionalization.
in the age-dependent model, the expressions for the loss rates of neofunctionalization and subfunctionalization remain unchanged , except that time t is replaced with the age t’ of the gene copy, i.e., μt'=αe−t'α for neofunctionalization and μt'=αeγ−t'1+eγ−t' for subfunctionalization.
moreover, it is assumed that the number of gene copies increases or decreases by  <dig> or remains the same during an infinitesimal interval  4a pnt+Δt=nt+1=ntλΔt+οΔt  4b pnt+Δt=nt−1=∑i=1ntμti'Δt+οΔt  4c pnt+Δt=nt=1−ntλ+∑i=1ntμti'Δt+οΔt
in , μti' is the loss rate of gene copy i at the age of ti' for i =  <dig>   <dig>  …, nt.
let ti <dig> be the duplication time of gene copy i. the age ti' of gene copy i is a random variable, because it is a function of the random duplication time ti <dig>  i.e., ti' = t − ti <dig>  therefore,  and  are integrated over all possible values of μti' with respect to the probability density function f of the age t ' of a gene copy.
the age-dependent loss rate μti' in  and  is replaced with its expectation eμti'.
since all ti' s have the same probability distribution, the loss rates of nt gene copies have the same expected values.
let t <dig> be the most recent duplication time of a gene copy that survives to time t. since t <dig> is the most recent duplication time, it indicates that no duplication or loss events have occurred between t <dig> and t on the gene copy.
it has been shown that the number of duplication or loss events follows the poisson distribution with mean ∫0tdx.
the probability of no duplication or loss events occurring within the time interval  is equal to e−∫0tλ+μxdx.
thus, the probability density of duplication time t <dig> is proportional to dt0e−∫0tλ+μxdx for 0 < t0 < t, in which dt <dig> is the duplication rate at time t <dig> and e−∫0tλ+μxdx is the probability that t <dig> is the most recent duplication time of the gene copy.
given that duplication occurs on a specific lineage, dt <dig> is equal to the duplication rate λ. thus, the probably density of the most recent duplication time t <dig> is  <dig> ft0=e−∫t0tλ+μxdx∫0te−∫t0tλ+μxdxdt <dig>
because the gene age t’ is equal to t – t <dig>  the probability density of age t’ for 0 < t0 < t is given by  <dig> ft'=e−∫t−t'tλ+μxdx∫0te−∫t−t'tλ+μxdxdt'
since the denominator in  is intractable, it is approximated by monte carlo simulation.
it follows that the mean loss rate at time t is ϕt=eμti'=∫0tμt'ft'dt'.
thus, the postulates in  and  become p = nϕtΔt + ο and p = 1 − ntΔt + ο. the loss rate in eq.
is replaced by the mean loss rate ϕt accordingly and p is modified as  <dig> pτt=1+∫τtϕteρτtdt− <dig>
finally, the joint and conditional probability density of duplication times  for the age dependent model remain unchanged, except that the loss rate μt in eq.
is replaced with the mean loss rate ϕt.
RESULTS
simulation for the time-dependent model
to evaluate the performance of the time-dependent birth-death model on simulated data where the true values of parameters are known, we generated duplication times of gene copies using the rejection-sampling algorithm with the conditional probability density function of duplication times in eq.
.
we found the maximum likelihood score for the conditional probability distribution using an optimization function optim in r. the maximum score was used as the upper bound in the rejection-sampling algorithm.
specifically, duplication times were simulated from eq.
with a fixed current time t =  <dig> and a fixed number of gene copies nt =  <dig> at time t. the first duplication time is set to  <dig>  i.e., t2 = 0; the second one is simulated conditional on the first one and so on so that additional  <dig> duplication times are generated sequentially.
duplication events were generated under each of  <dig> duplication mechanisms  with different parameterizations specified in table  <dig>  we set a constant duplication rate λ =  <dig>  for all simulations .
the loss rates were determined by the equations described previously for nonfunctionalization, neofunctionalization, and subfunctionalization models with parameters shown in table  <dig>  the values of parameters were selected such that three mechanisms have the same initial deletion rate.table  <dig> the values of parameters used in simulating duplication times under nonfunctionalization, neofunctionalization, and subfunctionalization are shown
for each mechanism, simulation was repeated  <dig> times.
the mean of simulated duplication times for each of three mechanisms are shown in fig.
1a.
duplication times simulated under different mechanisms show distinct patterns.
given the present time t and the number of gene copies nt, the overall duplication times for nonfunctionalization tend to be larger than those for neofunctionalization and subfunctionalization, and duplication times for neofunctionalization appear to be smaller than subfunctionalization.
the curves of duplication times for nonfunctionalization, neofunctionalization, and subfunctionalization are well separated , even though three mechanisms have the same duplication rate and the same starting deletion rate.
these results indicate that duplication times can be used to distinguish different mechanisms of gene retention, and to make inference about the underlying mechanism that generated the observed duplication times given the assumptions of the duplication models and their relationship to the underlying biology.
these results are consistent with the caveat that the time-dependent process uses a tree-dependent clock rather than the more biological situation of a duplication-event specific process.
the extension to the age-dependent birth-death model will be discussed below.
the joint probability density function in eq.
can be used to obtain the maximum likelihood estimates  of parameters in the time-dependent model, when duplication times are given as input data.
to visualize the divergence of the probability density functions of three mechanisms, we plotted the density curves of the first duplication time for nonfunctionalization, neofunctionalization, and subfunctionalization  with the values of parameters in table  <dig>  since each mechanism has a unique density curve, this result indicates that it is possible to distinguish three mechanisms using the time-dependent birth-death model.
moreover, we employed the akaike information criterion   <cit>  to evaluate the relative quality of the time-dependent models for nonfunctionalization, neofunctionalization, and subfunctionalization.
the data sets simulated from the time-dependent model were used as input data to calculate aic for nonfunctionalization, neofunctionalization, and subfunctionalization.
for each simulated data set, the mechanism with the lowest aic score was selected and compared with the true mechanism from which the data sets were generated.
we reported the percentage of the simulated data sets successfully identifying the true mechanism .
the overall average of the percentages of samples recovering the true mechanism is about 80 % .
in addition, subfunctionalization appears to be more difficult than neofunctionalization to distinguish from nonfunctionalization in this modeling framework .fig.
<dig> simulation results of the time-dependent model:  the means of duplication times simulated with  <dig> replicates for nonfunctionalization, neofunctionalization, and subfuncitonalization are shown;  the probability density curves of duplication times for nonfunctionalization, neofunctionalization, and subfunctionalization under the model are shown;  the percentage of samples identifying the true mechanism with aic
to examine the performance of maximum likelihood estimation, we use the simulated duplication times as data to estimate model parameters.
the sample size  ranges from  <dig> to  <dig>  the maximum likelihood estimates of parameters were obtained using metropolis-hastings markov chain monte carlo algorithm.
the standard errors of the maximum likelihood estimates are displayed in fig.
<dig>  for nonfunctionalization, the standard errors of the estimates of μ and λ decrease as the number of duplication times increases from  <dig> to  <dig>  similarly, the standard errors of the estimates of parameters for subfunctionalization and neofunctionalization decrease as the number of duplications grows.
however, the estimation of parameter α for neofunctionalization does not improve well with the increased number of gene copies , because duplication times in the simulated data are highly correlated and the auto-correlation between two adjacent duplication times increases as the number of duplication times increases.
as a result, when the number of highly correlated duplication times reaches a certain number, adding even more duplication times does not contribute more information for accurately estimating model parameters, especially for neofunctionalization where the loss rate quickly declines to a very low level.
similar results about biases and parameter estimates under constant and time-dependent birth-death processes have been obtained in  <cit> .
nevertheless, these results suggest that maximum likelihood methods can accurately estimate parameters in the time-dependent birth-death model when the sample size is large.fig.
<dig> the standard errors of the maximum likelihood estimates of parameters in the age-dependent models for nonfunctionalization, neofunctionliazation, and subfunctionalization
simulation for the age-dependent birth-death model
the simulation for the age-dependent model was conducted with the same parameterization and simulation procedure used for the time-dependent model.
we generated duplication times from the age-dependent models for nonfunctionalization, neofunctionalization, and subfunctionalization.
the mean duplication times given the current time t and gene copy number nt for the age-dependent models  appear to be less dispersed among nonfunctionalization, neofunctionalization, and subfunctionalization than those for the time-dependent models .
in addition, the density curve for subfunctionalization becomes closer to the nonfunctionalization curve under the age-dependent model , compared to the curves for the time-dependent model .
this is consistent with our expectation, because the age of a gene copy is less than the absolute time t and the beginning portion of the concavely declining loss rate of subfunctionalization is similar to the constant rate of nonfunctionalization.
in fig.
3b, the density curve for neofunctionalization is well separated from the density curves for nonfunctionalization and subfunctionalization.
in contrast, the loss rate of subfunctionalization is assumed to be a backwards-s-shaped logistic function of time, which is an intermediate state between the loss rates of nonfunctionalization and neofunctionalization.
if the loss curve of subfunctionalization moves to the right, it becomes closer to nonfunctionalization .
conversely, when the loss rate curve moves to the left, it gets closer to neofunctionalization .
although subfunctionalization is an intermediate state between nonfunctionalization and neofunctionalization, it is expected to be more similar to neofunctionalization, which can be tested in real data analysis.
the ultimate similarity comes with increasing time, as both neofunctionalization and subfunctionalization culminate in reduced hazard rates, unlike nonfunctionalization.
with a fixed duplication rate, these processes are expected to result in an increased number of copies.
conditional on the number of copies, subfunctionalization and neofunctionalization would be consistent with a reduced duplication rate and older duplication times.
the overall percentage of samples identifying the true mechanism increases as the number of gene copies grows .
the percentages of nonfunctionalization and neofunctionalization are significantly higher than the overall percentage.
although the performance of subfunctionalization is below average, the percentage of samples successfully identifying the true subfunctionalization increases to 60 % when the number of gene copies reaches  <dig>  moreover, the standard errors of the estimates of parameters in the age-dependent model appear to decrease as the number of gene copies grows, suggesting that maximum likelihood methods can accurately estimate parameters in the age-dependent model, when the sample size is large .fig.
<dig> simulation results of the age-dependent model:  the means of duplication times simulated with  <dig> replicates for nonfunctionalization, neofunctionalization, and subfuncitonalization are shown;  the probability density curves of duplication times for nonfunctionalization, neofunctionalization, and subfunctionalization under the model are shown;  the percentage of samples identifying the true mechanism with aic
fig.
<dig> the standard errors of maximum likelihood estimates of parameters in the age-dependent models for neofunctionliazation and subfunctionalization
discussion
summary of the gene family evolution model
we have derived the probability density function for the age-dependent birth-death model, in which the loss rate is a function of the ages of gene copies.
in addition, the conditional density function and a joint density function of duplication times with age-dependent loss rate have been developed in above age-dependent model, given the current time t and the number of gene copies at the time t. the conditional density function is used to efficiently simulate duplication times, and the simulation results suggest that maximum likelihood methods can accurately estimate model parameters in both time-dependent and age-dependent models.
in addition, the relative quality of various birth-death models was assessed with aic.
both time-dependent and age-dependent models can distinguish the three mechanisms  with high probabilities when the sample size is large.
these results indicate that the probabilistic models derived from the birth-death process with a time-dependent and age-dependent loss rates are useful for understanding the duplication and loss mechanisms of gene families that evolve over time in a single population with caveats discussed.
limitations and future study
as duplication times are often not observable, it is desirable to generalize the current model to dna sequence data.
we are currently working along this line to build a generalized model that includes two stochastic processes.
the birth and death process is used to derive the probability distribution of a gene family tree, while the mutation process is used to derive the probability distribution of dna sequence data given the gene family tree.
with this generalized model, we can estimate model parameters  from dna sequence data.
one of the limits of the current model is that it considers gene family evolution in a single population.
this model cannot be applied as currently implemented to understand the evolutionary process of gene families from multiple species.
to overcome this limit, the current model will be extended in the context of species trees, in which duplication process occurs along the lineages of species trees.
this generalization will certainly involve intensive computation, but such a model is quite useful for understanding gene family evolution in the context of the evolution of species.
another limitation of the current age-dependent model is that the likelihood is conditioned on observed extant duplicate copies and does not consider the full generative process including duplicates that were lost before the present.
future work will examine this in the context of approximate bayesian computation  <cit> .
further, the current model exists in the classes of interspecific models that treat all observations from a single individual from a species as fixed relative to observations from single individuals from other species.
recently, a correction for the effects of population dynamics has been introduced and can be considered in modeling efforts  <cit> .
missing data and genome assembly error are also not specifically addressed in the modeling framework and their impact on inference also needs to be addressed.
the gene loss models and their interpretations  make assumptions about the relationship between the accumulation of synonymous changes and of non-synonymous changes whereas there is information in the evolution of dn/ds vs. ds that can be taken advantage of in alternative formulations of the likelihood .
lastly, the models can be used to make predictions about functional evolution in the absence of actual functional data.
while such data does not currently exist in large scale, the future may bring data on the expression levels of protein duplicates compared to an ancestral state as well as binding and enzyme specificities , all of which can be integrated into a phylogenetic framework.
however, even with future comparative proteomic data, one still needs models that treat signals associated with selective pressures , as neutral changes in expression and functional properties would not lead to changes in retention profiles  and meaningful lineage-specific biology .
the model as currently developed also assumes that all duplicates in a gene family evolve under the same process.
a future opportunity is in examination of large gene family databases like ensembl  <cit> , hogenom  <cit> , or taed  <cit> , a mixture model of duplicate processes can be applied across all gene families and duplication events to enable a posteriori probabilistic identification of duplication retention mechanisms for individual gene duplication events.
the work presented in this manuscript, with a birth-death model in a phylogenetic context, brings this scale of modeling one step closer.
CONCLUSIONS
we develop a generalized birth-death model for probabilistic inference of the evolutionary mechanism for duplicate gene retention based upon the average rate of loss over time of the duplicate.
the time-dependent birth-death model assumes a molecular clock that starts ticking for all duplicates at the root.
the time-dependent model is then extended to the age-dependent model, which allows the gene loss rate dependent of duplication events.
simulation results indicate that the mechanisms of gene retentions  produce distinct likelihood functions, which can be used with comparative genomic data to quantitatively distinguish those mechanisms.
availability of supporting data
this study of a theoretical nature has not generated any novel supporting data.
abbreviations
bdbirth-death
mlemaximum likelihood estimate
competing interests
the authors declare that they have no competing interests.
authors’ contributions
the study was conceived by dal and ll.
mathematical derivations were performed by jz and ll.
simulations and programming were performed by jz and ait.
the manuscript was written by jz, ait, dal, and ll.
all authors have read and approved the final version of the manuscript.
this research is supported by national science foundation .

 a gene family is a group of genes with similar sequences that show evidence of descent from a common ancestor. this includes orthologs that originate through speciation as well as duplicates that can be found within a species or shared between species from an older duplication event that predated speciation. the large number of genes per family suggests that the newly arisen gene duplicates are potentially major contributors to evolutionary novelties. gene duplication can provide raw genetic material for evolutionary forces to act on. although a majority of duplicate genes may be silenced by degenerative mutations or lost due to population dynamics, some duplicated genes are able to evolve novel functions permanently preserved in the population. accurately estimating the timing and mode of gene duplications along the evolutionary history of species can provide invaluable information about underlying mechanisms by which the genomes of organisms evolved and the genes with novel functions arose. several biological models have been proposed to depict the mechanisms that lead to different evolutionary fates of a gene duplicate. nonfunctionalization refers to the process in which mutations occur on one of the gene duplicates and produce a non-functional protein. the neofunctionalization model assumes that duplication itself does not affect fitness. although a duplicate is most likely to be pseudogenized by degenerative mutation or lost due to population dynamics, the redundant copy may occasionally acquire a new beneficial function through mutation that will be preferentially preserved in the population. while this function may subsequently be optimized and accommodated within the genome structure by an evolutionary stokes shift, the initial event leading to retention is a single beneficial change. the waiting time for this single change gives rise to a convexly decaying hazard function when modeled together with non-functionalizing changes and is referred to as the neofunctionalization model. the duplication-degeneration-complementation model describes a so-called subfunctionalization mechanism in which two gene copies are partially damaged by degenerated mutations. both copies must be maintained in order to perform the original function of the gene. this model, called subfunctionalization, involves a waiting time for multiple events to occur as deleterious substitutions accumulate in both copies before the retaining mutation can occur. this waiting time for multiple changes gives rise to a switch from a convex to a concave hazard function when modeled together with non-functionalizing mutations. in addition to the processes acting on individual genes, large-scale gene duplication events may have occurred and produced multiple interacting genes together creating an additional retention mechanism. dosage balance promotes the retention of duplicated interaction networks, as loss of interaction stoichiometry can lead to declines in fitness. this gives rise to very different retention dynamics compared to neofunctionalization or subfunctionalization. the mechanistic models described for nonfunctionalization, subfunctionalization, neofunctionalization, and dosage balance represent one of many conceivable modeling frameworks for duplicate gene retention. the models here are used within a single population, reflecting a lineage of a phylogenetic tree, but the ultimate aim is to extend their use into an interspecific phylogenetic framework with the population genetic assumptions that accompany this. simpler models have already been incorporated into a fuller phylogenetic framework of this nature. accurately reconstructing the evolution of gene families requires informative datasets, powerful mathematical models, and efficient computational algorithms. advanced biotechnologies provide a vast amount of genetic data for understanding the evolution of gene families. meanwhile, probabilistic models, describing the process of gene family evolution, significantly enhance our ability to extract useful information from genetic data. the birth-death model, which has been broadly applied in analyzing species phylogenies, could also be adopted in phylogenetic analysis of gene families. in thompson introduced a phylogenetic model based on the birth-death process to understand the evolution of human populations. under the generalized birth-death model, nee et al. derived a reconstructed evolutionary process to estimate birth and death rates in a interspecific phylogenetic framework. rannala and yang developed a birth-death phylogenetic model for estimating phylogenetic trees from molecular sequence data. aldous and popovic proposed a continuous-time critical branching process conditioned on the number of species in the present, with the assumption that the birth and death rates are identical in macroevolution, which was later relaxed by gernhard to allow uncorrelated birth and death rates. with the assumption of constant birth and death rates, stadler derived the probability density function of a phylogenetic tree under the birth-death model. recently, time-dependent bd processes have attracted more attention as a mode of performing hypothesis-driven research. rabosky distinguished rate-variable models of diversification from rate-constant models by fitting bd models using likelihood methods. hohna and hallinan studied the reconstructed process with time-dependent rates in a more general setting by relaxing the assumptions about the number of species and the time of the process. the bd model was first adopted in and further extended by other researchers to reconcile gene and species trees. recently, boussau et al. established a bd phylogenetic model for co-estimating gene and species trees without the need of estimation of divergence times in species trees and duplication and loss rates. the current computational methods for analyzing gene family data suffer a variety of weakness that need to be addressed. there is no probabilistic model embedded in a birth-death phylogenetic modeling framework that can take into account the effects of different evolutionary mechanisms of gene retention when analyzing gene family data. it is desirable to build a stochastic model as a good approximation to the real biological process of gene duplication and loss. such probabilistic models can both add biological realism to improve the fit of the model to the data as well as enable mechanistic inference that is currently not possible. in this study, we integrate several evolutionary mechanisms of gene retention into the age-dependent bd model, in which the loss rate is a function of the ages of gene copies. moreover, we derive the probability density function of gene duplication times for each mechanism. the conditional density function of a duplication time given the previous duplication time is derived from the reconstructed process under the generalized birth-death model. the conditional density function can be utilized to calculate the joint density of duplication times, and to efficiently simulate duplication times under the generalized bd model. the simulation results suggest that the maximum likelihood approach can accurately estimate the parameters in the generalized bd model for different mechanisms of gene retention, and the proposed gene-retention model can be used to detect the underlying mechanism that drives the evolutionary process of duplicates within a gene family. methods modeling the loss rate for simplicity, we consider the process of gene duplication/loss in a single population. for a single population, we assume that a gene copy may duplicate or die at time t. the homogeneous birth-death model assumes that the rate of loss of a duplicated gene is constant through time. this expectation is consistent with the nonfunctionalization process, but does not take into account any of the processes of neofunctionalization and subfunctionalization, which can affect the loss rate of gene duplicates. the birth-death model for the fates of gene retention includes a time-dependent loss rate and a constant duplication rate λ. the time-dependent loss rates will be extended to age-dependent loss rates in the age-dependent birth-death model. the process starts at time and the number of gene copies at time is the process of gene duplication and loss occurs under the following postulates: the probability that a duplication will occur during an infinitesimal interval, while the probability that no duplication will occur is 1- ntλΔt + o, and the probability that a gene duplicate will be lost during an infinitesimal interval, while the probability that no loss will occur is 1- ntμt Δt + o, in which the loss rate μt is a function of time t. we introduce three formulas for the loss rate μt based on the processes of nonfunctionalization, neofunctionalization, and subfunctionalization, with assumptions about these processes made in the introduction and also described in. for nonfunctionalization, the loss rate μt is constant over time t, i.e., μt = μ. the neofunctionalization hazard rate declines with time. averaging across the probability of hitting a neofunctionalizing substitution, the nonfunctionalization probability for duplicate genes declines, leading to the overall decline of duplicate loss over long evolutionary time periods. this convexly declining loss rate has been described with a weibull hazard function to characterize the average process. we use an exponential function to model the loss rate of neofunctionlization, i.e., μt = αe− tα for 0 < α <  further, the subfunctionalization loss rate behavior has been characterized to be concavely declining based upon theoretical expectations of a waiting time for complementary mutations. konrad introduced an extended exponential hazard function to describe the instantaneous rate of loss. we adopt a generalized logistic function for the loss rate μt of subfunctionalization, i.e., μt=αeγ−t1+eγ−t, in which the scale parameter 0 < α <  and known location parameter γ >  the time-dependent birth-death model we are interested in the probability distribution of duplication times of the reconstructed lineages, because the phylogeny reconstructed from the sequences of contemporary species does not include the extinct lineages. the pure birth process of the reconstructed lineages can be derived from a generalized birth-death process. we use the following notations which are defined closely to nee et al. throughout this paper. let t2 =  be the first duplication time at the root of the tree, and t be the present time. let nt be the number of lineages at the present time t. let ni be the number of reconstructed lineages alive at ti that survive to the present. we use {ti | i =  …, nt} to denote the duplication times of nt lineages at the tips of a phylogenetic tree, and t2  < t3  < t4  < … < t. let p be the probability that one lineage at time τ leaves multiple descendants at the present time t, i.e., p = p, pτt=1+∫τtμteρτtdt− in eq., ρ = ∫τtds. since the integral ∫τtμteρ in eq. is analytically intractable, it is approximated by a monte carlo method. we define uij as the probability p that one lineage at time ti leaves multiple descendant reconstructed lineages at a later time tj. this probability has been derived under the time-dependent bd model, i.e., uij=pnj>1|ni=1=1−ptitjeρtitj in ). given the number nt of lineages at the present time t and the number n of lineages at time the probability density function of the duplication times t = {ti | i = n0 +  …, nt} is given by ft|nt,n t=∏i=n0+1nti−1λptit1−ηti− tii−1nt−1n0−11−η tn0η tnt−n in, ηij=1−1−uit1−ujt. the conditional probability distribution of duplication time ti, given its previous duplication time ti- t and nt, is given by fti|ti− nt,t=fti|ti−1pnt|nti,tpnt|nti− t. in eq., fti|ti−1=i−1λptit1−ηti− tii− and in ). with the conditional densities f of duplication times, the duplication events between times and t can be simulated recursively in forward direction. the conditional density in differs from the density of duplication times derived by hohna, in which the duplication events are treated as a random sample from a common probability distribution. the age-dependent birth-death model the time-dependent birth-death model described in the previous section starts with a single population corresponding to the lineage of a phylogenetic tree and assumes a molecular clock that starts ticking for all duplicates at the root. thus, in the time-dependent birth-death model, the loss rate μt of a gene copy is a function of time t. however, the loss rate μt should be a function of the ages of gene copies. in this section, the time-dependent birth-death process is extended to the age-dependent process, where the clock for each duplicate starts ticking at its birth. when the loss rate is constant, the age-dependent model is identical with the time-dependent model. thus, we only describe the age-dependent model for neofunctionalizaiton and subfunctionalization. in the age-dependent model, the expressions for the loss rates of neofunctionalization and subfunctionalization remain unchanged, except that time t is replaced with the age t’ of the gene copy, i.e., μt'=αe−t'α for neofunctionalization and μt'=αeγ−t'1+eγ−t' for subfunctionalization. moreover, it is assumed that the number of gene copies increases or decreases by or remains the same during an infinitesimal interval 4a pnt+Δt=nt+1=ntλΔt+οΔt 4b pnt+Δt=nt−1=∑i=1ntμti'Δt+οΔt 4c pnt+Δt=nt=1−ntλ+∑i=1ntμti'Δt+οΔt in, μti' is the loss rate of gene copy i at the age of ti' for i =  …, nt. let ti be the duplication time of gene copy i. the age ti' of gene copy i is a random variable, because it is a function of the random duplication time ti i.e., ti' = t − ti therefore, and are integrated over all possible values of μti' with respect to the probability density function f of the age t' of a gene copy. the age-dependent loss rate μti' in and is replaced with its expectation eμti'. since all ti' s have the same probability distribution, the loss rates of nt gene copies have the same expected values. let t be the most recent duplication time of a gene copy that survives to time t. since t is the most recent duplication time, it indicates that no duplication or loss events have occurred between t and t on the gene copy. it has been shown that the number of duplication or loss events follows the poisson distribution with mean ∫0tdx. the probability of no duplication or loss events occurring within the time interval is equal to e−∫0tλ+μxdx. thus, the probability density of duplication time t is proportional to dt0e−∫0tλ+μxdx for 0 < t0 < t, in which dt is the duplication rate at time t and e−∫0tλ+μxdx is the probability that t is the most recent duplication time of the gene copy. given that duplication occurs on a specific lineage, dt is equal to the duplication rate λ. thus, the probably density of the most recent duplication time t is ft0=e−∫t0tλ+μxdx∫0te−∫t0tλ+μxdxdt because the gene age t’ is equal to t – t the probability density of age t’ for 0 < t0 < t is given by ft'=e−∫t−t'tλ+μxdx∫0te−∫t−t'tλ+μxdxdt' since the denominator in is intractable, it is approximated by monte carlo simulation. it follows that the mean loss rate at time t is ϕt=eμti'=∫0tμt'ft'dt'. thus, the postulates in and become p = nϕtΔt + ο and p = 1 − ntΔt + ο. the loss rate in eq. is replaced by the mean loss rate ϕt accordingly and p is modified as pτt=1+∫τtϕteρτtdt− finally, the joint and conditional probability density of duplication times for the age dependent model remain unchanged, except that the loss rate μt in eq. is replaced with the mean loss rate ϕt.  simulation for the time-dependent model to evaluate the performance of the time-dependent birth-death model on simulated data where the true values of parameters are known, we generated duplication times of gene copies using the rejection-sampling algorithm with the conditional probability density function of duplication times in eq.. we found the maximum likelihood score for the conditional probability distribution using an optimization function optim in r. the maximum score was used as the upper bound in the rejection-sampling algorithm. specifically, duplication times were simulated from eq. with a fixed current time t =  and a fixed number of gene copies nt =  at time t. the first duplication time is set to i.e., t2 = 0; the second one is simulated conditional on the first one and so on so that additional duplication times are generated sequentially. duplication events were generated under each of duplication mechanisms with different parameterizations specified in table  we set a constant duplication rate λ =  for all simulations. the loss rates were determined by the equations described previously for nonfunctionalization, neofunctionalization, and subfunctionalization models with parameters shown in table  the values of parameters were selected such that three mechanisms have the same initial deletion rate.table the values of parameters used in simulating duplication times under nonfunctionalization, neofunctionalization, and subfunctionalization are shown for each mechanism, simulation was repeated times. the mean of simulated duplication times for each of three mechanisms are shown in fig. 1a. duplication times simulated under different mechanisms show distinct patterns. given the present time t and the number of gene copies nt, the overall duplication times for nonfunctionalization tend to be larger than those for neofunctionalization and subfunctionalization, and duplication times for neofunctionalization appear to be smaller than subfunctionalization. the curves of duplication times for nonfunctionalization, neofunctionalization, and subfunctionalization are well separated, even though three mechanisms have the same duplication rate and the same starting deletion rate. these results indicate that duplication times can be used to distinguish different mechanisms of gene retention, and to make inference about the underlying mechanism that generated the observed duplication times given the assumptions of the duplication models and their relationship to the underlying biology. these results are consistent with the caveat that the time-dependent process uses a tree-dependent clock rather than the more biological situation of a duplication-event specific process. the extension to the age-dependent birth-death model will be discussed below. the joint probability density function in eq. can be used to obtain the maximum likelihood estimates of parameters in the time-dependent model, when duplication times are given as input data. to visualize the divergence of the probability density functions of three mechanisms, we plotted the density curves of the first duplication time for nonfunctionalization, neofunctionalization, and subfunctionalization with the values of parameters in table  since each mechanism has a unique density curve, this result indicates that it is possible to distinguish three mechanisms using the time-dependent birth-death model. moreover, we employed the akaike information criterion to evaluate the relative quality of the time-dependent models for nonfunctionalization, neofunctionalization, and subfunctionalization. the data sets simulated from the time-dependent model were used as input data to calculate aic for nonfunctionalization, neofunctionalization, and subfunctionalization. for each simulated data set, the mechanism with the lowest aic score was selected and compared with the true mechanism from which the data sets were generated. we reported the percentage of the simulated data sets successfully identifying the true mechanism. the overall average of the percentages of samples recovering the true mechanism is about 80 %. in addition, subfunctionalization appears to be more difficult than neofunctionalization to distinguish from nonfunctionalization in this modeling framework.fig. simulation results of the time-dependent model: the means of duplication times simulated with replicates for nonfunctionalization, neofunctionalization, and subfuncitonalization are shown; the probability density curves of duplication times for nonfunctionalization, neofunctionalization, and subfunctionalization under the model are shown; the percentage of samples identifying the true mechanism with aic to examine the performance of maximum likelihood estimation, we use the simulated duplication times as data to estimate model parameters. the sample size ranges from to the maximum likelihood estimates of parameters were obtained using metropolis-hastings markov chain monte carlo algorithm. the standard errors of the maximum likelihood estimates are displayed in fig.  for nonfunctionalization, the standard errors of the estimates of μ and λ decrease as the number of duplication times increases from to similarly, the standard errors of the estimates of parameters for subfunctionalization and neofunctionalization decrease as the number of duplications grows. however, the estimation of parameter α for neofunctionalization does not improve well with the increased number of gene copies, because duplication times in the simulated data are highly correlated and the auto-correlation between two adjacent duplication times increases as the number of duplication times increases. as a result, when the number of highly correlated duplication times reaches a certain number, adding even more duplication times does not contribute more information for accurately estimating model parameters, especially for neofunctionalization where the loss rate quickly declines to a very low level. similar results about biases and parameter estimates under constant and time-dependent birth-death processes have been obtained in. nevertheless, these results suggest that maximum likelihood methods can accurately estimate parameters in the time-dependent birth-death model when the sample size is large.fig. the standard errors of the maximum likelihood estimates of parameters in the age-dependent models for nonfunctionalization, neofunctionliazation, and subfunctionalization simulation for the age-dependent birth-death model the simulation for the age-dependent model was conducted with the same parameterization and simulation procedure used for the time-dependent model. we generated duplication times from the age-dependent models for nonfunctionalization, neofunctionalization, and subfunctionalization. the mean duplication times given the current time t and gene copy number nt for the age-dependent models appear to be less dispersed among nonfunctionalization, neofunctionalization, and subfunctionalization than those for the time-dependent models. in addition, the density curve for subfunctionalization becomes closer to the nonfunctionalization curve under the age-dependent model, compared to the curves for the time-dependent model. this is consistent with our expectation, because the age of a gene copy is less than the absolute time t and the beginning portion of the concavely declining loss rate of subfunctionalization is similar to the constant rate of nonfunctionalization. in fig. 3b, the density curve for neofunctionalization is well separated from the density curves for nonfunctionalization and subfunctionalization. in contrast, the loss rate of subfunctionalization is assumed to be a backwards-s-shaped logistic function of time, which is an intermediate state between the loss rates of nonfunctionalization and neofunctionalization. if the loss curve of subfunctionalization moves to the right, it becomes closer to nonfunctionalization. conversely, when the loss rate curve moves to the left, it gets closer to neofunctionalization. although subfunctionalization is an intermediate state between nonfunctionalization and neofunctionalization, it is expected to be more similar to neofunctionalization, which can be tested in real data analysis. the ultimate similarity comes with increasing time, as both neofunctionalization and subfunctionalization culminate in reduced hazard rates, unlike nonfunctionalization. with a fixed duplication rate, these processes are expected to result in an increased number of copies. conditional on the number of copies, subfunctionalization and neofunctionalization would be consistent with a reduced duplication rate and older duplication times. the overall percentage of samples identifying the true mechanism increases as the number of gene copies grows. the percentages of nonfunctionalization and neofunctionalization are significantly higher than the overall percentage. although the performance of subfunctionalization is below average, the percentage of samples successfully identifying the true subfunctionalization increases to 60 % when the number of gene copies reaches moreover, the standard errors of the estimates of parameters in the age-dependent model appear to decrease as the number of gene copies grows, suggesting that maximum likelihood methods can accurately estimate parameters in the age-dependent model, when the sample size is large.fig. simulation results of the age-dependent model: the means of duplication times simulated with replicates for nonfunctionalization, neofunctionalization, and subfuncitonalization are shown; the probability density curves of duplication times for nonfunctionalization, neofunctionalization, and subfunctionalization under the model are shown; the percentage of samples identifying the true mechanism with aic fig. the standard errors of maximum likelihood estimates of parameters in the age-dependent models for neofunctionliazation and subfunctionalization discussion summary of the gene family evolution model we have derived the probability density function for the age-dependent birth-death model, in which the loss rate is a function of the ages of gene copies. in addition, the conditional density function and a joint density function of duplication times with age-dependent loss rate have been developed in above age-dependent model, given the current time t and the number of gene copies at the time t. the conditional density function is used to efficiently simulate duplication times, and the simulation results suggest that maximum likelihood methods can accurately estimate model parameters in both time-dependent and age-dependent models. in addition, the relative quality of various birth-death models was assessed with aic. both time-dependent and age-dependent models can distinguish the three mechanisms with high probabilities when the sample size is large. these results indicate that the probabilistic models derived from the birth-death process with a time-dependent and age-dependent loss rates are useful for understanding the duplication and loss mechanisms of gene families that evolve over time in a single population with caveats discussed. limitations and future study as duplication times are often not observable, it is desirable to generalize the current model to dna sequence data. we are currently working along this line to build a generalized model that includes two stochastic processes. the birth and death process is used to derive the probability distribution of a gene family tree, while the mutation process is used to derive the probability distribution of dna sequence data given the gene family tree. with this generalized model, we can estimate model parameters from dna sequence data. one of the limits of the current model is that it considers gene family evolution in a single population. this model cannot be applied as currently implemented to understand the evolutionary process of gene families from multiple species. to overcome this limit, the current model will be extended in the context of species trees, in which duplication process occurs along the lineages of species trees. this generalization will certainly involve intensive computation, but such a model is quite useful for understanding gene family evolution in the context of the evolution of species. another limitation of the current age-dependent model is that the likelihood is conditioned on observed extant duplicate copies and does not consider the full generative process including duplicates that were lost before the present. future work will examine this in the context of approximate bayesian computation. further, the current model exists in the classes of interspecific models that treat all observations from a single individual from a species as fixed relative to observations from single individuals from other species. recently, a correction for the effects of population dynamics has been introduced and can be considered in modeling efforts. missing data and genome assembly error are also not specifically addressed in the modeling framework and their impact on inference also needs to be addressed. the gene loss models and their interpretations make assumptions about the relationship between the accumulation of synonymous changes and of non-synonymous changes whereas there is information in the evolution of dn/ds vs. ds that can be taken advantage of in alternative formulations of the likelihood. lastly, the models can be used to make predictions about functional evolution in the absence of actual functional data. while such data does not currently exist in large scale, the future may bring data on the expression levels of protein duplicates compared to an ancestral state as well as binding and enzyme specificities, all of which can be integrated into a phylogenetic framework. however, even with future comparative proteomic data, one still needs models that treat signals associated with selective pressures, as neutral changes in expression and functional properties would not lead to changes in retention profiles and meaningful lineage-specific biology. the model as currently developed also assumes that all duplicates in a gene family evolve under the same process. a future opportunity is in examination of large gene family databases like ensembl, hogenom, or taed, a mixture model of duplicate processes can be applied across all gene families and duplication events to enable a posteriori probabilistic identification of duplication retention mechanisms for individual gene duplication events. the work presented in this manuscript, with a birth-death model in a phylogenetic context, brings this scale of modeling one step closer.  we develop a generalized birth-death model for probabilistic inference of the evolutionary mechanism for duplicate gene retention based upon the average rate of loss over time of the duplicate. the time-dependent birth-death model assumes a molecular clock that starts ticking for all duplicates at the root. the time-dependent model is then extended to the age-dependent model, which allows the gene loss rate dependent of duplication events. simulation results indicate that the mechanisms of gene retentions produce distinct likelihood functions, which can be used with comparative genomic data to quantitatively distinguish those mechanisms. availability of supporting data this study of a theoretical nature has not generated any novel supporting data. abbreviations bdbirth-death mlemaximum likelihood estimate competing interests the authors declare that they have no competing interests. authors’ contributions the study was conceived by dal and ll. mathematical derivations were performed by jz and ll. simulations and programming were performed by jz and ait. the manuscript was written by jz, ait, dal, and ll. all authors have read and approved the final version of the manuscript. this research is supported by national science foundation.
starting with a single population corresponding to the lineage of a phylogenetic tree and with an assumption of a clock that starts ticking for each duplicate at its birth, an age-dependent birth-death process is developed by extending the results from the time-dependent birth-death process. simulation results indicate that different mechanisms of gene retentions produce distinct likelihood functions, which can be used with genomic data to quantitatively distinguish those mechanisms. we develop an age-dependent birth-death model for understanding the mechanisms of gene retention, which allows a gene loss rate dependent on each duplication event. mechanistic models have previously been introduced that allow for probabilistic inference of the evolutionary mechanism for duplicate gene retention based upon the average rate of loss over time of the duplicate. use of mechanistic models in a phylogenetic framework requires an age-dependent birth-death process. the implementation of such models in a full phylogenetic framework is expected to enable large scale probabilistic analysis of duplicates in comparative genomic studies. in this study, we describe a generalized birth-death process for modeling the fates of gene duplication. gene duplicationphylogenetic methodsprobabilistic modelsbirth-death processesstochastic processeshttp://dx.doi.org/ /100000001national science foundation dms-1222745dms-1222940liberles david a. liu liang issue-copyright-statement© the author 2015 accurately estimating the timing and mode of gene duplications along the evolutionary history of species can provide invaluable information about underlying mechanisms by which the genomes of organisms evolved and the genes with novel functions arose. however, there is currently no probabilistic model embedded in a birth-death modeling framework that can take into account the effects of different evolutionary mechanisms of gene retention when analyzing gene family data.
 accurately estimating the timing and mode of gene duplications along the evolutionary history of species can provide invaluable information about underlying mechanisms by which the genomes of organisms evolved and the genes with novel functions arose. mechanistic models have previously been introduced that allow for probabilistic inference of the evolutionary mechanism for duplicate gene retention based upon the average rate of loss over time of the duplicate. however, there is currently no probabilistic model embedded in a birth-death modeling framework that can take into account the effects of different evolutionary mechanisms of gene retention when analyzing gene family data.  in this study, we describe a generalized birth-death process for modeling the fates of gene duplication. use of mechanistic models in a phylogenetic framework requires an age-dependent birth-death process. starting with a single population corresponding to the lineage of a phylogenetic tree and with an assumption of a clock that starts ticking for each duplicate at its birth, an age-dependent birth-death process is developed by extending the results from the time-dependent birth-death process. the implementation of such models in a full phylogenetic framework is expected to enable large scale probabilistic analysis of duplicates in comparative genomic studies.  we develop an age-dependent birth-death model for understanding the mechanisms of gene retention, which allows a gene loss rate dependent on each duplication event. simulation results indicate that different mechanisms of gene retentions produce distinct likelihood functions, which can be used with genomic data to quantitatively distinguish those mechanisms. keywords gene duplicationphylogenetic methodsprobabilistic modelsbirth-death processesstochastic processeshttp://dx.doi.org/ /100000001national science foundation dms-1222745dms-1222940liberles david a. liu liang issue-copyright-statement© the author 2015 