Parameterizing the genetic architecture under stabilizing selection

This paper derives a mechanistic, evolutionary-based linear mixed model that replaces the phenomenological $\alpha$-model by linking effect size frequency dependence to interpretable parameters like mutational variance and selection intensity, thereby enabling both the inference of evolutionary dynamics and improved genetic prediction via REML and BLUP.

The Gist

Imagine you are trying to understand why some people are very tall while others are short. You know that thousands of tiny genetic "switches" (variants) contribute to height. But here's the puzzle: rare switches (ones that only a few people have) tend to have a huge impact on height, while common switches (ones almost everyone has) tend to have a tiny impact.

For years, scientists have used a mathematical shortcut to describe this pattern. They called it the "$\alpha$-model." Think of this like a mechanic using a generic "rule of thumb" to fix a car engine. They know that turning the bolt a certain way makes the car run, but they don't really know why the engine was built that way in the first place. It works for prediction, but it's a bit of a black box.

This paper, written by Hanbin Lee and Jonathan Terhorst, says: "Let's stop guessing the rule of thumb. Let's look at the engine's blueprints."

Here is the story of their new approach, explained simply:

1. The Old Way: The "Rule of Thumb" ($\alpha$-model)

Imagine you are trying to predict the weather. The old method says, "If the barometer drops, it will rain." It's a statistical observation. It works, but it doesn't explain why the pressure dropped.
In genetics, the old model says: "If a genetic variant is rare, its effect size is big." It fits the data well, but it's just a curve drawn on a graph. It doesn't tell us about the evolutionary forces (like natural selection) that shaped these genes over thousands of years. Also, if a variant is extremely rare, this old math breaks down and gives impossible answers (like infinity).

2. The New Way: The "Evolutionary Blueprint"

The authors decided to build a model based on Fisher's Geometric Model.

• The Analogy: Imagine a hiker trying to reach the top of a mountain (the "perfect" health or trait). The mountain has a peak (optimal fitness).
• The Problem: If the hiker takes a giant step (a large genetic mutation), they are likely to fall off a cliff or get stuck in a ravine. They are far from the peak.
• The Solution: If the hiker takes tiny, careful steps (small mutations), they are more likely to stay on the path and reach the top.

Natural Selection acts like a strict coach. It constantly pushes the population back toward the peak.

• Large mutations are dangerous. The coach (selection) kicks them out of the population quickly. So, they stay rare.
• Small mutations are safe. The coach lets them stick around. So, they become common.

The authors derived a new formula that mathematically describes this "coach." Instead of just guessing the relationship between rarity and impact, they calculated it based on how hard the coach is pushing (selection intensity) and how much the hiker tends to stumble (mutation variance).

3. Why This Matters: The "Two-Legged Stool"

The new model connects two worlds that usually don't talk to each other:

1. Evolutionary Biology: The story of how nature shapes us over time.
2. Statistical Genetics: The math we use to predict traits (like height or disease risk) from DNA data.

The "Aha!" Moment:
The authors showed that you can use standard statistical tools (called REML and BLUP) to estimate these evolutionary "coach" parameters directly from real-world data.

• Before: We could predict who might get a disease, but we didn't know the evolutionary story behind the numbers.
• Now: We can predict the disease and tell you, "This trait is under strong evolutionary pressure," or "This trait is mostly neutral."

4. The Results: A Better Map

The authors tested their new model using computer simulations (like a video game of evolution).

• The Test: They simulated a population evolving under natural selection and then tried to predict traits using both the old "rule of thumb" and their new "blueprint."
• The Winner: The new model was more accurate. It didn't just predict better; it correctly identified the underlying "variance" (the true amount of genetic influence).
• The Old Model's Flaw: The old model often got confused. It would think a trait had less genetic influence than it actually did, just because it couldn't account for the "coach" pushing the rare variants.

Summary in a Nutshell

Think of the old model as a weather vane. It spins and points in the right direction, but it doesn't tell you about the wind currents.
The new model is the meteorologist's computer. It understands the physics of the wind (evolution), allowing it to not only predict the storm (the trait) but also explain why the storm is happening and how strong the winds are.

By grounding their math in the actual mechanics of evolution, Lee and Terhorst have given geneticists a more powerful, more honest, and more accurate tool to understand the complex tapestry of human traits.

Technical Summary

1. Problem Statement

In statistical genetics, a well-documented empirical observation is that genetic variants with larger effect sizes on complex traits tend to occur at lower allele frequencies. This pattern is widely interpreted as a signature of stabilizing selection (purifying selection), where large-effect mutations are deleterious and removed from the population.

To model this relationship, the field currently relies on the $\alpha$-model. This phenomenological model assumes that the variance of effect sizes is inversely proportional to allelic heterozygosity raised to a power $\alpha$ ($0 \le \alpha \le 1$):
$$\text{Var}(\beta_j) \propto [p_j(1-p_j)]^{-\alpha}$$
While useful for curve-fitting, the $\alpha$-model has significant theoretical limitations:

• Lack of Mechanism: It is a statistical heuristic without a direct derivation from population genetic theory.
• Singularity: As the allele frequency $p_j \to 0$, the variance term diverges to infinity, requiring ad-hoc fixes (e.g., binning rare variants).
• Interpretability: The parameters ($\alpha$ and heritability scaling) are tuning parameters without explicit biological meaning regarding mutation rates or selection intensity.

The authors aim to bridge this gap by deriving a mechanistic, evolution-based alternative that links fitness landscapes directly to linear mixed models (LMMs) used in Genome-Wide Association Studies (GWAS).

2. Methodology

The authors develop a framework that integrates Fisher's Geometric Model (FGM) of stabilizing selection with Linear Mixed Models (LMMs).

A. Theoretical Derivation

1. Fitness Landscape: They assume individual fitness follows FGM, where fitness decays quadratically with the distance from an optimal trait vector.
2. Diffusion Approximation: They derive the stationary distribution of allele frequencies ($x_j$) under stabilizing selection and mutation using a diffusion approximation. This yields a probability density function dependent on the squared effect size on fitness ($\|\alpha_j\|^2$) and the selection width ($W_S$).
3. Coupling Traits: They introduce a pleiotropic framework where a "focal trait" (the one being studied) is correlated with a "latent selected trait" (the one under direct stabilizing selection).
   • Let $\alpha_j$ be the effect size vector on the latent fitness trait.
   • Let $\beta_j$ be the effect size on the focal trait.
   • They assume a joint Gaussian distribution for $(\alpha_j, \beta_j)$ characterized by variances $\sigma_a^2, \sigma_b^2$ and correlation $\rho_{ab}$.
4. Deriving the Variance Component:
   • Using the law of total expectation, they derive the conditional variance of the focal effect size given the observed genotype data ($G$).
   • By integrating over the distribution of fitness effects and allele frequencies, they derive a closed-form expression for the variance of the random effects in an LMM:
     $$\text{Var}(\beta_j | G) = \sigma_b^2 \left( 1 - \rho_{ab}^2 \frac{2 \frac{\sigma_a^2}{W_S} p_j(1-p_j)}{1 + 2 \frac{\sigma_a^2}{W_S} p_j(1-p_j)} \right)$$
   • This formula naturally captures the frequency dependence: as selection ($\sigma_a^2/W_S$) increases, the variance of effect sizes for common alleles decreases relative to rare alleles.

B. Statistical Implementation

• Model Formulation: The derived variance structure is embedded into a standard Linear Mixed Model: $y = G\beta + \epsilon$.
• Estimation: The evolutionary parameters ($\sigma_a^2/W_S$, $\sigma_b^2$, and $\rho_{ab}$) are estimated using Restricted Maximum Likelihood (REML).
• Prediction: Once parameters are estimated, the model generates Best Linear Unbiased Predictors (BLUP) for genetic values.

C. Validation

• Forward Simulations: The authors used SLiM to simulate a panmictic population under stabilizing selection with pleiotropy. Unlike standard benchmarks that fix genotypes, these simulations evolve genotypes and phenotypes jointly under the specified evolutionary forces.
• Comparison: They compared their evolutionary model against the standard $\alpha$-model baselines ($\alpha = 0, 0.5, 1$).

3. Key Contributions

1. Mechanistic Derivation: The paper provides the first population-genetic derivation of the frequency-dependent effect size variance, replacing the heuristic $\alpha$-model with a formula grounded in evolutionary theory (mutation, selection, and pleiotropy).
2. Resolution of Pathologies: The new model avoids the singularity at $p_j \to 0$ found in the $\alpha$-model, as the variance term remains bounded.
3. Identifiable Parameters: The model introduces interpretable parameters:
   • $\sigma_b^2$: Overall scale of genetic variance.
   • $\sigma_a^2/W_S$: A composite parameter representing the ratio of mutational variance to selection intensity (strength of stabilizing selection).
   • $\rho_{ab}$: The degree of coupling between the focal trait and the selected trait.
4. Unified Framework: It unifies fitness-landscape modeling with standard mixed-model methodology, allowing simultaneous inference of evolutionary parameters and genetic prediction.

4. Results

A. Variance Component Recovery

• $\sigma_b^2$ (Overall Scale): REML accurately recovered the true $\sigma_b^2$ across all simulation scenarios, regardless of selection strength or coupling.
• $\sigma_a^2/W_S$ (Selection Sensitivity): The model successfully recovered the composite selection parameter. However, estimates were systematically lower than the theoretical single-locus expectation. The authors attribute this to the Bulmer effect (linkage disequilibrium reducing effective selection pressure), which they quantified and corrected for in the analysis.
• Coupling ($\rho_{ab}$): When $\rho_{ab} = 0$ (no coupling), the frequency dependence vanished, and the model correctly reduced to a frequency-independent prior.

B. Genetic Prediction (BLUP)

• Performance: The evolutionary model generally outperformed or matched the $\alpha$-model baselines in predicting held-out phenotypes (measured by $R^2$).
• Failure of $\alpha=1$: The standard $\alpha=1$ model (GCTA) performed poorly, especially under strong selection, failing to capture the architecture.
• The $\alpha=0$ Baseline: The evolutionary model performed similarly to the $\alpha=0$ model (constant variance) in many regimes. However, the authors note this is deceptive: while predictive accuracy was similar, the $\alpha=0$ model mischaracterized the underlying architecture by underestimating the total genetic variance ($\sigma_b^2$) to compensate for missing frequency dependence. The evolutionary model correctly identified the architecture without sacrificing predictive power.

5. Significance and Implications

• Theoretical Advancement: This work moves statistical genetics from "curve fitting" to "mechanistic modeling." It demonstrates that the frequency dependence of effect sizes is a natural consequence of stabilizing selection acting on pleiotropic traits.
• Practical Utility: The framework allows researchers to estimate evolutionary parameters (like selection strength) directly from GWAS summary statistics or genotype-phenotype data using standard REML software.
• Improved Inference: By correctly modeling the variance structure, the method provides unbiased estimates of genetic variance components, which are often biased in heuristic models.
• Future Directions: The authors acknowledge that the current derivation assumes linkage equilibrium (mean-field approximation). The systematic underestimation of selection strength suggests that future work should explicitly model Linkage Disequilibrium (LD) and the Bulmer effect to refine parameter estimates.

In summary, Lee and Terhorst provide a rigorous, evolutionarily grounded alternative to the $\alpha$-model, enabling more accurate inference of genetic architecture and evolutionary history from complex trait data.