Plotting >100 Ma of evolution

[LotharukpongJS]

Dear all,

I am very intrigued by the ideas proposed in this project and the papers motivating this study. I tried simulating the evolution for 1000 Ma (1 billion years) and found some issues with the software implementation.

While the simulation works with 100 million years & 1 million years per timestep, the simulation doesn't work with 1,000 years & 1 million years per timestep as well as with 100 million years & 10 million years per timestep. Incidentally, 10 million years & 0.01 million years per timestep works. Do you know what is the cause behind this? There is no error message in the R console. All results are displayed below.

I initially thought that the issue stemmed from too many steps in the simulation. Instead, it may be due to the (mis-)estimation of simulated values which breaks the simulation. This is because, with 100 million years & 10 million years per timestep, one would expect it to work if the issue was due to too many steps in the simulation. Using the same number of steps 10 million years & 0.01 million years per timestep (i.e. 10/0.01 = 1000 = 1000/1 from 1,000 years & 1 million years per timestep), the simulation works.

100 million years, 1 million years per timestep [works]

> simulateEvolution(simulationParams = SimulationParams(numSteps = 100, timePerStep = 1))
[1] "simulationOutput.out"
> data <- parseSimulationOutput("simulationOutput.out")
> summariseEvolution(data, "pairwise")
# A tibble: 10 × 8
# Groups:   realTime [10]
   realTime  step pairwiseMin pairwiseQ25 pairwiseMedian pairwiseQ75 pairwiseMax
      <dbl> <dbl>       <dbl>       <dbl>          <dbl>       <dbl>       <dbl>
 1       10    10           1       0.628          0.622       0.617       0.600
 2       20    20           1       0.700          0.575       0.562       0.508
 3       30    30           1       0.767          0.649       0.557       0.500
 4       40    40           1       0.858          0.633       0.567       0.5  
 5       50    50           1       0.840          0.699       0.556       0.521
 6       60    60           1       0.938          0.754       0.620       0.502
 7       70    70           1       0.801          0.698       0.598       0.516
 8       80    80           1       0.761          0.631       0.582       0.52 
 9       90    90           1       0.789          0.563       0.534       0.502
10      100   100           1       0.799          0.619       0.539       0.518
# ℹ 1 more variable: pairwiseMean <dbl>
> summariseEvolution(data, "initial")
# A tibble: 10 × 9
# Groups:   realTime [10]
   realTime  step numSequences divergenceMin divergenceQ25 divergenceMedian divergenceQ75
      <dbl> <dbl>        <int>         <dbl>         <dbl>            <dbl>         <dbl>
 1       10    10           50         0.776         0.771            0.766         0.761
 2       20    20           50         0.632         0.623            0.620         0.617
 3       30    30           50         0.547         0.538            0.534         0.530
 4       40    40           50         0.497         0.480            0.478         0.473
 5       50    50           50         0.448         0.441            0.436         0.432
 6       60    60           50         0.421         0.412            0.409         0.406
 7       70    70           50         0.400         0.390            0.388         0.383
 8       80    80           50         0.385         0.375            0.371         0.367
 9       90    90           50         0.374         0.363            0.358         0.352
10      100   100           50         0.359         0.348            0.344         0.342
# ℹ 2 more variables: divergenceMax <dbl>, divergenceMean <dbl>

1,000 million years, 1 million years per timestep [doesn't work]

> simulateEvolution(simulationParams = SimulationParams(numSteps = 1000, timePerStep = 1))
[1] "simulationOutput.out"
> data2 <- parseSimulationOutput("simulationOutput.out")
> summariseEvolution(data2, "pairwise")
# A tibble: 1 × 8
# Groups:   realTime [1]
  realTime  step pairwiseMin pairwiseQ25 pairwiseMedian pairwiseQ75 pairwiseMax
     <dbl> <dbl>       <dbl>       <dbl>          <dbl>       <dbl>       <dbl>
1      100   100           1       0.799          0.619       0.539       0.518
# ℹ 1 more variable: pairwiseMean <dbl>
> summariseEvolution(data2, "initial")
# A tibble: 1 × 9
# Groups:   realTime [1]
  realTime  step numSequences divergenceMin divergenceQ25 divergenceMedian divergenceQ75
     <dbl> <dbl>        <int>         <dbl>         <dbl>            <dbl>         <dbl>
1      100   100           50         0.359         0.348            0.344         0.342
# ℹ 2 more variables: divergenceMax <dbl>, divergenceMean <dbl>

100 million years, 10 million years per timestep [doesn't work]

> simulateEvolution(simulationParams = SimulationParams(numSteps = 100, timePerStep = 10))
[1] "simulationOutput.out"
> data3 <- parseSimulationOutput("simulationOutput.out")
> summariseEvolution(data3, "initial")
# A tibble: 1 × 9
# Groups:   realTime [1]
  realTime  step numSequences divergenceMin divergenceQ25 divergenceMedian divergenceQ75
     <dbl> <dbl>        <int>         <dbl>         <dbl>            <dbl>         <dbl>
1      100    10           42         0.357         0.349            0.344         0.338
# ℹ 2 more variables: divergenceMax <dbl>, divergenceMean <dbl>
> summariseEvolution(data3, "pairwise")
# A tibble: 1 × 8
# Groups:   realTime [1]
  realTime  step pairwiseMin pairwiseQ25 pairwiseMedian pairwiseQ75 pairwiseMax
     <dbl> <dbl>       <dbl>       <dbl>          <dbl>       <dbl>       <dbl>
1      100    10           1           1              1           1           1
# ℹ 1 more variable: pairwiseMean <dbl>

10 million years, 0.01 million years per timestep [works]

> simulateEvolution(simulationParams = SimulationParams(numSteps = 10, timePerStep = 0.01))
[1] "simulationOutput.out"
> data4 <- parseSimulationOutput("simulationOutput.out")
> summariseEvolution(data4, "pairwise")
# A tibble: 10 × 8
# Groups:   realTime [10]
   realTime  step pairwiseMin pairwiseQ25 pairwiseMedian pairwiseQ75 pairwiseMax
      <dbl> <dbl>       <dbl>       <dbl>          <dbl>       <dbl>       <dbl>
 1     0.01     1        1          1.00           0.999       0.999       0.998
 2     0.02     2        1.00       0.999          0.999       0.998       0.997
 3     0.03     3        1          0.999          0.998       0.998       0.997
 4     0.04     4        1          0.999          0.998       0.997       0.996
 5     0.05     5        1          0.998          0.997       0.997       0.995
 6     0.06     6        1          0.997          0.996       0.996       0.994
 7     0.07     7        1          0.997          0.996       0.995       0.993
 8     0.08     8        1.00       0.996          0.995       0.994       0.993
 9     0.09     9        1          0.995          0.994       0.994       0.992
10     0.1     10        1          0.995          0.994       0.993       0.992
# ℹ 1 more variable: pairwiseMean <dbl>
> summariseEvolution(data4, "initial")
# A tibble: 10 × 9
# Groups:   realTime [10]
   realTime  step numSequences divergenceMin divergenceQ25 divergenceMedian divergenceQ75
      <dbl> <dbl>        <int>         <dbl>         <dbl>            <dbl>         <dbl>
 1     0.01     1           25         1             1                1.00          0.999
 2     0.02     2           27         1             1.00             0.999         0.999
 3     0.03     3           31         1.00          0.999            0.999         0.999
 4     0.04     4           35         1.00          0.999            0.999         0.999
 5     0.05     5           39         1.00          0.999            0.999         0.998
 6     0.06     6           41         0.999         0.999            0.998         0.998
 7     0.07     7           47         0.999         0.998            0.998         0.997
 8     0.08     8           48         0.999         0.998            0.997         0.997
 9     0.09     9           50         0.999         0.998            0.997         0.997
10     0.1     10           50         0.999         0.997            0.997         0.996
# ℹ 2 more variables: divergenceMax <dbl>, divergenceMean <dbl>

Best,
Sodai