A model of compound heterozygous, loss-of-function alleles is broadly consistent with observations from complex-disease GWAS datasets

TL;DR: It is demonstrated that the model of gene action, relating genotype to phenotype, has a qualitative effect on several relevant aspects of the population genetic architecture of a complex trait, and the genetic model impacts genetic variance component partitioning across the allele frequency spectrum and the power of statistical tests.

Abstract: The genetic component of complex disease risk in humans remains largely unexplained. A corollary is that the allelic spectrum of genetic variants contributing to complex disease risk is unknown. Theoretical models that relate population genetic processes to the maintenance of genetic variation for quantitative traits may suggest profitable avenues for future experimental design. Here we use forward simulation to model a genomic region evolving under a balance between recurrent deleterious mutation and Gaussian stabilizing selection. We consider multiple genetic and demographic models, and several different methods for identifying genomic regions harboring variants associated with complex disease risk. We demonstrate that the model of gene action, relating genotype to phenotype, has a qualitative effect on several relevant aspects of the population genetic architecture of a complex trait. In particular, the genetic model impacts genetic variance component partitioning across the allele frequency spectrum and the power of statistical tests. Models with partial recessivity closely match the minor allele frequency distribution of significant hits from empirical genome-wide association studies without requiring homozygous effect-sizes to be small. We highlight a particular gene-based model of incomplete recessivity that is appealing from first principles. Under that model, deleterious mutations in a genomic region partially fail to complement one another. This model of gene-based recessivity predicts the empirically observed inconsistency between twin and SNP based estimated of dominance heritability. Furthermore, this model predicts considerable levels of unexplained variance associated with intralocus epistasis. Our results suggest a need for improved statistical tools for region based genetic association and heritability estimation.

Summary

Introduction

• Risk for complex diseases in humans, such as diabetes and hypertension, is highly heritable yet the causal DNA sequence variants responsible for that risk remain largely unknown.
• The authors have declared that no competing interests exist.
• This model assigns an effect size to a mutant allele, but formally makes no concrete statement regarding the molecular nature of the allele.
• When applied to molecular data, such as SNP genotypes in a GWAS, these models treat the SNPs themselves as the loci of interest.
• Further, the model of gene-based recessivity best explains the differences between estimates of additive and dominance variance components from SNP-based methods [27] and from twin studies [28] and is consistent with the distribution of frequencies of significant associations in GWAS [4, 26].

Results and Discussion

• The models As in [36],the authors simulate a 100 kilobase region of human genome, contributing to a complex disease phenotype and fitness.
• The expected fitness effect of a mutation is always deleterious because trait effects are sampled from an exponential distribution.
• Specifically, the authors studied three different genetic models and two different demographic models, holding the fitness model as a constant.
• Parameters are briefly described in Table 1.

Parameter Description

• This reflects the competing forces of increasing average genetic effect and decreasing average allele frequency which occurs as λ increases (S5 Fig).
• Yet, when simulating iMR model, the authors find that an intermediate degree of dominance, h = 0.25, results in distribution of significant hits which is similar to the GBR results (Fig 4).
• The reason the authors emphasize this feature of the data is to demonstrate that models with rare alleles of large effect do not necessarily imply a visual excess of rare significant GWAS hits.
• Of the models the authors explored, only the gene-based recessive model with intermediate to large effects is consistent with the difference between twin and SNP based estimates of dominance variance (Fig 2).

Materials and Methods

• Forward simulation Using the fwdpp template library v0.2.8 [87], the authors implemented a forward in time individualbased simulation of a Wright-Fisher population with mutation under the infinitely many sites model [88], recombination, and selection occurring each generation.
• For comparison, the authors calculated the true total heritability in the sample as H2sample ¼ ðVG;sampleÞ=ðVP;sampleÞ.
• Under MSGREMLd, many replicates resulted in numerical errors in GCTA.
• For dizygotic (DZ) twins the authors used two child gamete pairs, each with a unique environmental deviate.

S10 Fig. Distribution of significant hits under site based recessive models with incomplete

• Horizontal violin plots depict the distribution of minor allele frequencies (MAF) of the most strongly associated single marker in a GWAS.
• Moments were calculated using the boost C++ statistical accumulators library.
• Data are plotted as the mean across model replicates ± the standard error of the mean.
• Shown are the additive co-dominant (AC), gene-based (GBR) and complete multiplicative recessive (Mult. recessive (h = 0); cMR) models.
• (TIFF) S14 Fig. Regression Based Estimates of Genetic Variance.

S17 Fig. Additive genetic variance explained over allele frequency under the Tennessen

• Et al. [40] model for European demography.
• The left column of panels shows how VG changes over time under this model.
• The burden ratio [91] is calculated as the ratio of genetic load between simulations with only ancient growth and those with an additional recent bottleneck and growth.
• Here load is calculated as the average deviation from optimum fitness due to (left) fixed mutations, segregating mutations and all mutations.
• For large effect size models, under which there are relatively more mutations that experience strong selection, the authors see the characteristic drop in the burden ratio following the bottleneck and rebound following re-expansion [91].

S21 Fig. Non-parametric comparison between empirical and simulated GWAS hits. A

• Non-parametric comparison between distribution of allele frequencies between simulated and empirical GWAS hits.
• In cases where more than one marker was tied for the lowest p-value, one was chosen at random.
• Specific information regarding the empirical data can be obtained in S1 Table.
• The dashed lines show the analytical result and the solid curves are empirical cumulative distribution functions based on a sample of 500 mutation effects from an exponential distribution.

Acknowledgments

• The authors are thankful to Joseph Farran, Harry Mangalam, Adam Brenner, Garr Updegraff, and Edward Xia for administering the University of California, Irvine High Performance Computing cluster.
• The authors are thankful to Kirk Lohmueller for helpful detailed comments throughout this project.
• The authors would like to thank Peter Andolfatto, Bogdan Pasaniuc and Nick Mancuso for helpful discussion.

Figures

Fig 1. Variance explained over allele frequency. The cumulative additive and dominance genetic variance which can be explained by markers whose frequencies, q, are x. Each color represents a different value of λ: the mean effects size of a new deleterious mutation. Shown here are the gene-based (GBR), additive co-dominant (AC), incomplete multiplicative recessive (Mult. recessive (h = 0.25); iMR) and complete multiplicative recessive (Mult. recessive (h = 0);cMR) models. Solid lines show the additive variance alone and dotted lines show the combined additive and dominance variance. All data shown are for models where H2 * 0.08. These particular results are robust to changes H2 when VG is not changed, as is the case here. The additive and dominance genetic variance is estimated by the adjusted r2 of the regression of all markers (and their corresponding dominance encoding) with q x onto total genotypic value (see methods for details); data are displayed as the mean of 250 simulation replicates. The vertical dotted and dashed lines correspond to the q = 0.001 and q = 0.05, respectively. The curves under no growth appear to be truncated with respect to rapid growth because the range of the x-axis differs between growth and no growth (minimum q = 1/2N).

Fig 2. Heritability estimates compared to population heritability. Heritability estimates and population heritability as a function of λ: the mean effect size of a new deleterious mutation. Additive (A; orange) component of true heritability is calculated by multiplying the end point(q = 1) of the variance curves in Fig 1 by the broad-sense heritability values summarized in S1 Fig. HE-regression and GREMLd estimates were obtained from random population samples (n = 6000). GREMLd analysis was performed in GCTA using genotype data that was either unfiltered or filtered to remove variants with MAF<0.01. Twin study estimates are directly calculated using MZ and DZ twin correlations from 64 sets of twin studies. Each study consisted of pooling 2000 MZ twin pairs and 2000 DZ twin pairs from each of 8 model replicates for a total of 64,000 individual phenotypes. Data are plotted as the median across replicate sets ± half the interquartile range. Shown are the additive co-dominant (AC), gene-based (GBR) incomplete multiplicative recessive (Mult. recessive (h = 0.25); iMR) and complete multiplicative recessive (Mult. recessive (h = 0); cMR) models.

Fig 4. Distribution of significant GWAS hits. Horizontal violin plots depict the distribution of minor allele frequencies (MAF) of the most strongly associated single marker in a GWAS. Individual hits are plotted as translucent points and jittered to provide a sense of the total number and density of hits. Each panel contains simulated data pooled across model replicates for each value of λ with empirical data for comparison. Empirical data are described in Materials and Methods. In cases where more than one marker was tied for the lowest p-value, one was chosen at random. Shown here are the additive co-dominant (AC), gene-based (GBR), incomplete multiplicative recessive (Mult. recessive (h = 0.25); iMR) and complete multiplicative recessive (Mult. recessive (h = 0);cMR) models. All data shown are for models where H2 * 0.08, because single marker test power was too low under H2 * 0.04 to make informative density plots. To further increase the number simulated data points, we perfromed n = 1,250 replicates at each level for this figure. Simulated data were subjected to ascertainment sampling such that the MAF distribution of all markers on the simulated genotyping chip was uniform. Specific information regarding the empirical data can be obtained in S1 Table.

Table 1. Description of parameters used in the models.

Fig 3. Power of association tests. (A) The power of a single marker logistic regression, at significance threshold of α 10−8, as a function of λ: the mean effect size of a new deleterious mutation. For single marker tests we define power as the number of simulation replicates in which any single marker reaches genome wide significance. Two study designs were emulated. For the gene chip design only markers with MAF > 0.05 were considered and all markers were considered for the resequencing design. Genetic models shown here are the additive co-dominant (AC), gene-based (GBR), complete multiplicative site-based recessive (Mult. recessive (h = 0); cMR) and incomplete multiplicative site-based recessive models (Mult. recessive (h = 0.25); iMR) (B) The power of region-based rare variant association tests to detect association with the simulated causal gene region at significance threshold of α 10−6. For region-based tests, we define power as the percent of simulation replicates in which the p-value of the test was less than α. The p-values for the ESM, c-Alpha were evaluated using 2 × 106 permutations. SKAT p-values were determined by the SKAT R package and represent numerical approximations to the presumed analytical p-value.

Full text

RESEARCH ARTICLE
A Model of Compound Heterozygous, Loss-of-
Function Alleles Is Broadly Consistent with
Observations from Complex-Disease GWAS
Datasets
Jaleal S. Sanjak
1,2
*, Anthony D. Long
1,2
, Kevin R. Thornton
1,2
*
1 Department of Ecology and Evolutionary Biology, University of California, Irvine, Irvine, California, USA,
2 Center for Complex Biological Systems, University of California, Irvine, Irvine, California, USA
* jsanjak@uci.edu (JSS); krthornt@uci.edu (KRT)
Abstract
The genetic component of complex disease risk in humans remains largely unexplained. A
corollary is that the allelic spectrum of genetic variants contributing to complex disease risk
is unknown. Theoretical models that relate population genetic processes to the maintenance
of genetic variation for quantitative traits may suggest profitable avenues for future experi-
mental design. Here we use forward simulation to model a genomic region evolving under a
balance between recurrent deleterious mutation and Gaussian stabilizing selection. We
consider multiple genetic and demographic models, and several different methods for identi-
fying genomic regions harboring variants associated with complex disease risk. We demon-
strate that the model of gene action, relating genotype to phenotype, has a qualitative effect
on several relevant aspects of the population genetic architecture of a complex trait. In par-
ticular, the genetic model impacts genetic variance component partitioning across the allele
frequency spectrum and the power of statistical tests. Models with partial recessivity closely
match the minor allele frequency distribution of significant hits from empirical genome-wide
association studies without requiring homozygous effect sizes to be small. We highlight a
particular gene-based model of incomplete recessivity that is appealing from first principles.
Under that model, deleterious mutations in a genomic region partially fail to complement
one another. This model of gene-based recessivity predicts the empirically observed incon-
sistency between twin and SNP based estimated of dominance heritability. Furthermore,
this model predicts considerable levels of unexplained variance associated with intralocus
epistasis. Our results suggest a need for improved statistical tools for region based genetic
association and heritability estimation.
Author Summary
Gene action determines how mutations affect phenotype. When placed in an evolutionary
context, the details of the genotype-to-phenotype model can impact the maintenance of
genetic variation for complex traits. Likewise, non-equilibrium demographic history may
PLOS Genetics | DOI:10.1371/journal.pgen.1006573 January 19, 2017 1 / 30
a1111111111
a1111111111
a1111111111
a1111111111
a1111111111
OPEN ACCESS
Citation: Sanjak JS, Long AD, Thornton KR (2017)
A Model of Compound Heterozygous, Loss-of-
Function Alleles Is Broadly Consistent with
Observations from Complex-Disease GWAS
Datasets. PLoS Genet 13(1): e1006573.
doi:10.1371/journal.pgen.1006573
Editor: Simon Gravel, McGill University, CANADA
Received: April 18, 2016
Accepted: January 5, 2017
Published: January 19, 2017
Copyright: © 2017 Sanjak et al. This is an open
access article distributed under the terms of the
Creative Commons Attribution License, which
permits unrestricted use, distribution, and
reproduction in any medium, provided the original
author and source are credited.
Data Availability Statement: Our simulation code
and code for downstream analyses are freely
available at: http://github.com/ThorntonLab/
disease_sims, http://github.com/molpopgen/
buRden, http://github.com/molpopgen/fwdpy, and
http://github.com/molpopgen/TennessenEAonly.
Funding: This work was supported by NIH grant
R01-GM115564 to KRT. This work was supported
by NIH grant R01-GM115562 to ADL. This material
is based upon work supported by the National
Science Foundation Graduate Research Fellowship
Program under Grant No. DGE-1321846. Any

affect patterns of genetic variation. Here, we explore the impact of genetic model and pop-
ulation growth on distribution of genetic variance across the allele frequency spectrum
underlying risk for a complex disease. Using forward-in-time population genetic simula-
tions, we show that the genetic model has important impacts on the composition of
variation for complex disease risk in a population. We explicitly simulate genome-wide
association studies (GWAS) and perform heritability estimation on population samples. A
particular model of gene-based partial recessivity, based on allelic non-complementation,
aligns well with empirical results. This model is congruent with the dominance variance
estimates from both SNPs and twins, and the minor allele frequency distribution of
GWAS hits.
Introduction
Risk for complex diseases in humans, such as diabetes and hypertension, is highly heritable yet
the causal DNA sequence variants responsible for that risk remain largely unknown. Genome-
wide association studies (GWAS) have found many genetic markers associated with disease
risk [1]. However, follow-up studies have shown that these markers explain only a small por-
tion of the total heritability for most traits [2, 3].
There are many hypotheses which attempt to explain the ‘missing heritability’ problem [2–
5]. Genetic variance due to epistatic or gene-by-environment interactions is difficult to identify
statistically because of, among other reasons, increased multiple hypothesis testing burden [6,
7], and could artificially inflate estimates of broad-sense heritability [8]. Well-tagged interme-
diate frequency variants may not reach genome-wide significance in an association study if
they have smaller effect sizes [9, 10]. One appealing verbal hypothesis for this ‘missing herita-
bility’ is that there are rare causal alleles of large effect that are difficult to detect [4, 11, 12].
These hypotheses are not mutually exclusive, and it is probable that a combination of models
will be needed to explain all heritable disease risk [13].
The standard GWAS attempts to identify genetic polymorphisms that differ in frequency
between cases and controls. A complementary approach is to estimate the heritability
explained by genotyped (and imputed) markers (SNPs) under different population sampling
schemes [14, 15]. Stratifying markers by minor allele frequency (MAF) prior to performing
SNP-based heritability estimation allows the partitioning of genetic variation across the allele
frequency spectrum to be estimated [16], which is an important summary of the genetic archi-
tecture of a complex trait [16–23]. This approach has inferred a contribution of rare alleles to
genetic variance in both human height and body mass index (BMI) [16], consistent with theo-
retical work showing that rare alleles will have large effect sizes if fitness effects and trait effects
are correlated [18, 20–25]. Yet, simulations of causal loci harboring multiple rare variants with
large additive effects predict an excess of low-frequency significant markers relative to empiri-
cal findings [4, 26].
SNP-based heritability estimates have concluded that there is little missing heritability for
height and BMI, and that the causal loci simply have effect sizes that are too small to reach
genome-wide significance under current GWAS sample sizes [14, 16]. Further, extensions to
these methods decompose genetic variance into additive and dominance components and find
that dominance variance is approximately one fifth of the additive genetic variance on average
across seventy-nine complex traits [27]. When taken into account together with results from
GWAS, these observations can be interpreted as evidence that the genetic architecture of
human traits is best-explained by a model of small additive effects. However, a recent large
Compound Heterozygosity and Complex Traits
PLOS Genetics | DOI:10.1371/journal.pgen.1006573 January 19, 2017 2 / 30
opinions, findings, and conclusions or
recommendations expressed in this material are
those of the authors and do not necessarily reflect
the views of the National Science Foundation. The
funders had no role in study design, data collection
and analysis, decision to publish, or preparation of
the manuscript.
Competing Interests: The authors have declared
that no competing interests exist.

twin study found a substantial contribution of dominance variance for fourteen out of eighteen
traits [28]. The reason for this discrepancy in results remains unclear. One possibility is a sta-
tistical artifact; for example, twin studies may be prone to mistakenly infer non-additive effects
when none exist. Another possibility, which we return to later, is that this apparently contra-
dictory results are expected under a different model of gene action.
The design, analysis, and interpretation of GWAS are heavily influenced by the “standard
model” of quantitative genetics [29]. This model assigns an effect size to a mutant allele, but
formally makes no concrete statement regarding the molecular nature of the allele. Early appli-
cations of this model to the problem of human complex traits include Risch’s work on the
power to detect causal mutations [30, 31] and Pritchard’s work showing that rare alleles under
purifying selection may contribute to heritable variation in complex traits [17]. When applied
to molecular data, such as SNP genotypes in a GWAS, these models treat the SNPs themselves
as the loci of interest. For example, influential power studies informing the design of GWAS
assign effect sizes directly to SNPs and assume Risch’s model of multiplicative epistasis [32].
Similarly, the single-marker logistic regression used as the primary analysis of GWAS data
typically assumes an additive or recessive model at the level of individual SNPs [33]. Finally,
recent methods designed to estimate the heritability of a trait explained by genotyped markers
assigns additive and dominance effects directly to SNPs [14, 16, 27, 34]. Naturally, the results
of such analyses are interpreted in light of the assumed model of gene action.
A weakness of the multiplicative epistasis model [30, 31] when applied to SNPs is that the
concept of a gene, defined as a physical region where loss-of-function mutations have the same
phenotype [35], is lost. Specifically, under the standard model, the genetic concept of a failure
to complement is a property of SNPs and not “gene regions” (see [36] for a detailed discussion
of this issue). We have recently introduced an alternative model of gene action, one in which
risk mutations are unconditionally deleterious and fail to complement at the level of a “gene
region” [36]. This model, influenced by the standard operational definition of a gene [35],
gives rise to the sort of allelic heterogeneity typically observed for human Mendelian diseases
[37], and to a distribution of GWAS “hit” minor allele frequencies [4, 26] consistent with
empirical results [36]. In this article, we explore this “gene-based” model under more complex
demographic scenarios as well as its properties with respect to the estimation of variance com-
ponents using SNP-based approaches [34] and twin studies. We also compare this model to
the standard models of strictly additive co-dominant effects, and multiplicative epistasis with
dominance.
We further explore the power of several association tests to detect a causal gene region
under each genetic and demographic model. We find significant heterogeneity in the perfor-
mance of burden tests [36, 38, 39] across models of the trait and demographic history. We find
that population expansion reduces the power to detect causal gene-regions due to an increase
in rare variation, in agreement with work by [22, 23]. The behavior of the tests under different
models provides us with insight as to the circumstances in which each test is best suited.
In total, our results show that modeling gene action is key to modeling GWAS, and thus
plays an important role in both the design and interpretation of such studies. Further, the
model of gene-based recessivity best explains the differences between estimates of additive and
dominance variance components from SNP-based methods [27] and from twin studies [28]
and is consistent with the distribution of frequencies of significant associations in GWAS [4,
26]. Further, the genetic model plays a much more important role than the demographic
model, which is expected based on previous work on additive models showing that the genetic
load is approximately unaffected by changes in population size over time, [21, 22]. Consistent
with recent work by [23], we find that rapid population growth in the recent past increases the
contribution of rare variants to total genetic variance. However, we show here that different
Compound Heterozygosity and Complex Traits
PLOS Genetics | DOI:10.1371/journal.pgen.1006573 January 19, 2017 3 / 30

models of gene action are qualitatively different with respect to the partitioning of genetic vari-
ance across the allele frequency spectrum. We also show that these conclusions hold under the
more complex demographic models that have been proposed for human populations [21, 40].
Results and Discussion
The models
As in [36],we simulate a 100 kilobase region of human genome, contributing to a complex dis-
ease phenotype and fitness. The region evolves forward in time subject to neutral and deleteri-
ous mutation, recombination, selection, and drift. To perform genetic association and
heritability estimation studies in silico, we need to impose a trait onto simulated individuals. In
doing so, we introduce strong assumptions about the molecular underpinnings of a trait and
its evolutionary context.
How does the molecular genetic basis of a trait under natural selection influence population
genetic signatures in the genome? This question is very broad, and therefore it was necessary
to restrict ourselves to a small subset of molecular and evolutionary scenarios. We analyzed a
set of approaches to modeling a single gene region experiencing recurrent unconditionally-
deleterious mutation contributing to a quantitative trait subject to Gaussian stabilizing selec-
tion. The expected fitness effect of a mutation is always deleterious because trait effects are
sampled from an exponential distribution. Therefore, we do not allow for compensatory muta-
tions that may occur in more general models of stabilizing selection. Specifically, we studied
three different genetic models and two different demographic models, holding the fitness
model as a constant. Parameters are briefly described in Table 1.
We implemented three disease-trait models of the phenotypic form P = G + E. G is the
genetic component, and E ¼ Nð0; s
2
e
Þ is the environmental noise expressed as a Gaussian ran-
dom variable with mean 0 and variance s
2
e
. In this context, s
2
e
should be thought of as both the
contribution from the environment and from the remaining genetic variance at loci in linkage
equilibrium with the simulated 100kb region. The genetic models are named the additive co-
dominant (AC) model, multiplicative recessive (Mult. recessive; MR) model and the gene-
based recessive (GBR) model. The MR model has a parameter, h, that controls the degree of
Table 1. Description of parameters used in the models.
Parameter Description
N Population size
P Phenotype
P
opt
Optimum phenotype
G Genetic contribution to phenotype
E Environmental contribution to phenotype
λ Mean and standard deviation of trait effects
c
i
Specific trait effect of site i
h Dominance coefficient for trait effects
w Fitness, based on Gaussian function
s
2
s
The total inverse selection intensity
s
2
e
Environmental variance
V
A
Additive genetic variance
V
D
Dominance genetic variance
V
G
Genetic variance
V
A;q  x
Additive variance explained by variance below frequency q
doi:10.1371/journal.pgen.1006573.t001
Compound Heterozygosity and Complex Traits
PLOS Genetics | DOI:10.1371/journal.pgen.1006573 January 19, 2017 4 / 30

recessivity; we call this model the complete MR (cMR) when h = 0 and the incomplete MR
(iMR) when 0 h  1. Here, h = 1 corresponds to co-dominance, which is different from the
typical formulation used when modeling the fitness effects of mutations directly. It is also
important to note that here recessivity is being defined in terms of phenotypic effects; this may
be unusual for those more accustomed to dealing directly with recessivity for fitness effects.
An idealized relationship between dominance for fitness effects and trait effects of a mutation
on an unaffected genetic background is shown in S15 Fig.
The critical conceptual difference between recessive models is whether dominance is a
property of a locus (nucleotide/SNP) in a gene or the gene overall. Mathematically, this
amounts to whether one first determines diploid genotypes at sites (and then multiplies across
sites to get a total genetic effect) or calculates a score for each haplotype (the maternal and
paternal alleles). For completely co-dominant models, this distinction is irrelevant, however
for a model with arbitrary dominance one needs to be more specific. As an example, imagine a
compound heterozygote for two biallelic loci, i.e. genotype Ab/aB. In the case of traditional
multiplicative recessivity the compound heterozygote is wild type for both loci and therefore
wild-type over all; this implies that these loci are in different genes (or independent functional
units of the same gene) because the mutations are complementary. However, in the case of
gene-based recessivity [36], neither haplotype is wild-type and so the individual is not wild-
type; the failure of mutant alleles to complement defines these loci as being in the same gene
[35].
For a diploid with m
i
causative mutations on the i
th
haplotype, we may define the additive
model as
G
AC
¼
X
2
i¼1
X
m
i
j¼1
c
i;j
; ð1Þ
where c
i,j
is the effect size of the j
th
mutation on the i
th
haplotype. Each c
i,j
is sampled from an
exponential distribution with mean of λ, to reflect unconditionally deleterious mutation. In
other words, when a new mutation arises its effect c is drawn from an exponential distribution,
and remains constant throughout its entire sojourn in the population.
The GBR model is the geometric mean of the sum of effect sizes on each haplotype [36].
We sum the causal mutation effects on each allele (paternal and maternal) to obtain a haplo-
type score. We then take the square root of the product of the haplotype scores to determine
the total genetic value of the diploid.
G
GBR
¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
X
m
1
j¼1
c
1;j

X
m
2
j¼1
c
2;j
v
u
u
t
ð2Þ
Finally, the MR model depends on the number of positions for which a diploid is heterozy-
gous (m
Aa
) or homozygous (m
aa
) for causative mutations,
G
MR
¼
Y
m
Aa
j¼1
ð1 þhc
j
Þ
!
Y
m
aa
j¼1
ð1 þ 2c
j
Þ
!
 1: ð3Þ
Thus, h = 0 is a model of multiplicative epistasis with complete recessivity (cMR), and h = 1
closely approximates the additive model when effect sizes are small.
Here, phenotypes are subject to Gaussian stabilizing selection with an optimum at zero and
standard deviation of σ
s
= 1 such that the fitness, w, of a diploid is proportional to a Gaussian
Compound Heterozygosity and Complex Traits
PLOS Genetics | DOI:10.1371/journal.pgen.1006573 January 19, 2017 5 / 30