3. Calculate the lowest elevation path and barrier height between stable states

Jingmeng Cui

library(simlandr)

An important property of the states in a landscape is their (kinetic) stability, characterized by the barrier height between these states and other adjacent states. simlandr also provides tools to calculate the barrier heights from landscapes.

You can use the general function calculate_barrier() to calculate the barrier for most landscapes. There are also specific calculate_barrier_*() functions available. The output of these functions is a `barrier object.

The barrier objects contain the potential function and the position of both states and the saddle point. For 3D landscapes, the minimum energy path (MEP) is also provided. barriers can also be calculated for landscapes from multiple simulations. In this case, remember to set individual_landscape = TRUE in landscape construction functions.

The local minimums are searched in a square space around a given point. The point with the lowest potential value in the given region is set as the position of the stable state. If all the potential values in the region are equal to Umax (which represents ~Inf), the barrier calculation functions will expand the searching area automatically. Use expand = FALSE to disable this feature.

For landscapes from multiple simulations, the searching regions for their starting and ending points can be different. simlandr provides make_barrier_grid_2d() and make_barrier_grid_3d() functions to help you put these settings into a data frame with the correct format.

The barrier objects also provide a ggplot geom object that can be added to the landscape plots to show the starting (white), end (white), and saddle (red) points, as well as the MEP (white line, only for 3d landscapes). Use autolayer(b) to access those geoms.

Below are examples of different barrier calculations. See the help documents of those functions for further details.

Prepare data sets and landscapes (see vignette("landscape"))

single_test <- sim_fun_test(
  arg1 = list(ele1 = 1),
  arg2 = list(ele2 = 1, ele3 = 0)
)

l_single_2d <- make_2d_static(single_test, x = "out1")
l_single_3d <- make_3d_static(single_test, x = "out1", y = "out2")

batch_test <- new_arg_set()
batch_test <- batch_test %>%
  add_arg_ele("arg2", "ele3", 0.2, 0.5, 0.1)
batch_test_grid <- make_arg_grid(batch_test)
batch_test_result <- batch_simulation(batch_test_grid, sim_fun_test,
  default_list = list(
    arg1 = list(ele1 = 0),
    arg2 = list(ele2 = 0, ele3 = 0)
  ),
  bigmemory = FALSE
)
batch_test2 <- new_arg_set()
batch_test2 <- batch_test2 %>%
  add_arg_ele("arg1", "ele1", 0.2, 0.6, 0.2) %>%
  add_arg_ele("arg2", "ele2", 0.2, 0.6, 0.2)
batch_test_grid2 <- make_arg_grid(batch_test2)
batch_test_result2 <- batch_simulation(batch_test_grid2, sim_fun_test,
  default_list = list(
    arg1 = list(ele1 = 0),
    arg2 = list(ele2 = 0, ele3 = 0)
  ),
  bigmemory = FALSE
)

l_batch_3d_m1 <- make_3d_matrix(batch_test_result, x = "out1", y = "out2", cols = "ele3")
#> Wrangling the data...
#> Making the plots...
l_batch_2d_m2 <- make_2d_matrix(batch_test_result2, x = "out1", rows = "ele1", cols = "ele2", individual_landscape = TRUE)
#> Wrangling the data...
#> Making the plots...
l_batch_3d_m2 <- make_3d_matrix(batch_test_result2, x = "out1", y = "out2", rows = "ele1", cols = "ele2", Umax = 10, individual_landscape = TRUE)
#> Wrangling the data...
#> Making the plots...

Frequently used parameters for the family of barrier functions:

start_location_value,end_location_value: the initial position (in value) for searching the start/end point; start_r,end_r: the searching (L1) radius for searching the start/end point.

Barrier calculation for 2d single landscape

b_single_2d <- calculate_barrier(l_single_2d, start_location_value = -2, end_location_value = 2, start_r = 1, end_r = 1)

b_single_2d$local_min_start
#> $U
#> [1] 1.200247
#> 
#> $location
#>   x_index   x_value 
#> 22.000000 -2.677756
b_single_2d$local_min_end
#> $U
#> [1] 1.196954
#> 
#> $location
#>    x_index    x_value 
#> 179.000000   2.677764
b_single_2d$saddle_point
#> $U
#> [1] 2.186844
#> 
#> $location
#>      x_index      x_value 
#> 103.00000000   0.08528284

get_barrier_height(b_single_2d)
#> Warning: `summary()` was deprecated in simlandr 0.3.0.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
#> delta_U_start   delta_U_end 
#>     0.9865968     0.9898904

plot(l_single_2d) + autolayer(b_single_2d)
#> Warning in ggplot2::geom_point(mapping = ggplot2::aes(x = d$x[s_location_x_index], : All aesthetics have length 1, but the data has 200 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#>   a single row.
#> Warning in ggplot2::geom_point(mapping = ggplot2::aes(x = local_min_start$location["x_value"], : All aesthetics have length 1, but the data has 200 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#>   a single row.
#> Warning in ggplot2::geom_point(mapping = ggplot2::aes(x = local_min_end$location["x_value"], : All aesthetics have length 1, but the data has 200 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#>   a single row.

Barrier calculation for 3d single landscape

b_single_3d <- calculate_barrier(l_single_3d, start_location_value = c(-2.5, -2), end_location_value = c(2.5, 0), start_r = 0.3, end_r = 0.3)
plot(l_single_3d, 2) + autolayer(b_single_3d)
#> Warning in ggplot2::geom_point(ggplot2::aes(x = local_min_start$location["x_value"], : All aesthetics have length 1, but the data has 40000 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#>   a single row.
#> Warning in ggplot2::geom_point(ggplot2::aes(x = local_min_end$location["x_value"], : All aesthetics have length 1, but the data has 40000 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#>   a single row.
#> Warning in ggplot2::geom_point(ggplot2::aes(x = saddle_point$location["x_value"], : All aesthetics have length 1, but the data has 40000 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#>   a single row.

Barrier calculation for 2d batch landscape

b_batch_2d_m2 <- calculate_barrier(l_batch_2d_m2, start_location_value = -1, end_location_value = 1, start_r = 0.99, end_r = 0.99)
plot(l_batch_2d_m2) + autolayer(b_batch_2d_m2)

Barrier calculation for 3d batch landscape

b_batch_3d_m2 <- calculate_barrier(l_batch_3d_m2, start_location_value = c(-1, -1), end_location_value = c(1, 1), start_r = 0.9, end_r = 0.9)
plot(l_batch_3d_m2) + autolayer(b_batch_3d_m2)

Set the starting and end values for each landscape

b_batch_3d_m1 <- calculate_barrier(l_batch_3d_m1, start_location_value = c(0, 0), end_location_value = c(2, 1), start_r = 0.3, end_r = 0.6)
#> The U in this range is too high. Searching range expanded...
#> r = c(0.832073521958355,0.764100483821545)
plot(l_batch_3d_m1) + autolayer(b_batch_3d_m1)


## This barrier calculation doesn't find proper local minimums for several landscapes. Specify the searching parameters per landscape manually.
## First, print a template of the data format.

make_barrier_grid_3d(batch_test_grid, start_location_value = c(0, 0), end_location_value = c(2, 1), start_r = 0.3, end_r = 0.6, print_template = TRUE)
#> structure(list(start_location_value = list(c(0, 0), c(0, 0), 
#>     c(0, 0), c(0, 0)), start_r = list(c(0.3, 0.3), c(0.3, 0.3
#> ), c(0.3, 0.3), c(0.3, 0.3)), end_location_value = list(c(2, 
#> 1), c(2, 1), c(2, 1), c(2, 1)), end_r = list(c(0.6, 0.6), c(0.6, 
#> 0.6), c(0.6, 0.6), c(0.6, 0.6))), row.names = c(NA, -4L), class = c("arg_grid", 
#> "data.frame"))
#>   ele_list ele3 start_location_value  start_r end_location_value    end_r
#> 1      0.2  0.2                 0, 0 0.3, 0.3               2, 1 0.6, 0.6
#> 2      0.3  0.3                 0, 0 0.3, 0.3               2, 1 0.6, 0.6
#> 3      0.4  0.4                 0, 0 0.3, 0.3               2, 1 0.6, 0.6
#> 4      0.5  0.5                 0, 0 0.3, 0.3               2, 1 0.6, 0.6

## Then, modify the parameters as you want, and send this `barrier_grid` to the barrier calculation function.
b_batch_3d_m1 <- calculate_barrier(
  l_batch_3d_m1,
  make_barrier_grid_3d(batch_test_grid,
    df =
      structure(list(start_location_value = list(
        c(0, 0), c(0, 0), c(0, 0), c(0, 0)
      ), start_r = list(c(0.2, 0.2), c(0.3, 0.3), c(0.3, 0.3), c(0.3, 0.3)), end_location_value = list(c(1, 0.5), c(1.8, 0.8), c(2, 1), c(2, 1)), end_r = c(
        0.6, 0.6, 0.6, 0.6
      )), row.names = c(NA, -4L), class = c(
        "arg_grid",
        "data.frame"
      ))
  )
)
plot(l_batch_3d_m1) + autolayer(b_batch_3d_m1)


## Now it works well.