`library(MethEvolSIM)`

MethEvolSIM is an R package for simulating DNA methylation dynamics on different genomic structures along a given tree. It allows simulating methylation and demethylation events at loci within two different types of genomic structures: regions with low methylation or with high methylation. The data can be interpreted as:

1- CpG sites inside or outside CpG islands or

2- any type of methylation site in regions of higher or lower methylation or

3- CpG islands inside or outside genomic archipelagos (regions with many CpG islands).

In this document, unless specified otherwise, we use the first interpretation.

As a simple example with default parameter values we show a simulation of evolution of a simplified genomic structure on a tree:

```
# Define the tree:
<- "((a:1, b:1):2, c:2);"
tree # Define the genomic structure. Here 1 island ("U") containing 100 CpGs
# surrounded by 2 non-island structures ("M") each containing 10 CpGs
<- data.frame(n = c(10, 100, 10),
infoStr globalState = c("M", "U", "M"))
# Simulate 1 replicate of data at the tree tips
<- simulate_evolData(infoStr = infoStr, tree = tree)
output #> Using default parameter values
#> Simulating data at root and letting it evolve along given tree: ((a:1, b:1):2, c:2);
```

In the simulation output 0 stands for unmethylated, 1 for methylated and 0.5 for partially methylated or hemimethylated.

```
# In replicate 1, tip 2: name of the tip of the tree
$data[[1]][[2]]$name
output#> [1] "b"
# In replicate 1, tip 2: sequence of methylation states of left non-island structure
$data[[1]][[2]]$seq[[1]]
output#> [1] 1 1 1 1 1 1 1 1 1 0
# In replicate 1, tip 2: sequence of methylation states of island structure
$data[[1]][[2]]$seq[[2]]
output#> [1] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.0 0.5 0.0 0.0 0.0
#> [19] 0.0 0.5 0.5 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0
#> [37] 0.0 0.5 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5
#> [55] 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.5 0.5 0.5 0.0 0.0 0.5 0.0 0.0 0.5 0.5 0.5
#> [73] 0.5 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5
#> [91] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
# In replicate 1, tip 2: sequence of methylation states of right non-island structure
$data[[1]][[2]]$seq[[3]]
output#> [1] 1 1 1 1 1 1 1 1 1 0
```

Structures are categorized as CpG islands or non-islands (alternatively regions with low/high methylation or archipelagos/non-archipelagos). In the root of the genealogy we assume that each structure samples an equilibrium probability triple \((\pi_u,\pi_p, \pi_m)\) for the methylation states: unmethylated, partially-methylated and methylated, respectively. Current parameter values favor prevalent homogeneity in methylation levels within structures, most CpG islands predominantly unmethylated and a small proportion predominantly methylated, and most non-islands predominantly methylated and a small proportion predominantly unmethylated. The sampling of the equilibrium probability triple is done using two distinct Beta distributions to sample the corresponding methylation frequencies.

**Island parameterization: alpha_pI, beta_pI, alpha_mI,
beta_mI**

MethEvolSIM samples \(\pi_p\) from the first Beta distribution. Current parameterization uses as initial parameter values a relatively small \(\alpha_{I_p}\) and relatively large \(\beta_{I_p}\) to encourage minor assignment of partially methylated states.

The second Beta distribution, scaled by \(1 - \pi_p\), corresponds to the proportion of methylated and unmethylated states. Currently both parameter values (\(\alpha_{I_m}\) and \(\beta_{I_m}\)) set to small values while \({\alpha_{I_m} < \beta_{I_m}}\) to favour shifting the stochastic choice toward homogeneous states with a higher proportion of predominantly unmethylated islands, and assigned to \(\pi_m\). Finally, we set \(\pi_u\) as \(1 - \pi_p - \pi_m\).

**Non-island parameterization: alpha_pNI, beta_pNI, alpha_mNI,
beta_mNI**

Following the same rationale as for island parameterization, each non-island structure samples a common equilibrium probability triple from two Beta distributions with parameters \(\alpha_{NI_p}\), \(\beta_{NI_p}\), \(\alpha_{NI_m}\) and \(\beta_{NI_m}\). The parameter values are currently set to favor the characteristic high overall methylation level.

**Parameter ranges** The alpha and beta parameters of a
Beta distribution must be greater than 0, that is in (0, \(\infty\)).

DNA methylation state evolution happens by two structure-specific events: CpG single site events (SSEs) and CpG island wide events (IWEs).

IWEs change the island methylation probabilities and some of the CpG sites in the CpG island can simultaneously change their state at the same time. Different CpG islands, however, are assumed to evolve independently of each other. IWEs occur per CpG island at a rate denoted by \(\mu\). Consequently, over a branch of length \(l\), the expected number of IWEs per island is given by \(\mu \cdot l\).

**IWE parameterization: mu. Parameter range:** Should be
a small number, but it can take any value from 0, that is in \([0, \infty)\)

In addition to IWEs we allow for single site events (SSEs), which change the methylation states of single CpG sites within and between CpG islands. In our model SSEs account for two dynamics: independent of each other (SSEi), and assuming correlations between adjacent sites (SSEc).

As SSEc transition probabilities are dependent on adjacent sites we use Gillespie’s \(\tau\)-leap approximation. to allow for CpG sites to o update their neighbouring state and rates of change after short time steps. Therefore, we first discretize the branch intervals between IWEs in time intervals of default length \(dt =0.01\). For each CpG position and each time interval the rate matrix for the transitions between the states unmethylated, partially methylated and methylated results from the addition of the rate matrices of the two types of SSE events, so that \(\mathbb{E}(R_i + R_c)=1\), \(R_i\) and \(R_c\) representing the rate factors of independent and correlated processes. The independent process samples the new methylation state from the structure’s equilibrium frequencies and the correlated process copies the left or right neighbor state with equal probabilities. The proportion of each of the two SSE types is determined by the model parameter \(\iota\).

Different positions \(j\) have different SSEi rates \(R_{i,j}\) coming from a discretized gamma distribution with 3 categories, expected value \(\iota\) and shape parameters \(\alpha_{R_{i}}\) and \({\beta_{R_{i}}}=\frac{\alpha_{R_{i}}}{\iota}\). The probability to be in each respective rate category is 1/3. For the rate of collaborative SSEs \(R_c = 1-\iota\), the probability of considering left of right neighbour is equal.

**SSE process parameterization: iota, alpha_Ri**
**Parameter ranges** The alpha parameter of a Gamma
distribution must be greater than 0, that is (0, \(\infty\)). Iota must be between 0
(non-included) and 1, that is in (0,1].

The initial spatial distribution of CpG sites in islands and non-island structures is given by the user through a data frame with 2 columns named “n” and “globalState” , which encode the structural information. Each row provides the information of a different genomic structure.

**n**number of positions for each structure (e.g. number of CpG sites).**globalState**encodes island and non-island structures with the character value “U” for unmethylated regions (e.g. CpG islands) and “M” for methylated regions (e.g. non-CpG islands).

```
# Example 3 structures of length 100 each:
# one non-island containing 100 CpG sites,
# one island containing 100 CpG sites,
# and one non-island containing 100 CpG sites
<- data.frame(n = c(100, 100, 100),
infoStr globalState = c("M", "U", "M"))
infoStr#> n globalState
#> 1 100 M
#> 2 100 U
#> 3 100 M
```

```
# Example for 2 islands: first of length 1, second of length 15
<- data.frame(n = c(1, 15),
infoStr globalState = c("U", "U"))
infoStr#> n globalState
#> 1 1 U
#> 2 15 U
```

The initial methylation state distribution assigns a triple of methylation equilibrium frequencies() \(\pi_u\) for unmethylated, \(\pi_p\) for partially-methylated, and \(\pi_m\) for methylated) to each of the simulated genomic structures. This assignation can be done in two ways:

1 - Sampled as specified in the section “Model parameters: Initial
methylation state distribution” when the given `infoStr`

only
contains `n`

and `globalState`

columns.

2 - Determined by the user providing in `infoStr`

additional columns: `u_eqFreq`

for \(\pi_u\), `p_eqFreq`

for \(\pi_p\) and `m_eqFreq`

for \(\pi_m\). Note that when given, the triple
of frequency values for each of the structures needs to sum up to 1, so
that for each row \(i\), \(\pi_u(i)+\pi_p(i)+\pi_m(i)=1\).

```
# Example 3 structures of length 100 with customized initial methylation state
# equilibrium frequencies
<- data.frame(n = c(100, 100, 100),
infoStr globalState = c("M", "U", "M"),
u_eqFreq = c(0.1, 0.8, 0.1),
p_eqFreq = c(0.1, 0.1, 0.1),
m_eqFreq = c(0.8, 0.1, 0.8))
infoStr#> n globalState u_eqFreq p_eqFreq m_eqFreq
#> 1 100 M 0.1 0.1 0.8
#> 2 100 U 0.8 0.1 0.1
#> 3 100 M 0.1 0.1 0.8
```

The default parameter values used by MethEvolSIM can be obtained
using the function `get_parameterValues()`

, which returns a
data frame with the columns named as specified in the section “Model
parameters”:

```
<- get_parameterValues()
default_paramValues
default_paramValues#> alpha_pI beta_pI alpha_mI beta_mI alpha_pNI beta_pNI alpha_mNI beta_mNI mu
#> 1 0.1 1 0.1 0.5 0.1 1 0.5 0.1 0.1
#> alpha_Ri iota Ri1 Ri2 Ri3
#> 1 0.1 0.3 2.80506e-06 0.005770752 0.8942264
```

Additionally, it returns the three rates of change for SSEi
process,`Ri1`

, `Ri2`

, and `Ri3`

(see
section “Model Parameters: Single-site events (SSEs)”).

The section “Model parameters” provides further information about
each of the parameters, including the range of possible values for each
parameter. To simulate data with a different parameterization, the new
value of the chosen parameter(s) can be modified in the data frame, but
for the cases of `Ri1`

, `Ri2`

, and
`Ri3`

, which are modified through `iota`

and
`alpha_Ri`

.

```
# Example: modification of parameter iota to value 0.2
$iota <- 0.2
default_paramValues
default_paramValues#> alpha_pI beta_pI alpha_mI beta_mI alpha_pNI beta_pNI alpha_mNI beta_mNI mu
#> 1 0.1 1 0.1 0.5 0.1 1 0.5 0.1 0.1
#> alpha_Ri iota Ri1 Ri2 Ri3
#> 1 0.1 0.2 2.80506e-06 0.005770752 0.8942264
```

Additional details regarding the `get_parameterValues()`

function can be found in its documentation.

DNA Methylation Data is simulated using three R6 classes:
`singleStructureGenerator`

,
`combiStructureGenerator`

, and
`treeMultiRegionSimulator`

. Instances of the
`singleStructureGenerator`

class represent individual genomic
structures (e.g. islands or non-islands) and are encapsulated within
instances of the `combiStructureGenerator`

class, which, in
turn, represent the simulated data structure (e.g. a sequence consisting
of a non-island containing 100 CpGs, an island containing 100 CpGs, and
a non-island containing 100 CpGs). Finally, instances of the
`treeMultiRegionSimulator`

class contain a
`combiStructureGenerator`

object at each node of a given
tree, including its root and leaves.

Generally, simulations are conducted as explained in subsection
“Simulation of Data Evolution”, with a documented user-interface
function explaining the output structure:
`simulate_evolData()`

. Alternativelly, the user can simulate
initial data so that the root of the tree is fixed as explained in
subsection “Simulation of Initial Data”. The latter option generates
objects of class `combiStructureGenerator`

containing objects
of class `singleStructureGenerator`

. In the subsection
“Simulation of Initial Data: Exploration of Initial Data” an example to
interact with the classes is given, while further information about the
attributes and methods of each class is provided in their
documentation.

The initial data can be generated as an object of class
`combiStructureGenerator`

.

```
# Example with customized initial methylation frequencies and customized parameter values
<- data.frame(n = c(100, 100, 100),
custom_infoStr globalState = c("M", "U", "M"),
u_eqFreq = c(0.1, 0.8, 0.1),
p_eqFreq = c(0.1, 0.1, 0.1),
m_eqFreq = c(0.8, 0.1, 0.8))
<- get_parameterValues()
custom_params $mu <- 0.005
custom_params<- simulate_initialData(infoStr = custom_infoStr, params = custom_params)
initial_customD # Returns customized parameters and simulated data
$params
initial_customD#> alpha_pI beta_pI alpha_mI beta_mI alpha_pNI beta_pNI alpha_mNI beta_mNI mu
#> 1 0.1 1 0.1 0.5 0.1 1 0.5 0.1 0.005
#> alpha_Ri iota Ri1 Ri2 Ri3
#> 1 0.1 0.3 2.80506e-06 0.005770752 0.8942264
$data
initial_customD#> <combiStructureGenerator>
#> Public:
#> add_offspring_index: function (i)
#> branch_evol: function (branch_length, dt, testing = FALSE)
#> clone: function (deep = FALSE)
#> copy: function ()
#> get_island_index: function ()
#> get_island_number: function ()
#> get_IWE_events: function ()
#> get_mu: function ()
#> get_name: function ()
#> get_offspring_index: function ()
#> get_own_index: function ()
#> get_parent_index: function ()
#> get_singleStr: function (i)
#> get_singleStr_number: function ()
#> initialize: function (infoStr, params = NULL, testing = FALSE)
#> set_IWE_events: function (a)
#> set_name: function (a)
#> set_offspring_index: function (i)
#> set_own_index: function (i)
#> set_parent_index: function (i)
#> set_singleStr: function (singStrList)
#> Private:
#> interval_evol: function (interval_length, dt, testing = FALSE)
#> IWE_events: NULL
#> IWE_rate: 0.005
#> mu: 0.005
#> name: NULL
#> offspring_index: NULL
#> own_index: NULL
#> parent_index: NULL
#> set_IWE_rate: function ()
#> singleStr: list
#> singleStr_globalState: M U M
#> SSE_evol: function (dt, testing = FALSE)
```

```
# Example with sampled initial methylation frequencies and default parameter values
<- data.frame(n = c(100, 100, 100),
custom_infoStr globalState = c("M", "U", "M"))
<- simulate_initialData(infoStr = custom_infoStr)
initialD # Returns default parameters
$params
initialD#> alpha_pI beta_pI alpha_mI beta_mI alpha_pNI beta_pNI alpha_mNI beta_mNI mu
#> 1 0.1 1 0.1 0.5 0.1 1 0.5 0.1 0.1
#> alpha_Ri iota Ri1 Ri2 Ri3
#> 1 0.1 0.3 2.80506e-06 0.005770752 0.8942264
$data
initialD#> <combiStructureGenerator>
#> Public:
#> add_offspring_index: function (i)
#> branch_evol: function (branch_length, dt, testing = FALSE)
#> clone: function (deep = FALSE)
#> copy: function ()
#> get_island_index: function ()
#> get_island_number: function ()
#> get_IWE_events: function ()
#> get_mu: function ()
#> get_name: function ()
#> get_offspring_index: function ()
#> get_own_index: function ()
#> get_parent_index: function ()
#> get_singleStr: function (i)
#> get_singleStr_number: function ()
#> initialize: function (infoStr, params = NULL, testing = FALSE)
#> set_IWE_events: function (a)
#> set_name: function (a)
#> set_offspring_index: function (i)
#> set_own_index: function (i)
#> set_parent_index: function (i)
#> set_singleStr: function (singStrList)
#> Private:
#> interval_evol: function (interval_length, dt, testing = FALSE)
#> IWE_events: NULL
#> IWE_rate: 0.1
#> mu: 0.1
#> name: NULL
#> offspring_index: NULL
#> own_index: NULL
#> parent_index: NULL
#> set_IWE_rate: function ()
#> singleStr: list
#> singleStr_globalState: M U M
#> SSE_evol: function (dt, testing = FALSE)
```

To get the parameter values of a `combiStructureGenerator`

object, it can alternativelly be given as argument to the function
`get_parameterValues()`

.

```
<- initialD$data
combiStr_object get_parameterValues(rootData = combiStr_object)
#> alpha_pI beta_pI alpha_mI beta_mI alpha_pNI beta_pNI alpha_mNI beta_mNI mu
#> 1 0.1 1 0.1 0.5 0.1 1 0.5 0.1 0.1
#> alpha_Ri iota Ri1 Ri2 Ri3
#> 1 0.1 0.3 2.80506e-06 0.005770752 0.8942264
```

The methods of the class `combiStructureGenerator`

can be
used to access each of the contained structures:

```
# E.g. access fist structure (non-island)
$get_singleStr(1) combiStr_object
```

Each contained structure is an object of
`singleStructureGenerator`

and its methods can be used to get
the contained information.

```
# E.g. get the methylation equilibrium frequencies of the first singleStructureGenerator
$get_singleStr(1)$get_eqFreqs()
combiStr_object#> [1] 6.214339e-03 1.002724e-10 9.937857e-01
```

```
# E.g. get the sequence of methylation states of the second singleStructureGenerator
# Encoded as: 1 for unmethylated, 2 for partially-methylated, 3 for methylated
$get_singleStr(2)$get_seq()
combiStr_object#> [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
```

The evolution data is generated with the function
`simulate_evolData()`

. It returns the parameters used
(`$params`

), the length of the time step used for the
Gillespie’s \(\tau\)-leap approximation
(`$dt`

, default 0.01), the tree used (`$tree`

),
and the simulated data (`$data`

). The description of the
simulated data is available in the function’s documentation.

It takes as parameters:

`tree`

: A string representing a tree in Newick format.`infoStr`

: A dataframe with the information of the number of CpGs corresponding to each island and non-island region (columns named`n`

and`globalState`

). Additionally the initial methylation frequencies can be specified adding the columns`u_eqFreq`

,`p_eqFreq`

and`m_eqFreq`

.`dt`

: The length of time steps for Gillespie’s \(\tau\)-leap approximation. Numerical, default is 0.01.`n_rep`

: The number of replicates. Integer, default is 1.`only_tip`

: Logical indicating whether to extract data only for tips (TRUE, default) or to extract the information for all the tree branches (FALSE).`params`

: Default NULL for simulation with default parameter values. When customized parameter values are used the output dataframe of`get_parameterValues()`

can be modified and given to the function. E.g:

```
# Example with customized initial methylation frequencies, customized parameter values
# and default dt (0.01)
<- "((a:1, b:1):2, c:2, (d:3.7, (e:4, f:1):3):5);"
tree <- data.frame(n = c(100, 100, 100),
custom_infoStr globalState = c("M", "U", "M"),
u_eqFreq = c(0.1, 0.8, 0.1),
p_eqFreq = c(0.1, 0.1, 0.1),
m_eqFreq = c(0.8, 0.1, 0.8))
<- get_parameterValues()
custom_params $mu <- 0.005
custom_params<- simulate_evolData(infoStr = custom_infoStr, tree = tree, params = custom_params, n_rep = 3, only_tip = TRUE)
evolD #> Simulating data at root and letting it evolve along given tree: ((a:1, b:1):2, c:2, (d:3.7, (e:4, f:1):3):5);
#> Simulating data at root and letting it evolve along given tree: ((a:1, b:1):2, c:2, (d:3.7, (e:4, f:1):3):5);
#> Simulating data at root and letting it evolve along given tree: ((a:1, b:1):2, c:2, (d:3.7, (e:4, f:1):3):5);
```

```
# Returns customized parameters, tree used, time step length for SSE process used (dt)
# and simulated data
$data
evolD$dt
evolD$tree
evolD$params evolD
```

```
# The simulated data is returned as a list. Each element of the list corresponds to a
# simulation replicate:
length(evolD$data)
#> [1] 3
<- evolD$data[[1]]
rep1 <- evolD$data[[2]]
rep2 <- evolD$data[[3]] rep3
```

```
# When only_tip is TRUE, each replicate contains for each tip: its name and the
# sequence of methylation states. E.g.
<- rep1[[1]]
rep1_tip1 $name
rep1_tip1#> [1] "a"
# Seq is a list with the methylation states of each simulated structure.
# In this example, 3 structures: non-island, island and non-island ...
length(rep1_tip1$seq)
#> [1] 3
# Each with 100 CpGs
length(rep1_tip1$seq[[1]])
#> [1] 100
```

When simulations are conducted with `only_tip = FALSE`

the
`$data`

output contains more information.

```
<- "((a:1, b:1):2, c:2, (d:3.7, (e:4, f:1):3):5);"
tree <- data.frame(n = c(100, 100, 100),
custom_infoStr globalState = c("M", "U", "M"),
u_eqFreq = c(0.1, 0.8, 0.1),
p_eqFreq = c(0.1, 0.1, 0.1),
m_eqFreq = c(0.8, 0.1, 0.8))
<- get_parameterValues()
custom_params $mu <- 0.005
custom_params<- simulate_evolData(infoStr = custom_infoStr, tree = tree, params = custom_params, n_rep = 3, only_tip = FALSE)
evolD #> Simulating data at root and letting it evolve along given tree: ((a:1, b:1):2, c:2, (d:3.7, (e:4, f:1):3):5);
#> Simulating data at root and letting it evolve along given tree: ((a:1, b:1):2, c:2, (d:3.7, (e:4, f:1):3):5);
#> Simulating data at root and letting it evolve along given tree: ((a:1, b:1):2, c:2, (d:3.7, (e:4, f:1):3):5);
```

```
# When only_tip is FALSE. $data output contains the simulated data and the
# information on the relationship between tree branches.
names(evolD$data)
#> [1] "sim_data" "branchInTree"
```

```
# The information of the tree branches can be accessed using $branchInTree.
<- evolD$data$branchInTree
treeStr # It is a list in which each element index represents the index of the branch
<- treeStr[[1]] # First branch (or tree root)
root <- treeStr[[2]] # Second branch
b2 # ...
# Each branch contains information of the offspring and parent indexes.
$parent_index # The tree root does not have parent branches
root#> NULL
$offspring_index # The tree root has 3 daughter branches: 2, 5 and 6.
root#> [1] 2 5 6
$parent_index # Branch 2 therefore has as parent branch the root.
b2#> [1] 1
$offspring_index # Branch 2 also has 2 daughter branches: 3 and 4
b2#> [1] 3 4
```

```
# The simulated data can be accessed using $sim_data.
<- evolD$data$sim_data
sim_data # As before, each list element corresponds to a simulation replicate.
# E.g. replicate 3
<- sim_data[[3]]
rep3 # In each replicate, each element of the list corresponds to a tree branch
# the indexes correspond to the information in $branchInTree
<- rep3[[1]]
root <- rep3[[2]]
b2 <- rep3[[3]] b3
```

```
# Each tree branch contains:
# - the branch name (NULL for tree root and inner nodes and
# tip name for the tree tips):
$name
root#> NULL
$name
b2#> NULL
$name
b3#> [1] "a"
```

```
# - the information on IWE events that happened in that branch.
$IWE # The root always has NULL because its branch length is 0
root#> NULL
# The rest of the branches have FALSE when no IWE happened or a list containing
# $islands corresponds to the island index(es) to which the event(s) happened
# $times corresponds to the branch time in which the event happened.
$IWE
b2#> [1] FALSE
$IWE
b3#> [1] FALSE
```

```
# - A list with each element representing the sequence of methylation states of
# each simulated island and non-island
# Encoded as 0 for unmethylated, 0.5 for partially methylated and 1 for methylated
$seq[[1]] # First sequence (non-island)
root#> [1] 1.0 1.0 0.0 1.0 0.0 0.0 1.0 1.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0
#> [19] 1.0 1.0 1.0 1.0 0.5 1.0 1.0 1.0 1.0 0.5 0.0 1.0 1.0 1.0 0.5 1.0 1.0 1.0
#> [37] 1.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.5 1.0 1.0 1.0 0.5 0.5
#> [55] 1.0 1.0 1.0 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.5 1.0 1.0 1.0 1.0
#> [73] 1.0 0.5 1.0 1.0 1.0 1.0 1.0 1.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 0.5 1.0 1.0
#> [91] 1.0 1.0 0.5 1.0 0.5 1.0 1.0 1.0 1.0 0.0
$seq[[2]] # Second sequence (island)
root#> [1] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> [19] 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0
#> [37] 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.5 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.5 0.0
#> [55] 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> [73] 0.5 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> [91] 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0
```

```
# - A list with each element representing the methylation frequencies for the states: unmethylated, partially-methylated and methylated
$eqFreqs[[1]][1] # First structure of root (non-island) equilibrium frequency for unmethylated state
root#> [1] 0.1
$eqFreqs[[1]][2] # First structure of root (non-island) equilibrium frequency for partially-methylated state
root#> [1] 0.1
$eqFreqs[[1]][3] # First structure of root (non-island) equilibrium frequency for methylated state
root#> [1] 0.8
```

`rootData`

: Alternative argument to`infoStr`

(only one can be given). It enables simulating replicates of evolution along a tree with fixed data at root. When used, it takes the output of`simulate_initialData()$data`

. If customized parameter values are to be used, then the customized parameters are given to`simulate_initialData()`

instead to`simulate_evolData()`

.

```
# Example with customized initial methylation frequencies, customized parameter values and default dt (0.01)
<- "((a:1, b:1):2, c:2, (d:3.7, (e:4, f:1):3):5);"
tree <- data.frame(n = c(100, 100, 100),
custom_infoStr globalState = c("M", "U", "M"),
u_eqFreq = c(0.1, 0.8, 0.1),
p_eqFreq = c(0.1, 0.1, 0.1),
m_eqFreq = c(0.8, 0.1, 0.8))
<- get_parameterValues()
custom_params $mu <- 0.005
custom_params<- simulate_initialData(infoStr = custom_infoStr, params = custom_params)$data
initialD <- simulate_evolData(rootData =initialD, tree = tree)
evolD #> Parameter values set as in given rootData
#> Simulating evolution of given data at root along given tree: ((a:1, b:1):2, c:2, (d:3.7, (e:4, f:1):3):5);
```