A stochastic neutral baseline for G-matrix orientation

One ancestral population.
Many replicate lines.
Drift erodes genetic variation on average.
But replicate populations can diverge in the orientation of genetic variation.
Question: what neutral stochastic model should explain this pathwise divergence?
1. Biology
What is (G), and what does classical drift theory predict?
2. Framework
Derive stochastic moment dynamics from one population process.
3. Result
A neutral model predicts extreme pathwise genetic correlations.
The goal is not to replace multilocus or developmental theory. It is to supply a stochastic neutral baseline for (G)-matrix orientation.
For multivariate traits,
\[ \mathbf z=(z_1,\ldots,z_d)^\top . \]
The additive genetic covariance matrix is
\[ \mathbf G= \begin{pmatrix} G_{11} & G_{12} & \cdots \\ G_{21} & G_{22} & \cdots \\ \vdots & \vdots & \ddots \end{pmatrix}. \]
Diagonals: additive genetic variances.
Off-diagonals: additive genetic covariances.
Genetic correlation:
\[ \rho_{ij}(\mathbf G)= \frac{G_{ij}}{\sqrt{G_{ii}G_{jj}}}. \]

Random reproductive success
[ ]
random genetic drift

This intuition is pathwise: replicate populations do different things.
For quantitative traits, the standard neutral expectation is
\[ \mathbb E[\mathbf G_t]=\mathbf G_0 e^{-t/N_e}. \]
So drift erodes heritable variation.
But every entry shrinks proportionally:
\[ \rho\big(\mathbb E[\mathbf G_t]\big)=\rho(\mathbf G_0). \]


The expectation is clear. The neutral distribution of realized (G)-matrix orientations is not.

Drift is known to:
But for (G)-matrices, the classical result mostly tells us about (E[G_t]).
What about variation around that expectation?
A diffusion approximation separates robust population-level ingredients from microscopic turnover details.
Individual turnover
births, deaths, mutation, reproductive variance
Measure-valued diffusion
trait distribution evolves stochastically in (R^d)
Moment dynamics
(n), ({z}), (G), and derived quantities
Retained macroscopic ingredients:
\[ \pink{m(\nu,\mathbf z)}\quad\text{growth / selection}, \qquad \lime{\mathbf M}\quad\text{mutation covariance}, \qquad \violet{v}\quad\text{reproductive variance}. \]
For one useful scaling,
\[ k\left[W^{1/k}(\mathcal P_t^{(k)},\mathbf z)-1\right] \longrightarrow m(\nu_t,\mathbf z). \]
Different individual-based models can share the same limit.
The point is not to keep every microscopic detail; it is to identify what survives at the population scale.
The same stochastic population process gives consistent dynamics for abundance, mean traits and covariance.
\[ \mathrm dn=\pink{\bar m n}\,\mathrm dt+\violet{\sqrt{vn}\,\mathrm dB_n} \]
\[ \mathrm d\bar{\mathbf z}=\pink{\mathrm{Cov}(m,\mathbf z)}\,\mathrm dt +\violet{\sqrt{\frac vn\mathbf P}\,\mathrm d\mathbf B_{\bar{\mathbf z}}} \]
\[\begin{multline} \mathrm d\mathbf P=\Big[\lime{\mathbf M}+\pink{\mathrm{Cov}\big(m,(\mathbf z-\bar{\mathbf z})(\mathbf z-\bar{\mathbf z})^\top\big)}\violet{-\frac vn\mathbf P}\Big]\,\mathrm dt \\ +\violet{\sqrt{\frac vn(\mathbf K-\mathbf P\otimes\mathbf P)}:\mathrm d\mathbf B_{\mathbf P}}. \end{multline}\]
This hierarchy is not closed in general.
Assume additive genetic values are multivariate normal and traits decompose as
\[ \mathbf z=\mathbf g+\mathbf e, \qquad \mathrm{Cov}(\mathbf z)=\mathbf G+\mathbf E. \]
Then the (G) dynamics close as
\[\begin{multline} \mathrm d\mathbf G=\bigg[\lime{\mathbf M} +\pink{2\mathbf G(\nabla_{\mathbf G}\bar m-\overline{\nabla_{\mathbf G}m})\mathbf G} \violet{-\frac vn\mathbf G}\bigg]\,\mathrm dt \\ +\violet{\sqrt{\frac vn(\mathbf G\overline\otimes\mathbf G+ \mathbf G\underline\otimes\mathbf G)}:\mathrm d\mathbf B_{\mathbf G}}. \end{multline}\]
This is the part needed for a stochastic neutral model of (G).
\[ \mathrm d\mathbf G= \violet{-\frac{1}{N_e}\mathbf G\,\mathrm dt} + \violet{\sqrt{\frac{1}{N_e} (\mathbf G\underline\otimes\mathbf G+\mathbf G\overline\otimes\mathbf G)}:\mathrm d\mathbf B} \]
where
\[ \frac{1}{N_e}=\frac vn, \qquad N_e=\frac nv. \]
Assumptions for this neutral baseline:
Apply Itô’s formula to
\[ \rho=f(\mathbf G)=\frac{G_{ij}}{\sqrt{G_{ii}G_{jj}}}. \]
Component form:
\[ \mathrm dG_{ij}= -\frac vnG_{ij}\,\mathrm dt +\sqrt{\frac vn\left(G_{ii}G_{jj}+G_{ij}^2\right)}\,\mathrm dB_{G_{ij}}. \]
The Brownian terms covary. The martingale version keeps those covariances correct while combining stochastic terms:
\[ a\,\mathrm d\mathcal M(x)+b\,\mathrm d\mathcal M(y) =\mathrm d\mathcal M(ax+by). \]
\[ \mathrm d\rho= \pink{-\frac{1}{2N_e}\rho(1-\rho^2)}\,\mathrm dt + \violet{\sqrt{\frac{1}{N_e}(1-\rho^2)}\,\mathrm dB}. \]


Under drift alone in this isolated neutral baseline, realized populations concentrate near extreme genetic correlations.
The Itô drift term points back toward zero:
\[ -\frac{1}{2N_e}\rho(1-\rho^2). \]
But the diffusion coefficient vanishes at the boundaries:
\[ \sqrt{\frac{1}{N_e}(1-\rho^2)}. \]
So the process spends increasing time near (=).
Equivalently, with (u=()),
\[ \mathrm du = \sqrt{\frac{1}{N_e}}\,\frac{1}{\sqrt{1-\rho^2}}\,\mathrm dB. \]
Classical expectation:
\[ \mathbb E[\mathbf G_t]=\mathbf G_0e^{-t/N_e} \]
Realized neutral paths:

The expectation-level null is not wrong. It is incomplete as a pathwise null for orientation.


Explicit recombining extensions are still needed. But these data reinforce the need for a stochastic neutral orientation baseline.
Not just “moments are useful.”
The useful point is that (n), ({z}) and (G) are derived from one stochastic population process, so their stochastic terms and cross-covariances are consistent.
That makes it possible to derive new stochastic dynamics for quantities such as
\[ \rho(\mathbf G),\qquad \bar m(n,\bar{\mathbf z},\mathbf G),\qquad \alpha(\bar{\mathbf z}_1,\bar{\mathbf z}_2), \]
rather than specifying noise terms ad hoc.
Replicate populations can diverge in (G)-matrix orientation under drift.
The classical expectation for (E[G_t]) remains correct, but it does not describe typical realized orientations.
The drift-only stochastic neutral baseline predicts genetic correlations tending toward ().
A stochastic moment framework gives a reusable way to derive pathwise population-level models from underlying turnover.
Numerical implementations: github.com/bobweek/multi-mtgl
Steve Krone, Peter L. Ralph, Hinrich Schulenburg, Patrick C. Phillips, Arne Traulsen, Jonas Wickman, Brendan Bohannan




The framework is useful when a biological problem can be described through growth rates or fitness functions.
Ecological examples
\[ m(\nu,\mathbf z)=r-\tfrac12\mathbf z^\top\mathbf\Psi\mathbf z-cn \]
\[ m_i=r_i-\sum_j\alpha_{ij}n_j \]
Evolutionary examples
Numerical use:
github.com/bobweek/multi-mtgl.Analytical use:
\[\begin{align} \mathrm dn &= \bar m n\,\mathrm dt+\sqrt{vn}\,\mathrm dB_n,\\ \mathrm d\bar{\mathbf z} &= \mathbf G(\nabla_{\bar{\mathbf z}}\bar m-\overline{\nabla_{\bar{\mathbf z}}m})\,\mathrm dt +\sqrt{\frac vn\mathbf G}\,\mathrm d\mathbf B_{\bar{\mathbf z}},\\ \mathrm d\mathbf G &= \left[\mathbf M+2\mathbf G(\nabla_{\mathbf G}\bar m-\overline{\nabla_{\mathbf G}m})\mathbf G-\frac vn\mathbf G\right]\mathrm dt\\ &\qquad+\sqrt{\frac vn(\mathbf G\overline\otimes\mathbf G+\mathbf G\underline\otimes\mathbf G)}:\mathrm d\mathbf B_{\mathbf G}. \end{align}\]
The density animation is only a visualization.
The formal object is a measure-valued population process characterized by a martingale problem.
For a test function (x),
\[ \mathrm dN_t(x)=N_t(Lx)\,\mathrm dt+\mathrm d\mathcal M_t(x), \]
with quadratic covariation
\[ \mathrm d\mathcal M_t(x)\,\mathrm d\mathcal M_t(y) = v\,N_t(xy)\,\mathrm dt. \]
This is the source of the stochastic covariance structure in the moment equations.

(v=) reproductive variance
Higher (v) means stronger demographic stochasticity and stronger random genetic drift at fixed (n).
(v>0) random genetic drift
