traitEvolSim: Trait Evolution Simulator

In guenardg/MPSEM: Modelling Phylogenetic Signals using Eigenvector Maps

Trait Evolution Simulator

Description

Functions to simulate the evolution of traits as an Ornstein-Uhlenbeck process, optionally with trait optimal value evolving following a Markov process.

Usage

traitEvolSim(
  name,
  sigma = 1,
  step = 1,
  alpha = 0,
  optima = NULL,
  transition = NULL
)

## S3 method for class 'QTEM'
print(x, ...)

simulateTrait(x, tem, state, value, weighting = invDistWeighting, ...)

Arguments

name	Name of the trait whose evolution is to be created.
sigma	Variance parameter of the trait (default: sigma = 1).
step	Simulation step size (default: step = 1).
alpha	Trait selection rate (default: alpha = 0, meaning purely neutral trait evolution)
optima	Trait optimal value(s) (default: optima = NULL).
transition	A square, n-by-n, transition intensity matrix (default: transition = NULL, which is suitable for single trait optima).
x	A QTEM-class or graph-class object.
...	further arguments to be passed to other functions and methods
tem	A QTEM-class object from function traitEvolSim or a list of such objects (see details).
state	One or more initial trait state(s) (default: 1).
value	One or more initial trait value(s).
weighting	A weighting function that determine the contribution of the parent vertices as a function of the evolutionary distances (passed as its first argument. It may have any number of argument, each with a default value, and must accept arbitrary arguments (...). Default: invDistWeighting.

Details

The trait evolution simulator is based on a random walk process of one of categories: Weiner process and the Ornstein-Uhlenbeck process. The Weiner process, also known as 'Brownian motion' is a random neutral process with unbounded traits values. It is characterized using a single variance parameter called sigma that dictates to what extent the trait value fluctuates randomly in time.

The Ornstein-Uhlenbeck process is a random process with attractors (optima), with values varying about these optima (or single optimum). It is notable that whenever multiple optima are possible, only a single optimum is effective at any given time. The Ornstein-Uhlenbeck process involve two more values: the trait optimum effective at a particular time and the selection rate (alpha). The optima represents the value toward which the trait is attracted, whereas the value of alpha dictates to what extent such an attraction occurs. When alpha = 0, the optimum stops being an attractor and the Ornstein-Uhlenbeck process becomes purely a Weiner process, whereas when alpha is very large, the trait value detracts only slightly from its optimal value.

The trait simulation interface enabled to handle traits evolving neutrally (pure Weiner process), and according to an Ornstein-Uhlenbeck process with one or many optima. When many optima are used, transitions among the values are simulated as a Markov process. For that purpose, the simulator needs to be provided with a transition density matrix.

When multiple selection regimes are to be simulated, argument tem is given a list of quantitative trait evolution models instead of a single model. In that case, function simulateTrait will expect the graph (argument x) to have a edge property called ‘tem’. This property contains the indices of the trait evolution models for each edge of the graph. In practice, to have reasonable chances to resolve multiple evolution regimes, the number of selection regimes has to be small with respect to the number of vertices.

Value

Function traitEvolSim returns a QTEM-class object, which consists of 8 member functions:

$getName

It takes no argument and returns the name of the trait.

$getStep

It takes no argument and returns the size of the simulation time step.

$setStep

It takes the time step as an argument and sets, the simulation time step, including the recalculation of the transition intensity matrix in the case of a trait with multiple optima.

$getOptima

It has no argument and returns the trait optimum values.

$getTransition

It takes no argument and returns the transition density matrix.

$getProb

It takes no argument and returns the transition probability matrix.

$updateState

It takes the index of the actual trait optimum and returns the updated index (allows shifts in trait optimum).

updateValue

It takes the value of the trait and, optionally, the index of the optimum state, at returns the updated trait value.

$dumpConfig

It takes no argument and returns the trait evolution simulator parameters as a list.

$print

It takes no argument and prints the parameters of the trait simulator on the screen (returns NULL invisibly).

Function simulateTrait returns a two-column data.frame:

state

The simulated index of the trait optimum at the vertices of the graph.

value

The simulated trait value at the vertices of the graph.

Author(s)

Guillaume Guénard [aut, cre] (<https://orcid.org/0000-0003-0761-3072>), Pierre Legendre [ctb] (<https://orcid.org/0000-0002-3838-3305>)

Examples

## Setting the RNG (for example consistency):
set.seed(2182955)

## A list to contain the trait simulator objects:
tem <- list()

## The first simulator is non-neutral with three optima:
traitEvolSim(
  name = "Trait 1",
  sigma = 1.5,
  alpha = 0.15,
  optima = c(30,50,80),
  transition = matrix(c(NA,0.1,0.0,0.1,NA,0.1,0.0,0.1,NA), 3L, 3L)
) -> tem[[1]]
tem[[1]]

## The second simulator is neutral:
traitEvolSim(
  name = "Trait 2",
  sigma = 2.5,
  step = 1
) -> tem[[2]]
tem[[2]]

## The third simulator is non-neutral with a single optimum:
traitEvolSim(
  name = "Trait 3",
  alpha = 0.05,
  optima = 15
) -> tem[[3]]
tem[[3]]

## The fourth simulator is also non-neutral with a single optimum (but having
## a different value than the latter):
traitEvolSim(
  name = "Trait 4",
  sigma = 2.0,
  alpha = 0.25,
  optima = -25
) -> tem[[4]]
tem[[4]]

## Each simulator has a set of embedded member functions that can be called
## directly.

## Member $getName() returns the name of the trait as follows:
unlist(lapply(tem, function(x) x$getName()))

## Member $getStep() returns the step sizes for the traits as follows:
unlist(lapply(tem, function(x) x$getStep()))

## Member $setStep() sets the step sizes as follows:
lapply(tem, function(x) x$setStep(0.1))

## and returns the recalculated transition probability matrix of simulators
## having a transition intensity matrix (i.e., for non-neutral simulators
## with multiple trait optima).

## This is the modified step sizes:
unlist(lapply(tem, function(x) x$getStep()))

## Member $getOptima() returns the simulator's optimum (if available) or a
## vector of optima (if multiple optima are available) as follows:
lapply(tem, function(x) x$getOptima())

## Member $getTransition() returns the transition intensity matrix, when
## available (NULL otherwise) as follows:
lapply(tem, function(x) x$getTransition())

## Member $getProb() returns the transition intensity matrix, whenever
## available (NULL otherwise) as follows:
lapply(tem, function(x) x$getProb())

## When multiple optima are available, member $updateState() enables to
## simulate the transition from one optimal trait state (the one given as the
## argument) to another (the one which is returned by the function):
state <- 1
newstate <- tem[[1]]$updateState(state)
newstate

## Member $updateValue() simulates the evolution of the trait from one value
## to another. For a non-neutral with multiple optima, The trait state is
## provided using argument state (default: 1, which is the only applicable
## value for a single optimum).
oldvalue <- 31.5
newvalue <- tem[[1]]$updateValue(oldvalue, state=1)
newvalue

## Member $dumpConfig() returns the configuration list, which can be used

cfg <- lapply(tem, function(x) x$dumpConfig())

## The trait evolution simulators can be re-instantiated as follows:
lapply(
  cfg,
  function(x)
    traitEvolSim(
      name = x$name,
      sigma = x$sigma,
      step = x$step,
      alpha = x$alpha,
      optima = x$optima,
      transition = x$transition
    )
) -> tem_clone

## Clean up:
rm(cfg, tem_clone)

## Simulate trait evolution using the four simulators described previously:

## Set step size to 0.05
lapply(tem, function(x, a) x$setStep(a), a = 0.05)

## Results list:
res <- NULL
trNms <- lapply(tem, function(x) x$getName())
res$state <- matrix(NA, 1001L, length(tem), dimnames = list(NULL, trNms))
res$optim <- matrix(NA, 1001L, length(tem), dimnames = list(NULL, trNms))
res$trait <- matrix(NA, 1001L, length(tem), dimnames = list(NULL, trNms))
rm(trNms)

## res$state contains the trait state at the simulation time
## res$optim contains the trait's optimum value at the simulation time
## res$trait contains the trait value at the simulation time

## Setting the optimal state of the four traits at the beginning of the
## simulation period:
res$state[1L,] <- c(2,NA,1,1)  ## NBL trait #2 evolves neutrally.

## Getting the trait optima at the beginning of the simulation period:
unlist(
  lapply(
    tem,
    function(x)
      x$getOptima()[res$state[1L,x$getName()]]
  )
) -> res$optim[1L,]

## Setting the initial trait values:
res$trait[1L,] <- c(50,0,15,-25)

## The state of the simulation at the beginning of the simulation:
head(res$state)
head(res$optim)
head(res$trait)

## Setting RNG state to obtain 
set.seed(1234567)

## This loop simulates time steps #2 through #1001
for(i in 2L:1001L) {
  
  ## Simulate the evolution of the trait states (if relevant, which it is
  ## only for the first trait):
  unlist(
    lapply(
      tem,
      function(x)
        x$updateState(res$state[i - 1L,x$getName()])
    ) 
  )-> res$state[i,]
  
  ## Obtain the optimal trait value (relevant for all traits but the second,
  ## trait, which evolves neutrally):
  unlist(
    lapply(
      tem,
      function(x)
        x$getOptima()[res$state[i,x$getName()]]
    )
  ) -> res$optim[i,]
  
  ## Simulate the evolution of the trait value:
  unlist(
    lapply(
      tem,
      function(x)
        x$updateValue(
          state = res$state[i,x$getName()],
          value = res$trait[i - 1L,x$getName()]
        )
    )
  ) -> res$trait[i,]
}

## Plot the results:
par(mar=c(4,4,1,1))
plot(NA, xlim=c(0,0.05*1000), ylim=range(res$trait), xlab="Time",
     ylab="Trait value", las=1L)
lines(x=0.05*(0:1000), y=res$trait[,1L], col="black")
if(!is.na(tem[[1L]]$getOptima())[1L])
  lines(x=0.05*(0:1000), y=res$optim[,1L], col="black", lty=3L)
lines(x=0.05*(0:1000), y=res$trait[,2L], col="red")
if(!is.na(tem[[2L]]$getOptima())[1L])
  lines(x=0.05*(0:1000), y=res$optim[,2L], col="red", lty=3L)
lines(x=0.05*(0:1000), y=res$trait[,3L], col="blue")
if(!is.na(tem[[3L]]$getOptima())[1L])
  lines(x=0.05*(0:1000), y=res$optim[,3L], col="blue", lty=3L)
lines(x=0.05*(0:1000), y=res$trait[,4L], col="green")
if(!is.na(tem[[4L]]$getOptima())[1L])
  lines(x=0.05*(0:1000), y=res$optim[,4L], col="green", lty=3L)

## A linear evolutionary sequence with random edge lengths between 2 and 5:
randomGraph(
  NV = 100,
  NC = function(...) 1,
  NP = function(...) 1,
  timestep = function(ts_min, ts_max, ...) runif(1, ts_min, ts_max),
  maxDist = function(...) NULL,
  ts_min = 2,
  ts_max = 5
) -> gr_lin

## Simulate a trait (Ornstein Uhlenbeck):
simulateTrait(
  x = gr_lin,
  tem = tem[[1]],
  state = 2,
  value = 50,
  a = 1
) -> simTrait

## Showing the trait values with the optima (OU process):
x <- cumsum(c(0,attr(gr_lin,"edge")$distance))
plot(x=x, y=simTrait$value, type="l", las=1)
lines(x=x, y=c(30,50,80)[simTrait$state], lty=3)

## Simulate an other trait (Brownian motion):
simulateTrait(
  x = gr_lin,
  tem = tem[[2]],
  value = 10,
  a = 1
) -> simTrait

## Showing the trait values:
plot(x=x, y=simTrait$value, type="l", las=1)

## A distance-based network:
N <- 100
coords <- cbind(x=runif(N,-1,1), y=runif(N,-1,1))
rownames(coords) <- sprintf("N%d",1:N)
dst <- dist(coords)
gr_dst <- dstGraph(d=dst, th=0.35, origin=15)

## Simulate a trait (Brownian motion) on the distance-based network:
simulateTrait(
  x = gr_dst,
  tem = tem[[2]],
  value = 0,
  a = 1
) -> simTrait

## Showing the results of the simulation:
gp <- par(no.readonly = TRUE)
par(mar=c(5,5,1,1))

plot(NA, xlim=c(-1,1), ylim=c(-1,1), asp=1, las=1, xlab="X", ylab="Y")
points(x=coords[,"x"], y=coords[,"y"], pch=21, cex=abs(simTrait$value)/3,
       bg=gray(0.5 - 0.5*sign(simTrait$value)))

par(gp)

guenardg/MPSEM documentation built on April 14, 2025, 3:53 p.m.