Different effect sizes of mutations in different demes

Different effect sizes of mutations in different demes#

New in version 0.7.0.

This section describes how to allow mutations to have different effect sizes in different demes.

To allow for correlated effect sizes across demes, we used fwdpy11.mvDES to model multivariate distributions of effect sizes. This class inherits from fwdpy11.Sregion (see Distributions of effect sizes). In general, there is no standard general method for drawing random deviates from arbitrary multivariate distributions. The approach that we take is to use a multivariate Gaussian distribution as the underlying kernel.

First, we will describe the cases where the multivariate Gaussian kernel leads to output effect size distributions that have straightforward interpretations. Then, we will move onto the more general case allowing you to construct multivariate distributions from current fwdpy11.Sregion types.

The multivariate Gaussian#

We may model Gaussian effect sizes using the existing fwdpy11.MultivariateGaussianEffects in conjunction with fwdpy11.mvDES. Using fwdpy11.MultivariateGaussianEffects on its own is used to model pleiotropic effects on a trait. Here, we are re-using this type to model correlated effect sizes across demes.

At this time, it is probably best to look at an example. The following code models Gaussian stabilizing selection on a quantitative trait. The effects sizes within each deme are themselves given by Gaussian distributions and there is no correlation in the effect size in the two demes.

import demes
import fwdpy11
import numpy as np

yaml = """
description: two demes with symmetric migration
time_units: generations
demes:
 - name: deme0
   epochs:
    - start_size: 100
 - name: deme1
   epochs:
     - start_size: 100
migrations:
 - demes: [deme0, deme1]
   rate: 0.1
"""

graph = demes.loads(yaml)
demography = fwdpy11.ForwardDemesGraph.from_demes(graph, 1)

optima = [fwdpy11.Optimum(when=0, optimum=0.0, VS=10.0)]
pdict = {
    "nregions": [],
    "recregions": [],
    "sregions": [
        fwdpy11.mvDES(
            fwdpy11.MultivariateGaussianEffects(0, 1, 1, np.identity(2)), np.zeros(2)
        )
    ],
    "rates": (0, 1e-3, None),
    "demography": demography,
    "simlen": 100,
    "gvalue": fwdpy11.Additive(
        ndemes=2, scaling=2,
        gvalue_to_fitness=fwdpy11.GaussianStabilizingSelection.single_trait(optima)
    ),
    "prune_selected": False,
}

Most of the above is standard. Let’s dissect the new bits:

An instance of fwdpy11.mvDES is our only region with selected mutations.
This instance holds an instance of fwdpy11.MultivariateGaussianEffects that puts mutations on the interval \([0, 1)\) with weight 1 and an identity matrix specifies the correlation in effect sizes between demes 0 and 1. The identity matrix has the value zero for all off-diagonal elements, meaning no covariance in effect sizes across demes.
The final constructor argument specifies the mean of each marginal Gaussian distribution. The means are both zero.
Our genetic value type accepts an ndemes parameter, telling it that it has to look for deme-specific effect sizes. This value must be set to the maximum number of demes that will exist during a simulation.

Let’s evolve the model now:

params = fwdpy11.ModelParams(**pdict)
pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
rng = fwdpy11.GSLrng(1010)
fwdpy11.evolvets(rng, pop, params, 10)

Let’s extract the effect sizes from each deme:

for i in pop.tables.mutations:
    print(pop.mutations[i.key].esizes)

[ 0.09742853 -0.79068629]
[-0.46314495 -0.31686816]

Let’s look at another example where effect sizes covary negatively across demes and raise the mutation rate a bit:

vcv = np.array([1.0, -0.99, -0.99, 1.0]).reshape((2, 2))
pdict["sregions"] = [
    fwdpy11.mvDES(fwdpy11.MultivariateGaussianEffects(0, 1, 1, vcv), np.zeros(2))
]
pdict["rates"] = (0, 5e-3, None)
params = fwdpy11.ModelParams(**pdict)
pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
fwdpy11.evolvets(rng, pop, params, 10)
for i in pop.tables.mutations:
    print(pop.mutations[i.key].esizes)

[-0.65688902  0.67458945]
[ 1.44236014 -1.19184102]
[ 0.53515552 -0.49690106]
[0.11239562 0.03076654]
[ 0.0470034  -0.14227415]
[-0.44707712  0.39814511]
[-0.64860038  0.53073426]
[-0.12680578 -0.00521323]
[-1.702493    1.78230999]
[-0.12798166 -0.14205828]
[ 0.82686083 -0.68915806]

Now we see that the effect sizes often differ in sign between the two demes.

The multivariate lognormal#

If \(X\) is a multivariate Gaussian distribution, \(N(\mathbf{\mu}, \mathbf{\sum})\), where \(\mathbf{\mu}\) is a vector of mean values and \(\mathbf{\sum}\) is the covariance matrix, then \(Y = e^X\) is a multivariate lognormal random variable with mean \(E[Y]_i = e^{\mu_i + \frac{1}{2}\sum_{ii}}\) and covariance matrix \(Var[Y]_{i,j} = e^{\mu_i + \mu_j + \frac{1}{2}(\sum_{ii} + \sum_{jj})}(e^{\sum_{ij}}-1)\).

To specify a multivariate lognormal distribution of effect sizes, we use the static class method fwdpy11.LogNormalS.mv(). The following code constructs a distribution of effect sizes such that -2Ns (where N is the size of a single deme) is a multivariate lognormal with means zero and an identity matrix as a covariance matrix used to specify the multivate Gaussian kernel.

mvdes = fwdpy11.mvDES(
    fwdpy11.LogNormalS.mv(0, 1, 1, scaling=-200), np.zeros(2), np.identity(2)
)

Note

The lognormal distribution returns deviates \(> 0\). To model deleterious mutations/effect sizes < 0, use the scaling parameter with a negative value like we just did!

Let’s put it in a simulation and run it:

pdict = {
    "nregions": [],
    "recregions": [],
    "sregions": [mvdes],
    "rates": (0, 1e-3, None),
    "demography": demography,
    "simlen": 100,
    "gvalue": fwdpy11.Multiplicative(ndemes=2, scaling=2),
}
params = fwdpy11.ModelParams(**pdict)
pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
fwdpy11.evolvets(rng, pop, params, 10)
for i in pop.tables.mutations:
    print(pop.mutations[i.key].esizes)

[-0.00183755 -0.0009758 ]
[-0.00110216 -0.00049425]

“Custom” multivariate distributions#

The previous two sections cover cases where the methods for generating deviates from a multivariate distribution are straightforward and agreed upon.

In order to simulate multivariate distributions of effect sizes based on fwdpy11.Sregion types, we follow a fairly intuitive approach described in [XKS00]. Briefly, the multivariate Gaussian kernel is used to produce deviates. Then, the quantiles from the cummulative distribution of each marginal Gaussian are used to generate a deviate from the desired output distribution of interest.

For a simulation with n populations we need:

A list of n fwdpy11.Sregion objects
An array of n means for the multivariate Gaussian
An n-by-n covariance matrix for the multivariate Gaussian

The following generates exponentially distributed effect sizes in each deme with a high correlation across demes:

mvdes = fwdpy11.mvDES(
    [fwdpy11.ExpS(0, 1, 1, -0.5)] * 2,
    np.zeros(2),
    np.matrix([1, 0.9, 0.9, 1]).reshape((2, 2)),
)
pdict["sregions"] = [mvdes]
params = fwdpy11.ModelParams(**pdict)
pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
fwdpy11.evolvets(rng, pop, params, 10)
for i in pop.tables.mutations:
    print(pop.mutations[i.key].esizes)

[-0.05297616 -0.05994487]
[-0.00600096 -0.00063   ]

We can mix and match our distributions. Here, the distribution of effect sizes in deme 0 is exponential and the distribution in deme 1 is gamma. The two distributions have means with opposite signs and the magnitudes of the marginal deviates negatively covary:

mvdes = fwdpy11.mvDES(
    [fwdpy11.ExpS(0, 1, 1, -0.5), fwdpy11.GammaS(0, 1, 1, mean=0.1, shape_parameter=1)],
    np.zeros(2),
    np.matrix([1, -0.9, -0.9, 1]).reshape((2, 2)),
)
pdict["sregions"] = [mvdes]
params = fwdpy11.ModelParams(**pdict)
pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
fwdpy11.evolvets(rng, pop, params, 10)
for i in pop.tables.mutations:
    print(pop.mutations[i.key].esizes)

[-0.17037555  0.17324704]
[-0.22122016  0.22419496]
[-1.78012379  0.01713378]
[-0.36227095  0.10331541]
[-0.24379116  0.12145104]

The type fwdpy11.ConstantS has intuitive behavior:

mvdes = fwdpy11.mvDES(
    [fwdpy11.ExpS(0, 1, 1, -0.5), fwdpy11.ConstantS(0, 1, 1, -0.1)],
    np.zeros(2),
    np.matrix([1, -0.9, -0.9, 1]).reshape((2, 2)),
)
pdict["rates"] = (0, 5e-3, None)
pdict["sregions"] = [mvdes]
params = fwdpy11.ModelParams(**pdict)
pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
rng = fwdpy11.GSLrng(1010)
fwdpy11.evolvets(rng, pop, params, 10)
for i in pop.tables.mutations:
    print(pop.mutations[i.key].esizes)

[-0.99150137 -0.1       ]
[-0.22828726 -0.1       ]
[-0.10570464 -0.1       ]
[-0.06271967 -0.1       ]

Recipes#

Different signs in different demes#

Consider two demes. You want any beneficial mutation in one deme to be deleterious in the other and vice-versa.

For the multivariate Gaussian, use the covariance matrix as done above. Note that this approach only generates a tendency to different signs in different demes.

With the multivariate lognormal, the best we can do is to use negative correlations such that deleterious mutations in deme 0 are less deleterious in deme 1, etc.:

sregions = [
    fwdpy11.mvDES(
        fwdpy11.LogNormalS.mv(0, 1, 1, scaling=-200),
        np.zeros(2),
        np.matrix([1, -0.99, -0.99, 1]).reshape((2, 2)),
    )
]
sregions.append(
    fwdpy11.mvDES(
        fwdpy11.LogNormalS.mv(0, 1, 1, scaling=200),
        np.zeros(2),
        np.matrix([1, -0.99, -0.99, 1]).reshape((2, 2)),
    )
)
pdict["sregions"] = sregions
params = fwdpy11.ModelParams(**pdict)
pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
rng = fwdpy11.GSLrng(1010)
fwdpy11.evolvets(rng, pop, params, 10)
for i in pop.tables.mutations:
    print(pop.mutations[i.key].esizes)

[0.00039308 0.0693669 ]
[-0.01488498 -0.00173671]
[-0.00514551 -0.00519719]
[-0.00377736 -0.00592563]
[0.00069116 0.0453049 ]
[0.00355528 0.00643812]
[-0.00208255 -0.01409883]
[-0.01590584 -0.00175813]
[-0.00045263 -0.04933448]
[0.00227164 0.01008629]
[0.01128436 0.00198865]
[0.01109052 0.00222656]
[-0.00256051 -0.01002711]
[-0.00152782 -0.01807099]

In the output, we see that an effect size in deme i has a corresponding effect size in deme j that is a about an order of magnitude smaller in absolute value.

For the general approach, simply create a list of objects with the desired mean (or constant) effect sizes. For example:

sregions = [
    fwdpy11.mvDES(
        [fwdpy11.ExpS(0, 1, 1, -0.5), fwdpy11.ExpS(0, 1, 1, 0.1)],
        np.zeros(2),
        np.identity(2),
    )
]
sregions.append(
    fwdpy11.mvDES(
        [fwdpy11.ExpS(0, 1, 1, 0.1), fwdpy11.ExpS(0, 1, 1, -0.5)],
        np.zeros(2),
        np.identity(2),
    )
)
pdict["sregions"] = sregions
params = fwdpy11.ModelParams(**pdict)
pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
rng = fwdpy11.GSLrng(1010)
fwdpy11.evolvets(rng, pop, params, 10)
for i in pop.tables.mutations:
    print(pop.mutations[i.key].esizes)

[-0.99150137  0.08290456]
[ 0.04565745 -0.16002424]
[-0.10570464  0.21656547]
[ 0.41101789 -0.43686025]
[-0.88406416  0.16572598]
[ 0.25571088 -0.08508201]
[ 0.15472792 -0.292425  ]
[-0.14495801  0.08997944]
[-0.06271967  0.15351897]

Polygenic traits, multiple demes, correlated effect sizes, and different optima#

See Adaptation to different optima with multivariate distribution of effect sizes.