# Rheology of bidisperse non-Brownian suspensions

Abhinendra Singh<EMAIL_ADDRESS>Department of Macromolecular
Science and Engineering, Case Western Reserve University, Cleveland, OH,
10040, USA  James Franck Institute, University of Chicago, Chicago, Illinois
60637, USA  Pritzker School of Molecular Engineering, University of Chicago,
Chicago, Illinois 60637, USA  Christopher Ness  School of Engineering,
University of Edinburgh, Edinburgh EH9 3FG, United Kingdom  Abhishek K. Sharma
Pritzker School of Molecular Engineering, University of Chicago, Chicago,
Illinois 60637, USA  Juan J de Pablo  Pritzker School of Molecular
Engineering, University of Chicago, Chicago, Illinois 60637, USA  Materials
Science Division, Argonne National Laboratory, Lemont, Illinois 60439, USA
Heinrich M Jaeger  James Franck Institute, University of Chicago, Chicago,
Illinois 60637, USA  Department of Physics, The University of Chicago,
Chicago, Illinois 60637, USA

###### Abstract

We study the rheology of bidisperse non-Brownian suspensions using particle-
based simulation, mapping the viscosity as a function of the size ratio of the
species, their relative abundance, and the overall solid content. The
variation of the viscosity with applied stress exhibits paradigmatic shear
thickening phenomenology irrespective of composition, though the stress-
dependent limiting solids fraction governing the viscosity and its divergence
point are non-monotonic in the mixing ratio. Contact force data demonstrate an
exchange in dominant stress contribution from large-large to small-small
particle contacts as the mixing ratio of the species evolves. Combining a
prior model for shear thickening with one for composition-dependent jamming,
we obtain a full description of the rheology of bidisperse non-Brownian
suspensions capable of predicting effects such as the viscosity reduction
observed upon adding small particle fines to a suspension of large particles.

## Introduction.

Suspensions of small particles, radius $a\approx$
$100\text{\,}\mathrm{nm}$$-$$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$, form a
class of complex fluids abundant in nature and industry Ness  _et al._ (2022);
Morris (2020); Singh (2023). Their widespread use calls for detailed
constitutive characterization to enable reliable process design Roussel  _et
al._ (2010), especially in the dense regime where particles and fluid are
mixed roughly equally Brown and Jaeger (2014). Under external deformation,
these systems, which are apparently simple in composition, exhibit complex
rheology including yielding, shear thinning, shear thickening, and jamming
Brown and Jaeger (2014); Peters  _et al._ (2016); Singh  _et al._ (2019).
Recently this phenomenology has been linked to microscopic physics,
specifically constraints that control the relative translation and rotation of
interacting particle pairs Guy  _et al._ (2018); Singh  _et al._ (2020, 2022).
Shear thickening, for example, represents a crossover from unconstrained to
constrained tangential motion as the imposed particle stress $\sigma$ exceeds
a threshold set by the repulsive force $\sigma_{0}\sim F_{0}/a^{2}$ Seto  _et
al._ (2013). A mean-field approach by Wyart and Cates (2014) (WC) captures the
transition using stress-dependent jamming volume fraction $\phi_{J}(\sigma)$
interpolating between low ($\phi_{0}$) and high ($\phi_{m}$) stress limits,
reproducing (in some cases quantitatively) the steady state rheology Guy  _et
al._ (2018); Singh  _et al._ (2018).

The above conceptual framework was devised based on nearly monodisperse
suspensions, and most numerical and experimental works that seek to test it
reflect this Singh  _et al._ (2018); Guy  _et al._ (2015); Ness and Sun
(2016); Guy  _et al._ (2020); Pradeep  _et al._ (2021). As soon as significant
deviations from monodispersity are considered, however, complexity emerges
that is not captured by WC. In particular, adding small particles to a system
of large ones can reduce the viscosity under shear Roussel  _et al._ (2010);
Van Damme (2018); Flatt (2004), yet in the reverse case when larger particles
are added to small ones, only an increase in viscosity is reported Cwalina and
Wagner (2016); Madraki  _et al._ (2017, 2018). Extending the understanding of
constraint-controlled rheology to suspensions with size-disperse particles is
thus a key open problem. Our approach addresses this challenge in a model
system of just two species. Already in this minimal system, particle size
disparity can have a profound effect on the rheology, which, with the
exception of a few recent works Guy  _et al._ (2020); Malbranche  _et al._
(2023); Pednekar  _et al._ (2018); Monti and Rosti (2023), has remained
largely unexplored.

Figure 1:  Effect of mixing small and large particles on the rheology of a
model suspension with size ratio $\Delta=4$. Shown are (a) the relative
viscosity $\eta_{r}$ as a function of the volume ratio of small particles
$\alpha$ and the applied stress $\sigma/\sigma_{0}$, at $\phi=0.57$. Red and
green lines represent common data sets across each panel, showing,
respectively, $\eta_{r}(\alpha)$ for $\sigma/\sigma_{0}=100$, and
$\eta_{r}(\sigma)$ for $\alpha=0.25$. (b) $\eta_{r}$ as a function of
$\sigma/\sigma_{0}$ at fixed $\alpha=0.25$ for various $\phi$. (c) $\eta_{r}$
as a function of $\sigma/\sigma_{0}$ at fixed $\phi=0.57$ for various
$\alpha$. (d) $\eta_{r}$ as a function of $\alpha$ for various $\phi$ at fixed
$\sigma/\sigma_{0}=100$ (e) $\eta_{r}$ as a function of $\alpha$ for
$\phi=0.57$ for various $\sigma/\sigma_{0}$. (f)-(h) Snapshots of the
simulation at $\phi=0.57,\sigma/\sigma_{0}=100$, $\Delta=4$, and (top-to-
bottom) $\alpha=0.9,0.5,0.1$.

The rheology of bidisperse suspensions is governed by $\phi_{J}$, which is
controlled by the species size ratio $\Delta$, their volumetric mixing ratio
$\alpha=N_{S}a_{S}^{3}/(N_{S}a_{S}^{3}+N_{L}a_{L}^{3})$, and the particle
friction coefficient $\mu$ Maranzano and Wagner (2001a). Previous experimental
Maranzano and Wagner (2001a, b) and numerical Guy  _et al._ (2020); Malbranche
_et al._ (2023); Pednekar  _et al._ (2018) studies show that for a constant
solids fraction $\phi$, the relative suspension viscosity $\eta_{r}$ decreases
with $\Delta$ at fixed $\alpha$ while varying non-monotonically with $\alpha$
at fixed $\Delta$. This was explained on the basis of $\phi_{J}$ being non-
monotonic in $\alpha$, for which several models have been proposed Yu and
Standish (1991); Qi and Tanner (2011). In dry granular materials it is
understood that $\phi_{J}$ depends on $(\Delta,\alpha,\mu)$ Srivastava  _et
al._ (2021), and there have been attempts to relate this to the fraction of
smaller particles being rattlers Petit  _et al._ (2020); Srivastava  _et al._
(2021).

In this Letter, we explore systematically the role of bidispersity on the
rheology of dense, non-Brownian suspensions using an established simulation
scheme Seto  _et al._ (2013); Mari  _et al._ (2015). We explore $\Delta\leq
12$, $0<\alpha<1$ and $\phi$ close to the jamming point. We show that the non-
monotonic $\eta_{r}(\alpha)$ coincides with an exchange of stress contribution
dominance from large-large (LL) contacts at small $\alpha$ to small-small (SS)
contacts at large $\alpha$, while the stress carried by large-small (LS)
contacts is non-monotonic in $\alpha$ and vanishes at the extremes. Combining
an _ad hoc_ extension of WC with a geometric model for
$\phi_{J}(\Delta,\alpha)$ by Yu and Standish (1991) (YS), we obtain
qualitative predictions of the rheology for any bidisperse, shear thickening,
non-Brownian suspension. Using the model, we explore practical settings where
small additives are incorporated into suspensions of large particles,
rationalising the counter-intuitive finding that increasing $\phi$ can in some
circumstances reduce $\eta_{r}$. Meanwhile adding large particles to a
suspension of small ones always enhances $\eta_{r}$, corroborating
experimental findings Madraki  _et al._ (2017, 2018).

Figure 2:  Role of the particle size ratio $\Delta$. (a) Variation with
$\alpha$ of $\phi_{0}$, the jamming point at $\sigma/\sigma_{0}=0$, for a
range of $\Delta$. Solid lines with points represent simulation data; dashed
lines show model predictions Yu and Standish (1991); (b) Variation with
$\alpha$ of $\phi_{m}$, the jamming point at $\sigma/\sigma_{0}=\infty$, for a
range of $\Delta$; (c) Relative viscosity $\eta_{r}$ as a function of
$\sigma/\sigma_{0}$ at $\phi=0.57$ and $\alpha=0.5$ for various $\Delta$.
(d)-(f) Snapshots of simulations with $\phi=0.5$, $\mu=1$, $\alpha=0.5$ and
(left to right) $\Delta=2$, $3$, $4$, $6$, $12$.

## Simulations.

We model $6000$ spheres for $\Delta\leq 6$ ($12000$ for $\Delta=12$) dispersed
in density-matched Newtonian liquid under imposed shear stress $\sigma$, in a
constant volume periodic domain. The suspension flows with time-dependent
shear rate $\dot{\gamma}(t)$, and we compute the relative viscosity as
$\eta_{r}(t)=\sigma/\eta\dot{\gamma}(t)$, with $\eta$ the liquid viscosity.
Rheology data shown in the following are averages of $\eta_{r}$ over 5 strain
units at steady state, across 5 realizations. Particles interact through
short-range hydrodynamic lubrication and contact forces and torques. The
latter includes a tangential component regulated by a sliding friction
coefficient $\mu=1$ (with rolling friction Singh  _et al._ (2020) absent),
activated above a threshold normal force $F_{0}$ Seto  _et al._ (2013); Mari
_et al._ (2014) giving a characteristic stress scale $\sigma_{0}=F_{0}/a^{2}$.
We perform simulations for $\Delta=2,3,4,6,12$ and $\alpha=[0.05,0.9]$.

## Bidisperse suspension rheology.

Figure 1 shows the effect of $\alpha$ on the rheology for exemplar data with
$\Delta=4$. We present a map of $(\eta_{r},\alpha,\sigma)$ at $\phi=0.57$ in
Fig. 1(a) showing an overview of the behavior, with planar slices showing (b)
$\eta_{r}(\sigma/\sigma_{0},\phi)$ at $\alpha=0.25$; (c)
$\eta_{r}(\sigma/\sigma_{0},\alpha)$ at $\phi=0.57$; (d)
$\eta_{r}(\alpha,\phi)$ at $\sigma/\sigma_{0}=100$; (e)
$\eta_{r}(\alpha,\sigma/\sigma_{0})$ at $\phi=0.57$. In Figs. 1(f)-(h) are
snapshots of the simulation at $\alpha=0.9,0.5,0.1$. In Fig 1(b), one observes
canonical thickening behavior, qualitatively similar to quasi-monodisperse
systems showing a stress-mediated transition between two Newtonian plateaus
Mari  _et al._ (2014); Guy  _et al._ (2015); Ness and Sun (2016). The
viscosities of the plateaus increase with $\phi$ towards their respective
jamming points $\phi_{0}\equiv\phi_{J}(\sigma/\sigma_{0}\to 0)\approx 0.75$
and $\phi_{m}\equiv\phi_{J}(\sigma/\sigma_{0}\to\infty)\approx 0.65$ (these
numbers being sensitive to $\alpha$, $\Delta$ and $\mu$). In Fig. 1(c) we find
that $\sigma_{0}$ increases monotonically with $\alpha$ (since the former is
related to the particle size through $\sigma_{0}\sim F_{0}/a^{2}$), while the
frictionless and frictional viscosities measured at $\sigma/\sigma_{0}=0.01$
and $\sigma/\sigma_{0}=100$ respectively show non-monotonic dependence on
$\alpha$. The crossover of these flow curves is a result of the combined
effect of bidispersity on $\sigma_{0}$ and $\phi_{J}$ and is not predicted by
models that assume monodispersity. Next we present data in the
$\eta_{r}(\alpha)$ plane. Figure 1(d) shows $\eta_{r}$ for
$\sigma/\sigma_{0}=100$ as a function of $\alpha$ for various $\phi$. At fixed
$\phi$, $\eta_{r}$ first decreases with increasing $\alpha$, reaching a
minimum before increasing again so that
$\eta_{r}(\alpha=0)=\eta_{r}(\alpha=1)$. In the limits $\alpha=(0,1)$ the
suspension is monodisperse and exhibits identical rheology due to the size
invariance of non-Brownian suspensions (when $\sigma/\sigma_{0}$ is 0 or
$\infty$). The value of $\alpha$ that minimizes $\eta_{r}$ is insensitive to
$\phi$. For $\phi\geq 0.58$, the suspension near the extrema $\alpha=0,1$ is
jammed, and the window of $\alpha$ for which $\eta_{r}$ is finite decreases
with increasing $\phi$ so that no flow occurs at any $\alpha$ for $\phi\geq
0.65$. Figure 1(e) shows $\eta_{r}(\alpha)$ for various $\sigma/\sigma_{0}$ at
fixed $\phi=0.57$. For $\sigma/\sigma_{0}=100$, $\eta_{r}$ at the large and
small limits of $\alpha$ are equal, while for $\sigma/\sigma_{0}\leq 5$,
$\eta_{r}$ is higher for small $\alpha$ ($\alpha\to 0$) as compared to larger
values ($\alpha\to 1$). This can be explained based on the relation between
$\sigma_{0}$ and particle size $a$. Since $\sigma_{0}\sim 1/a^{2}$,
$\sigma_{0}$ increases with $\alpha$ (Fig. 1(c)) so that the jamming point is
governed by friction at lower stress when particles are large.

Figure 2 shows the effect of $\Delta$ on the jamming points $\phi_{0}$,
$\phi_{m}$, and on $\eta_{r}$. In Figs. 2(a) and (b) are the variation of
$\phi_{0}$ and $\phi_{m}$ with $\alpha$, for $\Delta=(2,3,4,6,12)$, obtained
by simulating the limits $\sigma/\sigma_{0}\to 0$ and
$\sigma/\sigma_{0}\to\infty$ for a range of $\phi$, then fitting the viscosity
as $\eta_{r}=(1-\phi/\phi_{m,0})^{-2}$ Krieger and Dougherty (1959). Both
$\phi_{0}$ and $\phi_{m}$ vary non-monotonically with $\alpha$, with the
dependence being more pronounced for increasing $\Delta$. Also plotted are
$\phi_{0,m}$ predictions as functions of $\alpha$ at various $\Delta$
following the model of Yu and Standish (1991), which produces jamming point
predictions based on geometry only. Interestingly the model works better for
the frictionless limit ($\phi_{0}$, $\sigma/\sigma_{0}=0$) as compared to the
frictional one ($\phi_{m}$, $\sigma/\sigma_{0}=\infty$), likely due to the
absence of friction and shear-induced structure in the theory. (Indeed,
understanding the disparity in the nature of jamming between frictionless and
frictional systems remains an open challenge Liu and Nagel (2010).) In the
frictionless case, the simulation data appear to be converging toward the
theory, though testing this for larger $\Delta$ rapidly becomes
computationally intractable. Plotting $\eta_{r}$ as a function of
$\sigma/\sigma_{0}$ for different $\Delta$ at fixed $\alpha=0.5,\phi=0.57$
(Fig. 2(c)), one finds the viscosity of both limiting states decreases with
increasing $\Delta$, as the proximity to jamming is decreased.

Figure 3:  Viscosity contributions from each contact type. (a) Relative
contribution $\eta_{r}^{XX}$ to the total viscosity $\eta_{r}$ of each contact
type (where SS, LS, and LL replace $XX$) as a function of $\alpha$, for
$\phi=0.57$ and $\sigma/\sigma_{0}=100$; (b)-(d) Total viscosity and its
contact contributions as a function of $\phi$ at $\sigma/\sigma_{0}=100$ for
$\alpha=$ (b) 0.15; (c) 0.5, and (d) 0.8;

## Contribution of different contact types.

We next address the microscopic underpinning of the non-monotonic
$\eta_{r}(\alpha)$ reported in Fig. 1 (see also Maranzano and Wagner (2001a);
Pednekar  _et al._ (2018); Malbranche  _et al._ (2023)), using simulation data
to separate the stress contributions of different types of contacts LL, LS,
and SS. Figure 3(a) shows the viscosity contribution of each contact type
scaled with $\eta_{r}$ for $\phi=0.57$ at $\sigma/\sigma_{0}=100$ as a
function of $\alpha$ for $\Delta=4$. At small $\alpha$, LL contacts provide
the dominant contribution to $\eta_{r}$, while SS take over at large $\alpha$
as expected. LL decay from 1 at $\alpha=0$ (by definition) to $\approx 0$ for
$\alpha\geq 0.6$ with the SS contribution being minimal for $\alpha\leq 0.2$
and increasing to 1 at $\alpha=1$. Notably the value of $\alpha$ at which SS
contributions begin to increase and LL vanishes are not symmetric with respect
to 0 and 1, and the contribution of LS is non-monotonic: by definition it
vanishes for $\alpha=0,1$, and is maximal around $\alpha=0.4$ where it
contributes $\approx~{}$75% of the overall viscosity. Generalizing these
findings for different values of $\phi$, Figs. 3(b)-(d) present $\eta_{r}$ as
a function of $\phi$ for $\alpha=$ 0.15, 0.5, and 0.8 respectively,
demonstrating that LL and SS are the dominant contributors to jamming at
$\alpha=0.15$ and 0.8, respectively, while for $\alpha=0.5$ the major
contribution to originates from LS.

Figure 4:  Constitutive model predictions. (a) Combining the WC model for
$\eta_{r}(\sigma/\sigma_{0})$ with the YS model for $\phi_{J}(\alpha,\Delta)$,
one obtains the rheology of bidisperse frictional suspensions. (b)-(d) Adding
small particles to a monodisperse large particle packing at
$\sigma/\sigma_{0}=\infty$: (b) Variation of $\phi$ and $\phi_{m}$ with
$\alpha$, and $\eta_{r}$ plotted as a function of (c) $\alpha$ and (d) $\phi$.
(e)-(g) Adding large particles to a monodisperse small particle packing at
$\sigma/\sigma_{0}=\infty$: (e) Variation of $\phi$ and $\phi_{m}$ with
$\alpha$, and $\eta_{r}$ plotted as a function of (f) $\alpha$ and (g) $\phi$.
(h)-(i) Adding small particles to a monodisperse large particle packing at
intermediate $\sigma/\sigma_{0}$, so that large particles are frictional and
small ones are frictionless. (h) Variation of $\phi$ and $\phi_{J}$ with
$\alpha$, and $\eta_{r}$ plotted as a function of (i) $\alpha$ and (j) $\phi$.
Also shown in (h) in gray are $\phi_{m}$ and $\phi_{0}$. Colors in (b)-(j)
represent different initial $\phi$; black lines represent $\phi_{m}(\alpha)$
and $\phi_{J}(\alpha,\sigma/\sigma_{0}$). Colored arrows indicate adding small
(pink) and large (green) particles.

## Constitutive model.

Given the qualitative agreement of the dependence of $\phi_{0}$ and $\phi_{m}$
on $\alpha$ and $\Delta$ with YS and of $\eta_{r}(\sigma/\sigma_{0})$ with WC,
we are motivated to construct a combined model to capture the full behavior.
Conventional WC writes
$\eta_{r}(\sigma/\sigma_{0})=(1-\phi/\phi_{J}(\sigma/\sigma_{0}))^{-2}$ with
$\phi_{J}(\sigma/\sigma_{0})=f\phi_{m}+(1-f)\phi_{0}$ and
$f=\exp(-\sigma_{0}/\sigma)$. Here we let $\phi_{0}$ and $\phi_{m}$ be
functions of $\alpha$ and $\Delta$ as plotted in Fig. 2(a)-(b), and we let
$\sigma_{0}$ interpolate in an _ad hoc_ way between small and large particle
limits as $\sigma_{0}=F_{0}/(\alpha a_{S}+(1-\alpha)a_{L})^{2}$. Doing so one
obtains a qualitative picture of the rheology, Fig. 4(a).

Using this model we can predict the effect on the rheology of compositional
changes. To demonstrate this, we consider in Fig. 4(b)-(d) the effect of
adding small particles to a suspension of monodisperse large ones, for
$\sigma/\sigma_{0}=\infty$ and $\Delta=4$. Counter-intuitively, for
sufficiently large initial $\phi$ (of large particles), the addition of small
particles leads to a decrease in $\eta_{r}$. The model predicts that the
increase in $\phi$ is slower than that of $\phi_{m}$, hence
$(1-\phi/\phi_{m})^{-2}$ decreases. Conversely, the asymmetry in
$\phi_{m}(\alpha)$ means adding large particles to a packing of small ones
leads, more intuitively, to an increase in both $\phi$ and $\eta_{r}$ (Fig.
4(e)-(g)). At intermediate $\sigma/\sigma_{0}$, $\phi_{J}$ interpolates
between frictionless and frictional values (Fig. 4(h)), and the range of
$\alpha$ for which viscosity reduction appears is broadened (Fig. 4(i)-(j)).
Recent experiments have shown that adding large non-Brownian particles to a
suspension of small ones leads to enhanced thickening behavior Madraki  _et
al._ (2018, 2017), consistent with our prediction (Fig. 4(e)) that doing so
will always move the system closer to its respective $\phi_{m}$.

## Concluding remarks.

We have shown that an existing constitutive relation for rate-dependent
rheology (WC, Wyart and Cates (2014)) can be combined in an _ad hoc_ way with
a geometric model by Yu and Standish (1991) to obtain a qualitative picture of
the rheology of bidisperse suspensions. Our ansatz for $\sigma_{0}(\alpha)$ is
not strictly in accordance with WC, in the sense that the order parameter $f$
can no longer represent the fraction of contacts that are frictional (see Ref
Guy  _et al._ (2020) for more details). Nonetheless our model predicts novelty
in the rheology that is present in reality but unaddressed by WC, specifically
that adding small particles to a system of large ones can counterintuitively
decrease the viscosity, whereas in the opposite case adding large particles to
small particles always leads to an increase in $\eta_{r}$. We show
additionally that the non-monotonic dependence of relative viscosity
$\eta_{r}$ on $\alpha$ for constant $\phi$ (as reported by Maranzano and
Wagner (2001a, b); Pednekar  _et al._ (2018); Malbranche  _et al._ (2023)) can
be understood by delineating the stress contributions of each type of contact,
and our results suggest that the rheology at $\Delta=12$ is already close to
the $\Delta\to\infty$ limit, so the predictions made here will be useful
across bidisperse systems of arbitrary size ratio. The consequences of these
results for more complex systems, especially polydisperse samples Mwasame  _et
al._ (2016) in which colloidal forces may become relevant Li  _et al._ (2023),
are broad. In particular, our results provide a direction towards the limiting
case of predicting the rheology of mixtures where larger additives are
included in a continuous background phase of much smaller particles Cwalina
and Wagner (2016). Although we have focused on the model systems of two
species, the results can be used to qualitatively predict the rheology of
polydisperse suspensions, which has many applications: industrial processing
of slurries and muds; manufacturing of amorphous solid dispersions; and
predicting the runout of geophysical flows comprising grains spanning many
orders of magnitude Breard  _et al._ (2023).

Codes and scripts necessary to reproduce the results reported in this article
are available on request. We acknowledge support from the Center for
Hierarchical Materials Design (CHiMaD) under award number 70NANB19H005 (US
Dept. Commerce), and from the Army Research Office under grants
W911NF-19-1-0245, W911NF-21-2-0146, and W911NF-21-1-0038. A. S. acknowledges
Case Western Reserve University for start-up funding. C.N. thanks Eric Breard
for pointing to Ref Yu and Standish (1991), and acknowledges support from the
Royal Academy of Engineering under the Research Fellowship scheme and from the
Leverhulme Trust under Research Project Grant RPG-2022-095. Part of this work
made use of the High Performance Computing Resource in the Core Facility for
Advanced Research Computing at Case Western Reserve University.

A.S. and C.N. contributed equally to this work.

## References

  * Ness  _et al._ (2022) C. Ness, R. Seto, and R. Mari, Annual Review of Condensed Matter Physics 13, 97 (2022).
  * Morris (2020) J. F. Morris, Annual Review of Fluid Mechanics 52, 121 (2020).
  * Singh (2023) A. Singh, MRS Communications , 1 (2023).
  * Roussel  _et al._ (2010) N. Roussel, A. Lemaître, R. J. Flatt, and P. Coussot, Cement and Concrete Research 40, 77 (2010).
  * Brown and Jaeger (2014) E. Brown and H. M. Jaeger, Reports on Progress in Physics 77, 046602 (2014).
  * Peters  _et al._ (2016) I. R. Peters, S. Majumdar, and H. M. Jaeger, Nature 532, 214 (2016).
  * Singh  _et al._ (2019) A. Singh, S. Pednekar, J. Chun, M. M. Denn, and J. F. Morris, Physical Review Letters 122, 098004 (2019).
  * Guy  _et al._ (2018) B. M. Guy, J. Richards, D. Hodgson, E. Blanco, and W. C. K. Poon, Physical Review Letters 121, 128001 (2018).
  * Singh  _et al._ (2020) A. Singh, C. Ness, R. Seto, J. J. de Pablo, and H. M. Jaeger, Physical Review Letters 124, 248005 (2020).
  * Singh  _et al._ (2022) A. Singh, G. L. Jackson, M. van der Naald, J. J. de Pablo, and H. M. Jaeger, Physical Review Fluids 7, 054302 (2022).
  * Seto  _et al._ (2013) R. Seto, R. Mari, J. F. Morris, and M. M. Denn, Physical Review Letters 111, 218301 (2013).
  * Wyart and Cates (2014) M. Wyart and M. E. Cates, Physical Review Letters 112, 098302 (2014).
  * Singh  _et al._ (2018) A. Singh, R. Mari, M. M. Denn, and J. F. Morris, Journal of Rheology 62, 457 (2018).
  * Guy  _et al._ (2015) B. M. Guy, M. Hermes, and W. C. K. Poon, Physical Review Letters 115, 088304 (2015).
  * Ness and Sun (2016) C. Ness and J. Sun, Soft Matter 12, 914 (2016).
  * Guy  _et al._ (2020) B. M. Guy, C. Ness, M. Hermes, L. J. Sawiak, J. Sun, and W. C. Poon, Soft Matter 16, 229 (2020).
  * Pradeep  _et al._ (2021) S. Pradeep, M. Nabizadeh, A. R. Jacob, S. Jamali, and L. C. Hsiao, Physical Review Letters 127, 158002 (2021).
  * Van Damme (2018) H. Van Damme, Cement and Concrete Research 112, 5 (2018).
  * Flatt (2004) R. Flatt, Materials and Structures 37, 289 (2004).
  * Cwalina and Wagner (2016) C. D. Cwalina and N. J. Wagner, Journal of Rheology 60, 47 (2016).
  * Madraki  _et al._ (2017) Y. Madraki, S. Hormozi, G. Ovarlez, E. Guazzelli, and O. Pouliquen, Physical Review Fluids 2, 033301 (2017).
  * Madraki  _et al._ (2018) Y. Madraki, G. Ovarlez, and S. Hormozi, Physical Review Letters 121, 108001 (2018).
  * Malbranche  _et al._ (2023) N. Malbranche, B. Chakraborty, and J. F. Morris, Journal of Rheology 67, 91 (2023).
  * Pednekar  _et al._ (2018) S. Pednekar, J. Chun, and J. F. Morris, Journal of Rheology 62, 513 (2018).
  * Monti and Rosti (2023) A. Monti and M. E. Rosti, Meccanica 58, 727 (2023).
  * Maranzano and Wagner (2001a) B. J. Maranzano and N. J. Wagner, Journal of Chemical Physics 114, 10514 (2001a).
  * Maranzano and Wagner (2001b) B. J. Maranzano and N. J. Wagner, Journal of Rheology 45, 1205 (2001b).
  * Yu and Standish (1991) A.-B. Yu and N. Standish, Industrial & Engineering Chemistry Research 30, 1372 (1991).
  * Qi and Tanner (2011) F. Qi and R. I. Tanner, Korea-Australia Rheology Journal 23, 105 (2011).
  * Srivastava  _et al._ (2021) I. Srivastava, S. A. Roberts, J. T. Clemmer, L. E. Silbert, J. B. Lechman, and G. S. Grest, Physical Review Research 3, L032042 (2021).
  * Petit  _et al._ (2020) J. C. Petit, N. Kumar, S. Luding, and M. Sperl, Physical Review Letters 125, 215501 (2020).
  * Mari  _et al._ (2015) R. Mari, R. Seto, J. F. Morris, and M. M. Denn, Physical Review E 91, 052302 (2015).
  * Mari  _et al._ (2014) R. Mari, R. Seto, J. F. Morris, and M. M. Denn, Journal of Rheology 58, 1693 (2014).
  * Krieger and Dougherty (1959) I. M. Krieger and T. J. Dougherty, Transactions of the Society of Rheology 3, 137 (1959).
  * Liu and Nagel (2010) A. J. Liu and S. R. Nagel, Annual Review of Condensed Matter Physics 1, 347 (2010).
  * Mwasame  _et al._ (2016) P. M. Mwasame, N. J. Wagner, and A. N. Beris, Journal of Rheology 60, 225 (2016).
  * Li  _et al._ (2023) X. Li, J. R. Royer, and C. Ness, arXiv preprint arXiv:2307.13802 (2023).
  * Breard  _et al._ (2023) E. C. Breard, J. Dufek, S. Charbonnier, V. Gueugneau, T. Giachetti, and B. Walsh, Nature Communications 14, 2079 (2023).


# A genealogical interpretation of fixed and random effects models of complex
traits

Hanbin Lee†,‡ †Department of Medicine, Seoul National University College of
Medicine
Seoul, Republic of Korea
‡Department of Mathematical Sciences, Seoul National University
Seoul, Republic of Korea<EMAIL_ADDRESS>
###### Abstract.

In the analysis of complex traits, genetic effects are frequently modelled as
either fixed or random effects. Such assumptions serve as a foundation of
defining heritability and relatedness using genome-wide single nucleotide
polymorphism (SNP) markers. In this work, I propose a genealogical framework
connecting the two assumptions conditional on the ancestral recombination
graph (ARG). It turns out that the reference time point in which the
probability is defined determines whether the effect of a variant to behave as
either a fixed or random effect. This lays a connection between the PC
regression and linear mixed model (LMM) for genetic association study. The
framework induces a genetic relatedness matrix (GRM) in which the elements are
a function of the time to the recent common ancestor. Subsequently, a novel
trait variance decomposition respect to the regions sharing a common genealogy
is followed. The variance decomposition then provides a natural means to
define heritability and a potential method for gene mapping.

##  1\. Introduction

The polygenic model of complex traits models the trait as a linearly additive
function of genetic variants. The model dates back to Fisher and played a
pivotal role in combining Mendelianism and Darwinism [Fisher, 1919]. The model
is now adopted in a wide range of applications including heritability
estimation [Yang et al., 2010], genome-wide associations study (GWAS) [Price
et al., 2006, Kang et al., 2010, Zhou and Stephens, 2012] and summary
statistics based methods [Bulik-Sullivan et al., 2015, Shi et al., 2016], to
name a few.

Although the functional form of the polygenic model is shared across studies,
additional modelling assumptions are frequently made. Among such assumptions,
random effects models of genetic effects is commonly adopted. The genetic
random effects model was first introduced in breeding literature by Henderson
[Henderson, 1974], and was later introduced to human genetics [Yang et al.,
2010]. Linear mixed model (LMM), which includes both fixed and random effects,
was independently adopted to gene mapping in model organisms, and subsequently
in humans [Kang et al., 2008, Kang et al., 2010].

The unique feature of genetic random effects is that the randomness is for the
features (or the columns) and not the samples (or the rows). Most applications
of random effects outside genetics employ random effects to model the row-wise
variability, for example, nested samples within groups [Gelman, 2006]. In such
a case, the source of randomness can be attributed to sampling variation.
However, it cannot explain the column-wise random effects in genetic
applications.

In this work, I show that the column-wise variability of genetic random
effects can be attributed to randomly occurring mutations in the genome. For a
given reference time point in the past [Weir and Cockerham, 1984, Thompson,
2013], polymorphic loci appears randomly through mutations which can be
thought as sampling trait-associated sites randomly. This serves as the basis
of column-wise randomness of genetic random effects. In turn, mutations that
have emerged prior to the reference time point behave as fixed effects. This
observation suggests a putative connection between PC regression and LMM used
in genetic association studies.

The resulting genetic relatedness matrix (GRM) closely resembles the version
of McVean [McVean, 2009] which was first developed to interpret principal
component analysis (PCA) applied to the genotype matrix. Under neutral
evolution, the elements of the GRM is a function of pairwise coalescence time
which is a natural measure of relatedness. The trait variance can be
decomposed into local regions sharing the same genealogical history that gives
a natural way of defining local heritability and a potential method to perform
quantitative trait loci (QTL) mapping.

##  2\. Results

###  2.1. The ancestral recombination graph

The ancestral recombination graph (ARG) is a collection of evolving local
coalescent trees as we navigate through the genome [Hudson, 1992, Griffiths
and Marjoram, 1996]. The tree at a fixed position is called the marginal tree.
Nearby trees are correlated due to recombination events. Let
$l=1,\ldots,n_{l}$ be the index of loci across the genome. $n_{l}$ is the
number of regions with different marginal tree topology. Sites indexed by $s$
are nested within a locus ($s\in l$) and we denote the number of sites within
locus $l$ as $n^{l}$.

The marginal tree at locus $l$ is denoted as $T^{l}$. We index the branches of
a tree as $b$ and write its length as $L(b)$. The total length of $T^{l}$ is
$L^{l}=\sum_{b\in T^{l}}L(b)$. We override the the notation to write $b\in
T^{l}$ which means the the branch is within $T^{l}$. The mutation rate of a
locus is denoted as $u^{l}$. As there are $n^{l}$ sites within the locus, the
per-site mutation rate is given by $u^{l}/n^{l}$. The number of total
mutations on $T^{l}$ follows a Poisson distribution with mean $u^{l}\sum_{b\in
T^{l}}L(b)=u^{l}L^{l}$. This implicitly assumes that the mutation is neutral
so that the mutation does not alter the topology of the tree.

We make the infinite sites assumption so that every time a mutation occurs, it
appears at a new site [Kimura, 1969]. Therefore, back mutation is ignored. In
this work, the ARG is given fixed and mutations occur at random. The general
theory of this setting has been studied by Ralph [Ralph, 2019].

###  2.2. An interpretation of the genetic random effects

Let $i$ be the index of haploids and $s$ be the index of sites. We write the
row indices as subscripts and the column indices as superscripts.

Let $\mathbb{I}^{s}$ be the indicator function of site $s$ being polymorphic.
The total number of sites is denoted as $N_{s}$. Assuming the infinite sites
model, a mutation occurs only once at a site with probability proportional to
the total length of the tree at the locus [Kimura, 1969]. $G_{i}^{s}=0,1$ is
the genetic value of individual $i$ at site $s$. Then the polygenic model can
be written as

$\displaystyle y_{i}=\sum_{l}\sum_{s\in
l}G_{i}^{s}\beta^{s}\mathbb{I}^{s}+\beta^{0}+\epsilon_{i}$ (1)

where $\beta^{0}$ is the intercept and $\epsilon_{i}$ is the random error.

We consider $n$ samples. The trait vector, the genotype matrix, the error
vector of the sample are $\mathbf{y}=(y_{1},\ldots,y_{n})^{T}$,
$\mathbf{G}^{s}=(G_{1}^{s},\ldots,G_{n}^{s})^{T}$ and
$\boldsymbol{\epsilon}=(\epsilon_{1},\ldots,\epsilon_{n})^{T}$. Equation 1 can
be written as

$\displaystyle\mathbf{y}=\sum_{l}\sum_{s\in
l}\mathbf{G}^{s}\beta^{s}\mathbb{I}^{s}+\beta^{0}\mathbf{1}_{n}+\boldsymbol{\epsilon}$
(2)

where $\mathbf{1}_{n}=(\overbrace{1,\ldots,1}^{n\text{ copies}}$$)^{T}$.

Here, each $\beta^{s}$ are assumed to be constant as the variant at a fixed
position in the genome is likely to have the same effect. Instead,
$\mathbb{I}^{s}$ and $G_{i}^{s}$ are random variables respect to the
mutational process emerging on the ARG. Each evolutionary trial beginning at a
reference time point from the past, we observe a different set of variants
corresponding to their effect size because mutations occur at random
positions. The distribution function of this random variable $\beta$ is
therefore

$\displaystyle\begin{aligned} \mathbb{P}\left(\beta\leq
x\right)&=\sum_{s:\beta^{s}\leq x}\mathbb{P}\left(\mathbb{I}^{s}=1\right)\\\
&=\sum_{l}\sum_{s:\beta^{s}\leq x,s\in
l}\frac{u^{l}L^{l}}{n^{l}}e^{-u^{l}L^{l}/n^{l}}\\\
&\approx\sum_{l}u^{l}L^{l}\cdot\frac{\\#\\{\text{sites $s$ such that
$\beta^{s}\leq x$}\\}}{n^{l}}\end{aligned}$ (3)

The effect size distribution is often termed as the genetic architecture
[Timpson et al., 2017]. According to Equation 3, the genetic architecture of a
trait depends on the mutation rate and the topology of the ARG.

###  2.3. The reference time point

Let the current time $t_{0}=0$ and the time to the most recent common ancestor
of the population as $-t_{\mathrm{MRCA}}<0$.
$t_{\mathrm{ref}}\in[-t_{\mathrm{MRCA}},t_{0}]$ is the reference time point.
Time to MRCA can be different across loci. We add a superscript to denote this
difference $t_{\mathrm{MRCA}}^{l}$. The value is constant within a locus
sharing the local tree.

At $t_{\mathrm{MRCA}}$, all sites are monomorphic so $\mathbb{I}^{s}=0$ for
all $s$. Since mutations accumulate over time, the number of polymorphic sites
$n_{s}$ increases with time. Therefore, $\mathbb{I}^{s}$ increases over time,
and is a monotone stochastic process on $[-t_{\mathrm{MRCA}},t_{0}]$. I denote
the associated filteration as $\mathcal{F}(t)$. Then the conditional
expectation is defined as
$\mathbb{E}_{t}\left[\;\cdot\;\right]:=\mathbb{E}\left[\;\cdot\mid\mathcal{F}(t)\right]$.
Once the mutation has arrived at time $t$, $\mathbb{I}^{s}$
($s=1,\ldots,N_{s}$) and $G_{i}^{s}$ ($i=1,\ldots,N$) remain constant.

###  2.4. The genetic relatedness matrix and the variance decomposition

We now assume that $G_{i}^{s}$ is independent from $\epsilon_{i}$. Since
$G_{i}^{s}$ is determined by the mutational process holding the ARG fixed,
this means the mutational process is independent from $\epsilon_{i}$. Let
$S(t)$ be the set of sites in which mutations have arrived before time $t$.
Also, $h_{t}^{l}=\mathrm{min}\\{t_{\mathrm{MRCA}}^{l},t\\}$ and
$h_{t,ij}^{l}=\mathrm{min}\\{t_{\mathrm{MRCA},ij}^{l},t\\}$ where
$t_{\mathrm{MRCA},ij}^{l}$ is the coalescence time of individual $i$ and $j$
at locus $l$.

######  Theorem 1.

The trait variance respect to the reference time point $t_{\mathrm{ref}}$ is

$\displaystyle\mathrm{Var}_{t_{\mathrm{ref}}}\left(\mathbf{y}\right)=\sum_{l}(\sigma_{t_{\mathrm{ref}}}^{l})^{2}\mathbf{A}_{t_{\mathrm{ref}}}^{l}+\sigma_{e}^{2}\mathbf{I}_{n\times
n}$

where

$\displaystyle\begin{gathered}\left(\sigma_{t}^{l}\right)^{2}=\frac{u^{l}}{n_{l}}\sum_{s\notin
S(t)\&s\in l}\left(\beta^{s}\right)^{2}\\\
\left[\mathbf{A}_{t}^{l}\right]_{ij}=h_{t}^{l}-h_{t,ij}^{l}\end{gathered}$

Note that the diagonal elements of $\mathbf{A}_{t}^{l}$ are
$h_{t}^{l}-h_{t,ij}^{l}=h_{t}^{l}$ because $h_{t,ij}^{l}=0$ for $i=j$.
Therefore, one might choose to normalize the matrix with $h_{t}^{l}$ to write
as

$\displaystyle\begin{gathered}\left(\widetilde{\sigma}_{t}^{l}\right)^{2}=\frac{h_{t}^{l}u^{l}}{n_{l}}\sum_{s\notin
S(t)\&s\in l}\left(\beta^{s}\right)^{2}\\\
\left[\widetilde{\mathbf{A}}_{t}^{l}\right]_{ij}=\frac{h_{t}^{l}-h_{t,ij}^{l}}{h_{t}^{l}}\end{gathered}$

Then the overall and the local heritability can be defined as

$\displaystyle\begin{gathered}h_{t}^{2}=\frac{\sum_{m}\left(\widetilde{\sigma}_{t_{\mathrm{ref}}}^{m}\right)^{2}}{\sum_{m}\left(\widetilde{\sigma}_{t_{\mathrm{ref}}}^{m}\right)^{2}+\sigma_{e}^{2}}\\\
\left(h_{t}^{l}\right)^{2}=\frac{\left(\widetilde{\sigma}_{t_{\mathrm{ref}}}^{l}\right)^{2}}{\sum_{m}\left(\widetilde{\sigma}_{t_{\mathrm{ref}}}^{m}\right)^{2}+\sigma_{e}^{2}}\end{gathered}$
(6)

because the diagonal elements are set to $1$.

Theorem 1 gives a variance decomposition respect to genetic regions $l$
sharing the same tree topology $T^{l}$. As $\sigma_{t}^{l}$ is proportional to
the sum of squares of $\beta^{s}$ within the region, a non-zero value
indicates a QTL. $\sigma_{t}^{l}$ does not depend on the topology of the tree.
Therefore, the value is the same to the one that would have been obtained from
the whole population. However, the normalized version
$\widetilde{\sigma}_{t}^{l}$ depends on $h_{t}^{l}$ which depends on the tree
topology of the sample. Therefore, the heritability in Equation 6 should be
understood as a sample-specific quantity.

###  2.5. The fixed-random transition of genetic effects

Let $S(t)$ be the set of sites where mutations have arrived at time $t$. After
the arrival of the mutation, $\mathbb{I}^{s}$ and $G_{i}^{s}$ all behave as
constants. $\mathbb{I}^{s}=1$ and $G_{i}^{s}$ will take either $0$ or $1$
based on the position of the mutation on the local tree.

######  Theorem 2.

The trait mean respect to the reference time point $t_{\mathrm{ref}}$ is

$\displaystyle\mathbb{E}_{t_{\mathrm{ref}}}\left[y_{i}\right]=\sum_{l}\sum_{s\in
l,S(t_{\mathrm{ref}})}\beta^{s}G_{i}^{s}+\sum_{l}\left(\frac{u^{l}h_{t_{\mathrm{ref}}}^{l}}{n^{l}}\cdot\sum_{s\in
l,s\notin S(t_{\mathrm{ref}})}\beta^{s}\right)+\beta^{0}$

which also can be written in vector notation as

$\displaystyle\mathbb{E}_{t_{\mathrm{ref}}}\left[\mathbf{y}\right]=\sum_{l}\sum_{s\in
l,S(t_{\mathrm{ref}})}\beta^{s}\mathbf{G}^{s}+\sum_{l}\left(\frac{u^{l}h_{t_{\mathrm{ref}}}^{l}}{n^{l}}\cdot\sum_{s\in
l,s\notin
S(t_{\mathrm{ref}})}\beta^{s}\right)\mathbf{1}_{n}+\beta^{0}\mathbf{1}_{n}$

Looking at Theorem 1 and Theorem 2, at $t_{\mathrm{ref}}$, we can see that the
time when the site became polymorphic determines whether the genetic effect of
the variant at the site behaves as either a fixed or a random effect. As one
moves $t_{\mathrm{ref}}$ from the past to the present, genetic effects evolve
gradually from random effects to fixed effects.

Sites that have become polymorphic prior to $t_{\mathrm{ref}}$ appears in the
expectation as seen in Theorem 2. However, it does not contribute to the
pairwise relatedness in the GRM as shown in Theorem 1. The exact opposite
holds for variants that are yet polymorphic at time $t_{\mathrm{ref}}$.
Although these variants do appear in the second term of Theorem 2, the value
is constant for all individuals so it has no discernible effect.

This result reveals a putative connection between PC regression and LMM. PC
regression can be thought as a penalized regression for high-dimensional
problems [Hastie et al., 2009]. The regression equation including all sites in
$S(t_{\mathrm{ref}})$ fulfills the condition. In this perspective, PC
regression is a penalized approximation of the regression in Theorem 2.
Therefore, a pure fixed-effects PC regression can be viewed as a special case
at $t_{\mathrm{ref}}=t_{0}$ and the pure random-effects LMM can be viewed as a
special case at
$t_{\mathrm{ref}}=\max_{l=1,\ldots,n_{l}}\\{t_{\mathrm{MRCA}}^{l}\\}$, the
grand TMCRA.

###  2.6. Quantitative trait loci (QTL) mapping

The variance decomposition of Theorem 1 gives a natural mean to perform
quantitative trait loci (QTL) mapping with arbitrary range up to the
resolution given by the ARG. As the true ARG is unknown, the working method is
based on the inferred ARG [Speidel et al., 2019, Kelleher et al., 2019, Wohns
et al., 2022].

The likelihood estimator we fit is based on

$\displaystyle\mathbf{y}\sim\mathcal{N}\left(\sum_{l}\sum_{s\in
l,S(t_{\mathrm{ref}})}\beta^{s}\mathbf{G}^{s}+\beta_{t_{\mathrm{ref}}}^{0}\mathbf{1}_{n},\sum_{l}(\sigma_{t_{\mathrm{ref}}}^{l})^{2}\mathbf{A}_{t_{\mathrm{ref}}}^{l}+\sigma_{e}^{2}\mathbf{I}_{n\times
n}\right)$ (7)

where $\beta_{t}^{0}=\sum_{l}\left(\frac{u^{l}h_{t}^{l}}{n^{l}}\cdot\sum_{s\in
l,s\notin S(t)}\beta^{s}\right)+\mathbf{1}_{n}$. Setting $t_{\mathrm{ref}}$ to
the grand TMRCA $\max_{l=1,\ldots,n_{l}}\\{t_{\mathrm{MRCA}}^{l}\\}$ will
measure the effect of all variants that emerged since the grand TMRCA, and
moving $t_{\mathrm{ref}}$ closer to the present will capture more recent
variants. Let
$\mathbf{V}_{n}=\sum_{l}(\sigma_{t_{\mathrm{ref}}}^{l})^{2}\mathbf{A}_{t_{\mathrm{ref}}}^{l}+\sigma_{e}^{2}\mathbf{I}_{n\times
n}$. Then the minus log-likelihood according to Equation 7 is

$\displaystyle\frac{1}{2}\log|\mathbf{V}_{n}|+\frac{1}{2}\left(\mathbf{y}-\sum_{l}\sum_{s\in
l,S(t_{\mathrm{ref}})}\beta^{s}\mathbf{G}^{s}-\beta_{t_{\mathrm{ref}}}^{0}\mathbf{1}_{n}\right)^{T}\mathbf{V}_{n}^{-1}\left(\mathbf{y}-\sum_{l}\sum_{s\in
l,S(t_{\mathrm{ref}})}\beta^{s}\mathbf{G}^{s}-\beta_{t_{\mathrm{ref}}}^{0}\mathbf{1}_{n}\right)$
(8)

The goal is to find the parameter minimizing the above objective function. A
non-zero $\sigma_{t_{\mathrm{ref}}}^{l}$ can be thought as an evidence for a
genetic signal at region $l$.

A critical drawback of estimating Equation 7 from real data is that the number
of loci is far larger than the number of samples available. A penalization
would be required to properly estimate the parameters [Tibshirani, 1996,
Hastie et al., 2009]. Sparse penalization would have an additional advantage
of loci selection.

##  3\. Discussion

This work presents a genealogical interpretation of the genetic fixed and
random effects model. Under the neutral evolution assumption, the sample
variance can be expressed in terms of local GRMs and the associated variances.
This gives a natural way to define time-specific genealogy-based heritability
and a mean to detect QTLs using genealogies. Within the framework, genetic
effects transform from random to fixed effects as the reference time point
moves closer to the present.

In this work, the ARG is given fixed and mutations occur uniformly on the
branches of the ARG following a Poisson process [Ralph, 2019]. As a result,
genetic variation is dominated by the mutational process. The remaining
variability comes from the non-genetic factor which was assumed to be
independent from the mutational process. To place mutations uniformly on a
fixed ARG, it requires the evolutionary process to be neutral. Since non-
neutral mutations subject to selection alter the topology of the tree, the
distribution of mutations depends on the position at the tree under selection.
This will subsequently change the GRM derived in this work.

Similar GRMs have been studied to date. The GRM of this study resembles that
of McVean in a sense that the genotypes were not normalized [McVean, 2009].
The core difference is that the McVean’s GRM allows the genealogy to vary. As
a result, the elements appear as a function of expected branch lengths rather
than the branch lengths themselves. My GRM conditions on the ARG so the branch
lengths are given constant which does not require taking expectations. Other
approaches have conditioned on ARG resulting GRMs as a direct function of
branch lengths rather then their expectations [Wang et al., 2021, Zhang et
al., 2021, Fan et al., 2022, Link et al., 2023]. The difference is that
genotypes are normalized with an associated factor of $\alpha$ in these works
[Speed and Balding, 2014]. My GRM is identical to the $\alpha=0$ case up to a
multiplicative constant.

This work suggests a novel form of narrow sense heritability. The advantage of
the novel definition is that it naturally breaks down into local regions
without making assumptions about the linkage disequilibrium (LD). As mutation
at different sites are assumed to be independent, conditioning the ARG removes
the correlation between sites. Therefore, LD can be ignored which gives local
genetic variances that simply adds up to the overall genetic variance. This is
in contrast with previously proposed SNP-based heritability where the total
genetic variance is not an additive function of local variances due to LD
[Yang et al., 2010, Bulik-Sullivan et al., 2015, Shi et al., 2016].

Although the paper proposes a potential method for gene mapping, the model
requires a method to estimate variance components of a high-dimensional linear
model making empirical implementation difficult. Penalized variance component
methods have been proposed but does not scale to biobank-scale data [Schaid et
al., 2020, Kim et al., 2021]. An alternative strategy would be testing local
GRMs marginally as proposed by Link and colleagues [Link et al., 2023].
However, the current theory provides very little information on the
consequence of testing local GRMs separately rather than as a whole as
suggested by the theory. Further work is warranted to understand the behavior
of marginal testing within the framework.

One notable contribution is that it hints a connection between PC regression
and LMM. Existing theoretical analyses have largely focused on the statistical
perspectives such as the number of PCs being included [Hoffman, 2013, Zhang
and Pan, 2014]. This work provides an orthogonal viewpoint based on population
genetics. PC regression accounts genetic effects through fixed effects while
LMM does so by random effects. Under my framework, a genetic effect is either
fixed or random given $t_{\mathrm{ref}}$ so PC regression and LMM can be
thought as two extremes at the same continuum. This might explain hybrid
methods that includes both PC and GRMs [Loh et al., 2015, Jiang et al., 2019].

This work has several limitations. First, the theory is based on the neutral
evolution assumption. As background and balancing selection have been found to
be pervasive, further work is required to incorporate more realistic
evolutionary scenarios [Zeng et al., 2018, O'Connor et al., 2019, Simons et
al., 2022, Veller and Coop, 2023]. Second, a realistic algorithm is required
to estimate the heritability decomposition proposed by the theory. Third, the
GRM proposed by the theory does not exactly match the ones that are used in
practice. Therefore, the analysis of previous methods in the light of the
current theory should not be taken as decisive. Finally, the model considered
is purely additive so non-linear effects such as gene-environment interaction
are ignored.

## Appendix: Proofs

The following proofs are based on the quantitative trait model and I restate
the model for convenience.

$\displaystyle\mathbf{y}=\sum_{l}\sum_{s\in
l}\mathbf{G}^{s}\beta^{s}\mathbb{I}^{s}+\beta^{0}\mathbf{1}_{n}+\boldsymbol{\epsilon}$
(A1)

###### Proof of Theorem 1.

First, expand the variance.

$\displaystyle\begin{aligned}
\mathrm{Var}_{t}\left(\mathbf{y}\right)&=\mathrm{Var}_{t}\left(\sum_{l}\sum_{s\in
l}\mathbf{G}^{s}\beta^{s}\mathbb{I}^{s}+\beta^{0}\mathbf{1}_{n}+\boldsymbol{\epsilon}\right)\\\
&=\sum_{l}\sum_{s\in
l}\mathrm{Var}_{t}\left(\mathbf{G}^{s}\beta^{s}\mathbb{I}^{s}\right)+\mathrm{Var}_{t}\left(\boldsymbol{\epsilon}\right)\\\
&=\sum_{l}\sum_{s\in
l}\left(\beta^{s}\right)^{2}\mathrm{Var}_{t}\left(\mathbf{G}^{s}\mathbb{I}^{s}\right)+\mathrm{Var}_{t}\left(\boldsymbol{\epsilon}\right)\\\
\end{aligned}$ (A2)

The second line follows from the independence of the mutational process across
sites.

Next, we compute each element of
$\mathrm{Var}_{t}\left(\mathbf{G}^{s}\mathbb{I}^{s}\right)$.

$\displaystyle\begin{aligned}
\mathrm{Cov}_{t}\left(G_{i}^{s}\mathbb{I}^{s},G_{j}^{s}\mathbb{I}^{s}\right)&=\mathbb{E}_{t}\left[G_{i}^{s}G_{j}^{s}\mathbb{I}^{s}\right]-\mathbb{E}_{t}\left[G_{i}^{s}\mathbb{I}^{s}\right]\mathbb{E}_{t}\left[G_{j}^{s}\mathbb{I}^{s}\right]\\\
&=\frac{h_{t}^{l}-h_{t,ij}^{s}}{L^{l}}\mathbb{E}_{t}\left[\mathbb{I}^{s}\right]-\frac{h_{t}^{l}}{L^{l}}\mathbb{E}_{t}\left[\mathbb{I}^{s}\right]\cdot\frac{h_{t}^{l}}{L^{l}}\mathbb{E}_{t}\left[\mathbb{I}^{s}\right]\\\
&=\frac{h_{t}^{l}-h_{t,ij}^{s}}{L^{l}}\cdot\frac{u^{l}}{n^{l}}L^{l}e^{-u^{l}L^{l}/n^{l}}-\frac{h_{t}^{l}}{L^{l}}\cdot\frac{u^{l}}{n^{l}}L^{l}e^{-u^{l}L^{l}/n^{l}}\cdot\frac{h_{t}^{l}}{L^{l}}\cdot\frac{u^{l}}{n^{l}}L^{l}e^{-u^{l}L^{l}/n^{l}}\\\
&\approx(h_{t}^{l}-h_{t,ij}^{s})\cdot\frac{u^{l}}{n^{l}}-h_{t}^{l}\cdot\frac{u^{l}}{n^{l}}\cdot
h_{t}^{l}\cdot\frac{u^{l}}{n^{l}}\\\
&\approx(h_{t}^{l}-h_{t,ij}^{s})\cdot\frac{u^{l}}{n^{l}}\quad(\because\mbox{$u^{l}/n^{l}$
is small})\\\
&=(h_{t}^{l}-h_{t,ij}^{l})\cdot\frac{u^{l}}{n^{l}}\quad(\because\mbox{the
topology is same within $l$})\end{aligned}$ (A3)

Substituting Equation A3 to Equation A2 gives the desired result.

∎

###### Proof of Theorem 2.

First, expand the expectation.

$\displaystyle\begin{aligned}
\mathbb{E}_{t}\left[\mathbf{y}\right]&=\mathbb{E}_{t}\left[\sum_{l}\sum_{s\in
l}\mathbf{G}^{s}\beta^{s}\mathbb{I}^{s}+\beta^{0}\mathbf{1}_{n}+\boldsymbol{\epsilon}\right]\\\
&=\sum_{l}\sum_{s\in
l}\mathbb{E}_{t}\left[\mathbf{G}^{s}\beta^{s}\mathbb{I}^{s}\right]+\beta^{0}\mathbf{1}_{n}\\\
&=\sum_{l}\sum_{s\in l,s\in
S(t)}\mathbb{E}_{t}\left[\mathbf{G}^{s}\beta^{s}\mathbb{I}^{s}\right]+\sum_{l}\sum_{s\in
l,s\notin
S(t)}\mathbb{E}_{t}\left[\mathbf{G}^{s}\beta^{s}\mathbb{I}^{s}\right]+\beta^{0}\mathbf{1}_{n}\\\
&=\sum_{l}\sum_{s\in l,s\in
S(t)}\beta^{s}\mathbb{E}_{t}\left[\mathbf{G}^{s}\beta^{s}\mathbb{I}^{s}\right]+\sum_{l}\sum_{s\in
l,s\notin
S(t)}\beta^{s}\mathbb{E}_{t}\left[\mathbf{G}^{s}\mathbb{I}^{s}\right]+\beta^{0}\mathbf{1}_{n}\\\
&=\sum_{l}\sum_{s\in l,s\in S(t)}\beta^{s}\mathbf{G}^{s}+\sum_{l}\sum_{s\in
l,s\notin
S(t)}\beta^{s}\mathbb{E}_{t}\left[\mathbf{G}^{s}\mathbb{I}^{s}\right]+\beta^{0}\mathbf{1}_{n}\\\
\end{aligned}$ (A4)

∎

How we compute each element of
$\mathbb{E}_{t}\left[\mathbf{G}^{s}\mathbb{I}^{s}\right]$ for $s\notin S(t)$.

$\displaystyle\begin{aligned}
\mathbb{E}_{t}\left[G_{i}^{s}\mathbb{I}^{s}\right]&=\frac{h_{t}^{l}}{L^{l}}\mathbb{E}_{t}\left[\mathbb{I}^{s}\right]\\\
&=\frac{h_{t}^{l}}{L^{l}}\cdot\frac{u^{l}}{n^{l}}L^{l}e^{-u^{l}L^{l}/n^{l}}\\\
&\approx\frac{h_{t}^{l}u^{l}}{n^{l}}\quad(\because\mbox{$u^{l}/n^{l}$ is
small})\end{aligned}$ (A5)

Substituting Equation A4 to Equation A5 gives the desired result.

## References

  * [Bulik-Sullivan et al., 2015] Bulik-Sullivan, B. K., , Loh, P.-R., Finucane, H. K., Ripke, S., Yang, J., Patterson, N., Daly, M. J., Price, A. L., and Neale, B. M. (2015).  LD score regression distinguishes confounding from polygenicity in genome-wide association studies.  Nat Genet, 47(3):291–295.
  * [Fan et al., 2022] Fan, C., Mancuso, N., and Chiang, C. W. (2022).  A genealogical estimate of genetic relationships.  The American Journal of Human Genetics, 109(5):812–824.
  * [Fisher, 1919] Fisher, R. A. (1919).  XV.—the correlation between relatives on the supposition of mendelian inheritance.  Trans. R. Soc. Edinb., 52(2):399–433.
  * [Gelman, 2006] Gelman, A. (2006).  Multilevel (hierarchical) modeling: What it can and cannot do.  Technometrics, 48(3):432–435.
  * [Griffiths and Marjoram, 1996] Griffiths, R. and Marjoram, P. (1996).  Ancestral inference from samples of DNA sequences with recombination.  Journal of Computational Biology, 3(4):479–502.
  * [Hastie et al., 2009] Hastie, T., Tibshirani, R., and Friedman, J. (2009).  The Elements of Statistical Learning.  Springer New York.
  * [Henderson, 1974] Henderson, C. (1974).  General flexibility of linear model techniques for sire evaluation.  Journal of Dairy Science, 57(8):963–972.
  * [Hoffman, 2013] Hoffman, G. E. (2013).  Correcting for population structure and kinship using the linear mixed model: Theory and extensions.  PLoS ONE, 8(10):e75707.
  * [Hudson, 1992] Hudson, R. R. (1992).  Gene genealogies and the coalescent process, volume 7 of Oxford Surveys in Evolutionary Biology.  Oxford University Press, New York, NY.
  * [Jiang et al., 2019] Jiang, L., Zheng, Z., Qi, T., Kemper, K. E., Wray, N. R., Visscher, P. M., and Yang, J. (2019).  A resource-efficient tool for mixed model association analysis of large-scale data.  Nat Genet, 51(12):1749–1755.
  * [Kang et al., 2010] Kang, H. M., Sul, J. H., Service, S. K., Zaitlen, N. A., yee Kong, S., Freimer, N. B., Sabatti, C., and Eskin, E. (2010).  Variance component model to account for sample structure in genome-wide association studies.  Nat Genet, 42(4):348–354.
  * [Kang et al., 2008] Kang, H. M., Zaitlen, N. A., Wade, C. M., Kirby, A., Heckerman, D., Daly, M. J., and Eskin, E. (2008).  Efficient control of population structure in model organism association mapping.  Genetics, 178(3):1709–1723.
  * [Kelleher et al., 2019] Kelleher, J., Wong, Y., Wohns, A. W., Fadil, C., Albers, P. K., and McVean, G. (2019).  Inferring whole-genome histories in large population datasets.  Nat Genet, 51(9):1330–1338.
  * [Kim et al., 2021] Kim, J., Shen, J., Wang, A., Mehrotra, D. V., Ko, S., Zhou, J. J., and Zhou, H. (2021).  VCSEL: Prioritizing SNP-set by penalized variance component selection.  Ann. Appl. Stat., 15(4).
  * [Kimura, 1969] Kimura, M. (1969).  The number of heterozygous nucleotide sites maintained in a finite population due to steady flux of mutations.  Genetics, 61(4):893–903.
  * [Link et al., 2023] Link, V., Schraiber, J. G., Fan, C., Dinh, B., Mancuso, N., Chaing, C. W., and Edge, M. D. (2023).  Tree-based qtl mapping with expected local genetic relatedness matrices.  Unpublished work.
  * [Loh et al., 2015] Loh, P.-R., Tucker, G., Bulik-Sullivan, B. K., Vilhjálmsson, B. J., Finucane, H. K., Salem, R. M., Chasman, D. I., Ridker, P. M., Neale, B. M., Berger, B., Patterson, N., and Price, A. L. (2015).  Efficient bayesian mixed-model analysis increases association power in large cohorts.  Nat Genet, 47(3):284–290.
  * [McVean, 2009] McVean, G. (2009).  A genealogical interpretation of principal components analysis.  PLoS Genetics, 5(10):e1000686.
  * [O'Connor et al., 2019] O'Connor, L. J., Schoech, A. P., Hormozdiari, F., Gazal, S., Patterson, N., and Price, A. L. (2019).  Extreme polygenicity of complex traits is explained by negative selection.  The American Journal of Human Genetics, 105(3):456–476.
  * [Price et al., 2006] Price, A. L., Patterson, N. J., Plenge, R. M., Weinblatt, M. E., Shadick, N. A., and Reich, D. (2006).  Principal components analysis corrects for stratification in genome-wide association studies.  Nat Genet, 38(8):904–909.
  * [Ralph, 2019] Ralph, P. L. (2019).  An empirical approach to demographic inference with genomic data.  Theoretical Population Biology, 127:91–101.
  * [Schaid et al., 2020] Schaid, D. J., Sinnwell, J. P., Larson, N. B., and Chen, J. (2020).  Penalized variance components for association of multiple genes with traits.  Genetic Epidemiology, 44(7):665–675.
  * [Shi et al., 2016] Shi, H., Kichaev, G., and Pasaniuc, B. (2016).  Contrasting the genetic architecture of 30 complex traits from summary association data.  The American Journal of Human Genetics, 99(1):139–153.
  * [Simons et al., 2022] Simons, Y. B., Mostafavi, H., Smith, C. J., Pritchard, J. K., and Sella, G. (2022).  Simple scaling laws control the genetic architectures of human complex traits.  bioRxiv.
  * [Speed and Balding, 2014] Speed, D. and Balding, D. J. (2014).  Relatedness in the post-genomic era: is it still useful?  Nat Rev Genet, 16(1):33–44.
  * [Speidel et al., 2019] Speidel, L., Forest, M., Shi, S., and Myers, S. R. (2019).  A method for genome-wide genealogy estimation for thousands of samples.  Nat Genet, 51(9):1321–1329.
  * [Thompson, 2013] Thompson, E. A. (2013).  Identity by descent: Variation in meiosis, across genomes, and in populations.  Genetics, 194(2):301–326.
  * [Tibshirani, 1996] Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288.
  * [Timpson et al., 2017] Timpson, N. J., Greenwood, C. M. T., Soranzo, N., Lawson, D. J., and Richards, J. B. (2017).  Genetic architecture: the shape of the genetic contribution to human traits and disease.  Nat Rev Genet, 19(2):110–124.
  * [Veller and Coop, 2023] Veller, C. and Coop, G. (2023).  Interpreting population and family-based genome-wide association studies in the presence of confounding.  bioRxiv.
  * [Wang et al., 2021] Wang, S., Ge, S., Colijn, C., Biller, P., Wang, L., and Elliott, L. T. (2021).  Estimating genetic similarity matrices using phylogenies.  Journal of Computational Biology, 28(6):587–600.
  * [Weir and Cockerham, 1984] Weir, B. S. and Cockerham, C. C. (1984).  Estimating $F$ -statistics for the analysis of population structure.  Evolution, 38(6):1358–1370.
  * [Wohns et al., 2022] Wohns, A. W., Wong, Y., Jeffery, B., Akbari, A., Mallick, S., Pinhasi, R., Patterson, N., Reich, D., Kelleher, J., and McVean, G. (2022).  A unified genealogy of modern and ancient genomes.  Science, 375(6583).
  * [Yang et al., 2010] Yang, J., Benyamin, B., McEvoy, B. P., Gordon, S., Henders, A. K., Nyholt, D. R., Madden, P. A., Heath, A. C., Martin, N. G., Montgomery, G. W., Goddard, M. E., and Visscher, P. M. (2010).  Common SNPs explain a large proportion of the heritability for human height.  Nat Genet, 42(7):565–569.
  * [Zeng et al., 2018] Zeng, J., de Vlaming, R., Wu, Y., Robinson, M. R., Lloyd-Jones, L. R., Yengo, L., Yap, C. X., Xue, A., Sidorenko, J., McRae, A. F., Powell, J. E., Montgomery, G. W., Metspalu, A., Esko, T., Gibson, G., Wray, N. R., Visscher, P. M., and Yang, J. (2018).  Signatures of negative selection in the genetic architecture of human complex traits.  Nat Genet, 50(5):746–753.
  * [Zhang et al., 2021] Zhang, B. C., Biddanda, A., and Palamara, P. F. (2021).  Biobank-scale inference of ancestral recombination graphs enables genealogy-based mixed model association of complex traits.  bioRxiv.
  * [Zhang and Pan, 2014] Zhang, Y. and Pan, W. (2014).  Principal component regression and linear mixed model in association analysis of structured samples: Competitors or complements?  Genet. Epidemiol., 39(3):149–155.
  * [Zhou and Stephens, 2012] Zhou, X. and Stephens, M. (2012).  Genome-wide efficient mixed-model analysis for association studies.  Nat Genet, 44(7):821–824.