Stochastic Eco-Evolutionary Dynamics of Multivariate Traits
A stochastic neutral baseline for G-matrix orientation
Bob Week
KiTE Postdoc · Schulenburg Group · Kiel University
bobweek.github.io · @bweek.bsky.social
The \(\mathbf G\)-matrix: heritable variation with a direction
\(\mathbf G\) = additive genetic (co)variances of \(d\) traits.
- diagonal — heritable variance of each trait
- off-diagonal — heritable covariance between traits
We care about its direction: the breeder’s equation
\[
\Delta\bar{\mathbf z}=\mathbf G\,\boldsymbol\beta
\]
bends the response to selection along \(\mathbf G\), not along the selection gradient \(\boldsymbol\beta\).
Drift: the familiar allele-frequency picture
Random reproductive success
[ ]
random genetic drift
This intuition is pathwise: replicate populations do different things.
Drift scatters the orientation of \(\mathbf G\)
One ancestor, many replicate lines. Drift erodes variation on average — yet each line’s \(\mathbf G\) can end up pointing a different way.
Question: what neutral model predicts this pathwise divergence in orientation?
Roadmap
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.
Classical drift prediction for (G)
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.
The missing null
Drift is known to:
- erode heritable variation;
- generate random associations among loci;
- make replicate populations diverge.
But for (G)-matrices, the classical result mostly tells us about (E[G_t]).
What about variation around that expectation?
Strategy: derive (G) dynamics from a population process
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}.
\]
Heuristic picture of the diffusion limit
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.
Moment equations from one population process
● growth / selection ● mutation ● drift (noise)
Abundance, mean trait and covariance fall out of the same stochastic process:
\[\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{\tfrac vn\mathbf P}\,\mathrm d\mathbf B_{\bar{\mathbf z}}}\]
\[
\begin{aligned}
\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{-\tfrac vn\mathbf P}\,\Big]\mathrm dt\\
&+\violet{\sqrt{\tfrac vn(\mathbf K-\mathbf P\otimes\mathbf P)}:\mathrm d\mathbf B_{\mathbf P}}
\end{aligned}
\]
Not closed: \(\mathbf P\) needs the kurtosis \(\mathbf K\), and so on up the hierarchy.
Gaussian closure: a usable \(\mathbf G\) equation
Assume additive genetic values are multivariate normal, with \(\mathbf z=\mathbf g+\mathbf e\) and \(\mathrm{Cov}(\mathbf z)=\mathbf G+\mathbf E\). Then the kurtosis closes, \(\mathbf K\to\mathbf G\obar\mathbf G+\mathbf G\ubar\mathbf G\), and \(\mathbf G\) follows
\[
\begin{aligned}
\mathrm d\mathbf G=\Big[\,&\lime{\mathbf M}+\pink{2\mathbf G(\nabla_{\mathbf G}\bar m-\overline{\nabla_{\mathbf G}m})\mathbf G}\;\violet{-\tfrac vn\mathbf G}\,\Big]\mathrm dt\\
&+\violet{\sqrt{\tfrac vn(\mathbf G\obar\mathbf G+\mathbf G\ubar\mathbf G)}:\mathrm d\mathbf B_{\mathbf G}}
\end{aligned}
\]
Everything now closes in \(n,\bar{\mathbf z},\mathbf G\) — the objects quantitative geneticists already use.
Neutral (G)-matrix dynamics under drift
\[
\mathrm d\mathbf G=
\violet{-\frac{1}{N_e}\mathbf G\,\mathrm dt}
+
\violet{\sqrt{\frac{1}{N_e}
(\mathbf G\ubar\mathbf G+\mathbf G\obar\mathbf G)}:\mathrm d\mathbf B}
\]
where
\[
\frac{1}{N_e}=\frac vn,
\qquad
N_e=\frac nv.
\]
Assumptions for this neutral baseline:
- (m(,z));
- (M=);
- constant (n) and (v);
- no explicit recombination, selection or new mutation.
From (G) to genetic correlations
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).
\]
Genetic correlations tend towards ()
\[
\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.
Why the drift term does not contradict this
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.
\]
Same expectation, different pathwise reality
Classical expectation:
\[
\mathbb E[\mathbf G_t]=\mathbf G_0e^{-t/N_e}
\]
Realized neutral paths:
- erode heritable variation;
- alter (G) orientation;
- push genetic correlations toward extremes.
The expectation-level null is not wrong. It is incomplete as a pathwise null for orientation.
Recombination opposes the drift effect
The baseline is drift-only. Recombination breaks up allelic associations, pushing genetic covariances toward zero — it counteracts the pull to \(\pm1\).
Yet replicate \(\mathbf G\) orientations diverged even in outcrossing C. elegans (Mallard et al. 2023):
Open problem: the drift–recombination balance of genetic correlations.
What the framework buys
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.
Takeaways
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.
More in the paper
![]()
Numerical implementations: github.com/bobweek/multi-mtgl
Thanks & Questions
Steve Krone, Peter L. Ralph, Hinrich Schulenburg, Patrick C. Phillips, Arne Traulsen, Jonas Wickman, Brendan Bohannan
Backup: how to choose (m)
The framework is useful when a biological problem can be described through growth rates or fitness functions.
Ecological examples
- Logistic growth + stabilizing selection
\[
m(\nu,\mathbf z)=r-\tfrac12\mathbf z^\top\mathbf\Psi\mathbf z-cn
\]
- Lotka–Volterra interactions
\[
m_i=r_i-\sum_j\alpha_{ij}n_j
\]
Evolutionary examples
- coevolution;
- evolutionary rescue;
- host–parasite matching;
- frequency-dependent interactions;
- density-dependent eco-evolutionary feedbacks.
Backup: analytical and numerical use
Numerical use:
- use the Brownian-motion form;
- integrate with Euler–Maruyama or SDE solvers;
- see
github.com/bobweek/multi-mtgl.
Analytical use:
- use Itô’s formula for (df(n,{z},G));
- use martingale rules to combine correlated stochastic terms;
- avoid inventing independent noise processes where the population process says they covary.
Backup: full Gaussian closure
\[
\begin{aligned}
\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\obar\mathbf G+\mathbf G\ubar\mathbf G)}:\mathrm d\mathbf B_{\mathbf G}.
\end{aligned}
\]
Backup: reproductive variance
(v=) reproductive variance
Higher (v) means stronger demographic stochasticity and stronger random genetic drift at fixed (n).
Backup: fixation intuition
(v>0) random genetic drift
