diff --git a/DESCRIPTION b/DESCRIPTION
index 5be3bff..34beac9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -28,6 +28,6 @@ Suggests:
     knitr,
     testthat
 VignetteBuilder: knitr
-RoxygenNote: 6.1.1
+RoxygenNote: 7.1.1
 LazyData: true
 Encoding: UTF-8
diff --git a/R/mod.acts.R b/R/mod.acts.R
index 4bdb142..7dd462a 100644
--- a/R/mod.acts.R
+++ b/R/mod.acts.R
@@ -23,6 +23,8 @@ acts_msm <- function(dat, at) {
   status <- dat$attr$status
   diag.status <- dat$attr$diag.status
   race <- dat$attr$race
+  geog.lvl <- dat$param$netstats$geog.lvl
+  race.flag <- dat$param$netstats$race
   age <- dat$attr$age
   stage <- dat$attr$stage
   vl <- dat$attr$vl
@@ -33,7 +35,8 @@ acts_msm <- function(dat, at) {
   # Parameters
   acts.mod <- dat$param$epistats$acts.mod
   acts.aids.vl <- dat$param$acts.aids.vl
-  acts.scale <- dat$param$acts.scale
+  acts.scale.main <- dat$param$acts.scale.main
+  acts.scale.casl <- dat$param$acts.scale.casl
 
   # Construct edgelist
   el <- rbind(dat$el[[1]], dat$el[[2]], dat$el[[3]])
@@ -72,14 +75,43 @@ acts_msm <- function(dat, at) {
   hiv.concord.pos[cp] <- 1
 
   # Model predictions
-  x <- data.frame(ptype = el.mc[, "ptype"],
-                  duration = durations,
-                  race.combo = race.combo,
-                  comb.age = comb.age,
-                  hiv.concord.pos = hiv.concord.pos,
-                  city = 1)
-  rates <- unname(predict(acts.mod, newdata = x, type = "response"))/52
-  rates <- rates * acts.scale
+  if (!is.null(geog.lvl)) {
+    if (race.flag == TRUE) {
+      x <- data.frame(ptype = el.mc[, "ptype"],
+                      duration = durations,
+                      race.combo = race.combo,
+                      comb.age = comb.age,
+                      hiv.concord.pos = hiv.concord.pos,
+                      geogYN = 1)
+      rates <- unname(predict(acts.mod, newdata = x, type = "response"))/52
+    } else {
+      x <- data.frame(ptype = el.mc[, "ptype"],
+                      duration = durations,
+                      comb.age = comb.age,
+                      hiv.concord.pos = hiv.concord.pos,
+                      geogYN = 1)
+      rates <- unname(predict(acts.mod, newdata = x, type = "response"))/52
+    }
+  } else {
+    if (race.flag == TRUE) {
+      x <- data.frame(ptype = el.mc[, "ptype"],
+                      duration = durations,
+                      race.combo = race.combo,
+                      comb.age = comb.age,
+                      hiv.concord.pos = hiv.concord.pos)
+      rates <- unname(predict(acts.mod, newdata = x, type = "response"))/52
+    } else {
+      x <- data.frame(ptype = el.mc[, "ptype"],
+                      duration = durations,
+                      comb.age = comb.age,
+                      hiv.concord.pos = hiv.concord.pos)
+      rates <- unname(predict(acts.mod, newdata = x, type = "response"))/52
+    }
+  }
+
+  rates[x$ptype == 1] <- rates[x$ptype == 1] * acts.scale.main
+  rates[x$ptype == 2] <- rates[x$ptype == 2] * acts.scale.casl
+
   ai <- rpois(length(rates), rates)
   el.mc <- cbind(el.mc, durations, ai)
 
diff --git a/R/mod.condoms.R b/R/mod.condoms.R
index 70353f7..11d8afa 100644
--- a/R/mod.condoms.R
+++ b/R/mod.condoms.R
@@ -24,6 +24,8 @@ condoms_msm <- function(dat, at) {
 
   # Attributes
   race <- dat$attr$race
+  race.flag <- dat$param$netstats$race
+  geog.lvl <- dat$param$netstats$race
   age <- dat$attr$age
   diag.status <- dat$attr$diag.status
   prepStat <- dat$attr$prepStat
@@ -55,32 +57,95 @@ condoms_msm <- function(dat, at) {
 
   ## Main/casual partnerships ##
   mc.parts <- which(el[, "ptype"] != 3)
-  el.mc <- el[mc.parts, ]
-
-  x <- data.frame(ptype = el.mc[, "ptype"],
-                  duration = el.mc[, "durations"],
-                  race.combo = race.combo[mc.parts],
-                  comb.age = comb.age[mc.parts],
-                  hiv.concord.pos = hiv.concord.pos[mc.parts],
-                  prep = any.prep[mc.parts],
-                  city = 1)
-  cond.prob <- unname(predict(cond.mc.mod, newdata = x, type = "response"))
-  el.mc <- cbind(el.mc, cond.prob)
+  el.mc <- el[mc.parts, , drop = FALSE]
+
+  if (nrow(el.mc) > 0) {
+    if (!is.null(geog.lvl)) {
+      if (race.flag == TRUE){
+        x <- data.frame(ptype = el.mc[, "ptype"],
+                        duration = el.mc[, "durations"],
+                        race.combo = race.combo[mc.parts],
+                        comb.age = comb.age[mc.parts],
+                        hiv.concord.pos = hiv.concord.pos[mc.parts],
+                        prep = any.prep[mc.parts],
+                        geogYN = 1)
+        cond.prob <- unname(predict(cond.mc.mod, newdata = x, type = "response"))
+      } else {
+        x <- data.frame(ptype = el.mc[, "ptype"],
+                        duration = el.mc[, "durations"],
+                        comb.age = comb.age[mc.parts],
+                        hiv.concord.pos = hiv.concord.pos[mc.parts],
+                        prep = any.prep[mc.parts],
+                        geogYN = 1)
+        cond.prob <- unname(predict(cond.mc.mod, newdata = x, type = "response"))
+      }
+    } else {
+      if (race.flag == TRUE){
+        x <- data.frame(ptype = el.mc[, "ptype"],
+                        duration = el.mc[, "durations"],
+                        race.combo = race.combo[mc.parts],
+                        comb.age = comb.age[mc.parts],
+                        hiv.concord.pos = hiv.concord.pos[mc.parts],
+                        prep = any.prep[mc.parts])
+        cond.prob <- unname(predict(cond.mc.mod, newdata = x, type = "response"))
+      } else {
+        x <- data.frame(ptype = el.mc[, "ptype"],
+                        duration = el.mc[, "durations"],
+                        comb.age = comb.age[mc.parts],
+                        hiv.concord.pos = hiv.concord.pos[mc.parts],
+                        prep = any.prep[mc.parts])
+        cond.prob <- unname(predict(cond.mc.mod, newdata = x, type = "response"))
+      }
+    }
+    el.mc <- cbind(el.mc, cond.prob)
+  }
 
   ## One-off partnerships ##
   oo.parts <- which(el[, "ptype"] == 3)
-  el.oo <- el[oo.parts, ]
-
-  x <- data.frame(race.combo = race.combo[oo.parts],
-                  comb.age = comb.age[oo.parts],
-                  hiv.concord.pos = hiv.concord.pos[oo.parts],
-                  prep = any.prep[oo.parts],
-                  city = 1)
-  cond.prob <- unname(predict(cond.oo.mod, newdata = x, type = "response"))
-  el.oo <- cbind(el.oo, cond.prob)
-
-  ## Bind el together
-  el <- rbind(el.mc, el.oo)
+  el.oo <- el[oo.parts, , drop = FALSE]
+
+  if (nrow(el.oo) > 0) {
+    if (!is.null(geog.lvl)) {
+      if (race.flag == TRUE){
+        x <- data.frame(race.combo = race.combo[oo.parts],
+                        comb.age = comb.age[oo.parts],
+                        hiv.concord.pos = hiv.concord.pos[oo.parts],
+                        prep = any.prep[oo.parts],
+                        geogYN = 1)
+        cond.prob <- unname(predict(cond.oo.mod, newdata = x, type = "response"))
+        el.oo <- cbind(el.oo, cond.prob)
+      } else {
+        x <- data.frame(comb.age = comb.age[oo.parts],
+                        hiv.concord.pos = hiv.concord.pos[oo.parts],
+                        prep = any.prep[oo.parts],
+                        geogYN = 1)
+        cond.prob <- unname(predict(cond.oo.mod, newdata = x, type = "response"))
+        el.oo <- cbind(el.oo, cond.prob)
+      }
+    } else {
+      if (race.flag == TRUE){
+        x <- data.frame(race.combo = race.combo[oo.parts],
+                        comb.age = comb.age[oo.parts],
+                        hiv.concord.pos = hiv.concord.pos[oo.parts],
+                        prep = any.prep[oo.parts],
+                        geogYN = 1)
+        cond.prob <- unname(predict(cond.oo.mod, newdata = x, type = "response"))
+        el.oo <- cbind(el.oo, cond.prob)
+      } else {
+        x <- data.frame(comb.age = comb.age[oo.parts],
+                        hiv.concord.pos = hiv.concord.pos[oo.parts],
+                        prep = any.prep[oo.parts],
+                        geogYN = 1)
+        cond.prob <- unname(predict(cond.oo.mod, newdata = x, type = "response"))
+        el.oo <- cbind(el.oo, cond.prob)
+      }
+    }
+
+    ## Bind el together
+    el <- rbind(el.mc, el.oo)
+  } else {
+    el <- el.mc
+  }
 
   # Acts
   ai.vec <- el[, "ai"]
diff --git a/R/mod.dat_updater.R b/R/mod.dat_updater.R
new file mode 100644
index 0000000..4f16b15
--- /dev/null
+++ b/R/mod.dat_updater.R
@@ -0,0 +1,86 @@
+#' Update list `x` using the elements of list `new_x`
+#'
+#' @param x a list
+#' @param new_x a list
+#'
+#' @return the full `x` list with the modifications added by `new_x`
+#'
+#' @details
+#' This function updates list `x` by name. If `x` and `new_x` elements are not
+#' named, the function will not work properly.
+#' If a function is provided to replace an element that was originaly not a
+#' function, this function will be applied to the original value.
+update_list <- function(x, new_x) {
+  for (n in names(new_x)) {
+    if (is.list(new_x[[n]])) {
+      x[[n]] <- update_list(x[[n]], new_x[[n]])
+    } else if (is.function(new_x[[n]]) && ! is.function(x[[n]])) {
+      x[[n]] <- new_x[[n]](x[[n]])
+    } else {
+      x[[n]] <- new_x[[n]]
+    }
+  }
+
+  return(x)
+}
+
+#' Module to modify the `param` list at specified time steps during the simulation
+#'
+#' @inheritParams aging_msm
+#'
+#' @details
+#' if a list `dat$param$param_updaters` is present, this module will update the
+#' `param` list with new values at given timesteps.
+#' `dat$control$param_updaters` is a list containing `updaters`. An updater is a
+#' list containing an `at` element telling when the changes will happend, an
+#' optional `verbose` boolean controlling whether to output a message when a
+#' change is made (default = TRUE) and a `param` list similar
+#' to the simulation's `dat$param` list.
+#' if the updater is a function but not the value to replace, the
+#' function will be applied to the current element (see example) .
+#'
+#' @examples
+#'  ## example of a `param_updaters` list
+#'  list(
+#'    list(
+#'      at = 4860,
+#'      param = list(
+#'        hiv.test.rate = rep(0.0128, 3),
+#'        trans.scale = c(1.61, 0.836, 0.622)
+#'      )
+#'    ),
+#'    list(
+#'      at = 5160,
+#'      verbose = FALSE,
+#'      param = list(
+#'        hiv.test.rate = function(x) x * 3,
+#'        trans.scale = function(x) x^2 / 3
+#'      )
+#'    )
+#'  )
+#'
+param_updater <- function(dat, at) {
+  if (is.null(dat$param$param_updaters))
+    return(dat)
+
+  param_updaters <- dat$param$param_updaters
+
+  for (i in seq_along(param_updaters)) {
+    if (at == param_updaters[[i]][["at"]]) {
+      verbose <- param_updaters[[i]][["verbose"]]
+      verbose <- if (is.null(verbose)) TRUE else verbose
+
+      new_params <- param_updaters[[i]][["param"]]
+
+      if (verbose) {
+        message(paste0(
+          "\n\nAt time step = ", at, " the `param` list was modified: \n"))
+        message(str(new_params))
+      }
+
+      dat$param <- update_list(dat$param, new_params)
+    }
+  }
+
+  return(dat)
+}
diff --git a/R/mod.departure.R b/R/mod.departure.R
index f40c620..6a3f95c 100644
--- a/R/mod.departure.R
+++ b/R/mod.departure.R
@@ -33,6 +33,10 @@ departure_msm <- function(dat, at) {
 
   aids.mr <- dat$param$aids.mr
   asmr <- dat$param$netstats$demog$asmr
+  asmr.row <- apply(as.matrix(age), 1,
+                      FUN = function(x) which(asmr[,"age"] == floor(x)))
+  asmr.temp <- asmr[asmr.row,]
+
 
   idsElig <- which(active == 1)
   rates <- rep(NA, length(idsElig))
@@ -40,7 +44,7 @@ departure_msm <- function(dat, at) {
   races <- sort(unique(race))
   for (i in seq_along(races)) {
     ids.race <- which(race == races[i])
-    rates[ids.race] <- asmr[age[ids.race], i + 1]
+    rates[ids.race] <- asmr.temp[ids.race, i + 1]
   }
   idsDep <- idsElig[rbinom(length(rates), 1, rates) == 1]
 
diff --git a/R/mod.hivtest.R b/R/mod.hivtest.R
index 3b25522..adbde52 100644
--- a/R/mod.hivtest.R
+++ b/R/mod.hivtest.R
@@ -87,14 +87,6 @@ hivtest_msm <- function(dat, at) {
 
   dat$epi$tot.neg.tests[at] <- length(tstNeg)
 
-  # number of new diagnoses by timing
-  dat$epi$newDx[at] <- length(tstPos)
-  diag.time <- dat$attr$diag.time
-  dat$epi$newDx45[at] <- length(intersect(tstPos, which(diag.time - inf.time <= 45/7)))
-  dat$epi$newDx140[at] <- length(intersect(tstPos, which(diag.time - inf.time <= 140/7)))
-  dat$epi$newDx200[at] <- length(intersect(tstPos, which(diag.time - inf.time <= 200/7)))
-  dat$epi$newDx2y[at] <- length(intersect(tstPos, which(diag.time - inf.time > 104)))
-
   return(dat)
 }
 
diff --git a/R/mod.hivtrans.R b/R/mod.hivtrans.R
index 18b122e..351e7dd 100644
--- a/R/mod.hivtrans.R
+++ b/R/mod.hivtrans.R
@@ -236,14 +236,7 @@ hivtrans_msm <- function(dat, at) {
     dat$attr$cuml.time.on.tx[infected] <- 0
     dat$attr$cuml.time.off.tx[infected] <- 0
 
-    # Attributes of transmitter
-    transmitter <- as.numeric(c(disc.ip[trans.ip == 1, 1],
-                                disc.rp[trans.rp == 1, 2]))
-    tab.trans <- table(transmitter)
-    uni.trans <- as.numeric(names(tab.trans))
-    dat$attr$count.trans[uni.trans] <- dat$attr$count.trans[uni.trans] +
-                                       as.numeric(tab.trans)
-   }
+  }
 
   # Summary Output
   dat$epi$incid[at] <- length(infected)
@@ -251,23 +244,6 @@ hivtrans_msm <- function(dat, at) {
   dat$epi$incid.H[at] <- sum(dat$attr$race[infected] == 2)
   dat$epi$incid.W[at] <- sum(dat$attr$race[infected] == 3)
 
-  if (length(infected) > 0) {
-    dat$epi$incid.undx[at] <- sum(dat$attr$diag.status[transmitter] == 0)
-    dat$epi$incid.dx[at] <- sum(dat$attr$diag.status[transmitter] == 1 &
-                                dat$attr$cuml.time.on.tx[transmitter] == 0)
-    dat$epi$incid.linked[at] <- sum(dat$attr$diag.status[transmitter] == 1 &
-                                    dat$attr$cuml.time.on.tx[transmitter] > 0 &
-                                    dat$attr$vl[transmitter] > log10(200))
-    dat$epi$incid.vsupp[at] <- sum(dat$attr$diag.status[transmitter] == 1 &
-                                   dat$attr$cuml.time.on.tx[transmitter] > 0 &
-                                   dat$attr$vl[transmitter] <= log10(200))
-  } else {
-    dat$epi$incid.undx[at] <- 0
-    dat$epi$incid.dx[at] <- 0
-    dat$epi$incid.linked[at] <- 0
-    dat$epi$incid.vsupp[at] <- 0
-  }
-
   return(dat)
 }
 
diff --git a/R/mod.prevalence.R b/R/mod.prevalence.R
index d2c0461..5c40512 100644
--- a/R/mod.prevalence.R
+++ b/R/mod.prevalence.R
@@ -27,7 +27,6 @@ prevalence_msm <- function(dat, at) {
   active <- dat$attr$active
   status <- dat$attr$status
   diag.status <- dat$attr$diag.status
-  diag.stage <- dat$attr$diag.stage
   diag.time <- dat$attr$diag.time
   aids.time <- dat$attr$aids.time
   inf.time <- dat$attr$inf.time
@@ -81,26 +80,17 @@ prevalence_msm <- function(dat, at) {
 
   # Care continuum stats (primary)
   dat$epi$cc.dx[at] <- sum(diag.status == 1 & inf.time >= 2, na.rm = TRUE) /
-                       sum(status == 1 & inf.time >= 2, na.rm = TRUE)
+    sum(status == 1 & inf.time >= 2, na.rm = TRUE)
   dat$epi$cc.dx.B[at] <- sum(diag.status == 1 & inf.time >= 2 & race == 1, na.rm = TRUE) /
-                         sum(status == 1 & inf.time >= 2 & race == 1, na.rm = TRUE)
+    sum(status == 1 & inf.time >= 2 & race == 1, na.rm = TRUE)
   dat$epi$cc.dx.H[at] <- sum(diag.status == 1 & inf.time >= 2 & race == 2, na.rm = TRUE) /
-                         sum(status == 1 & inf.time >= 2 & race == 2, na.rm = TRUE)
+    sum(status == 1 & inf.time >= 2 & race == 2, na.rm = TRUE)
   dat$epi$cc.dx.W[at] <- sum(diag.status == 1 & inf.time >= 2 & race == 3, na.rm = TRUE) /
-                         sum(status == 1 & inf.time >= 2 & race == 3, na.rm = TRUE)
+    sum(status == 1 & inf.time >= 2 & race == 3, na.rm = TRUE)
 
   dat$epi$cc.dx.aids[at] <- sum(diag.status == 1 & stage == 4 & inf.time >= 2 &
                                   aids.time - diag.time <= 52, na.rm = TRUE) /
                             sum(diag.status == 1 & inf.time >= 2, na.rm = TRUE)
-  dat$epi$cc.dx.aids.B[at] <- sum(diag.status == 1 & stage == 4 & inf.time >= 2 &
-                                    aids.time - diag.time <= 52 & race == 1, na.rm = TRUE) /
-                              sum(diag.status == 1 & inf.time >= 2 & race == 1, na.rm = TRUE)
-  dat$epi$cc.dx.aids.H[at] <- sum(diag.status == 1 & stage == 4 &  inf.time >= 2 &
-                                    aids.time - diag.time <= 52 & race == 2, na.rm = TRUE) /
-                              sum(diag.status == 1 & inf.time >= 2 & race == 2, na.rm = TRUE)
-  dat$epi$cc.dx.aids.W[at] <- sum(diag.status == 1 & stage == 4 & inf.time >= 2 &
-                                    aids.time - diag.time <= 52 & race == 3, na.rm = TRUE) /
-                              sum(diag.status == 1 & inf.time >= 2 & race == 3, na.rm = TRUE)
 
   dat$epi$cc.linked1m[at] <- sum(tx.init.time - diag.time <= 4 & diag.time >= 2, na.rm = TRUE) /
                              sum(dat$attr$diag.status == 1 & diag.time >= 2, na.rm = TRUE)
@@ -123,13 +113,13 @@ prevalence_msm <- function(dat, at) {
   dat$epi$cc.vsupp[at] <- sum(vl <= log10(200) & diag.status == 1 & diag.time >= 2, na.rm = TRUE) /
                           sum(diag.status == 1 & diag.time >= 2, na.rm = TRUE)
   dat$epi$cc.vsupp.B[at] <- sum(vl <= log10(200) & diag.status == 1 &
-                                diag.time >= 2 & race == 1, na.rm = TRUE) /
+                                  diag.time >= 2 & race == 1, na.rm = TRUE) /
                             sum(diag.status == 1 & diag.time >= 2 & race == 1, na.rm = TRUE)
   dat$epi$cc.vsupp.H[at] <- sum(vl <= log10(200) & diag.status == 1 &
-                                diag.time >= 2 & race == 2, na.rm = TRUE) /
+                                  diag.time >= 2 & race == 2, na.rm = TRUE) /
                             sum(diag.status == 1 & diag.time >= 2 & race == 2, na.rm = TRUE)
   dat$epi$cc.vsupp.W[at] <- sum(vl <= log10(200) & diag.status == 1 &
-                                diag.time >= 2 & race == 3, na.rm = TRUE) /
+                                  diag.time >= 2 & race == 3, na.rm = TRUE) /
                             sum(diag.status == 1 & diag.time >= 2 & race == 3, na.rm = TRUE)
 
   dat$epi$cc.vsupp.all[at] <- sum(vl <= log10(200) & status == 1 & inf.time >= 2, na.rm = TRUE) /
@@ -141,103 +131,8 @@ prevalence_msm <- function(dat, at) {
   dat$epi$cc.vsupp.all.W[at] <- sum(vl <= log10(200) & status == 1 & inf.time >= 2 & race == 3, na.rm = TRUE) /
                                 sum(status == 1 & inf.time >= 2 & race == 3, na.rm = TRUE)
 
-  dat$epi$cc.vsupp.dur1y[at] <- 1 - (sum((at - vl.last.usupp) <= 52 & diag.time >= 2 &
-                                         diag.status == 1, na.rm = TRUE) /
-                                     sum(diag.status == 1 & diag.time >= 2, na.rm = TRUE))
-  dat$epi$cc.vsupp.dur1y.B[at] <- 1 - (sum((at - vl.last.usupp) <= 52 & diag.time >= 2 &
-                                         diag.status == 1 & race == 1, na.rm = TRUE) /
-                                     sum(diag.status == 1 & diag.time >= 2 & race == 1, na.rm = TRUE))
-  dat$epi$cc.vsupp.dur1y.H[at] <- 1 - (sum((at - vl.last.usupp) <= 52 & diag.time >= 2 &
-                                           diag.status == 1 & race == 2, na.rm = TRUE) /
-                                       sum(diag.status == 1 & diag.time >= 2 & race == 2, na.rm = TRUE))
-  dat$epi$cc.vsupp.dur1y.W[at] <- 1 - (sum((at - vl.last.usupp) <= 52 & diag.time >= 2 &
-                                           diag.status == 1 & race == 3, na.rm = TRUE) /
-                                       sum(diag.status == 1 & diag.time >= 2 & race == 3, na.rm = TRUE))
-
   dat$epi$cc.HIV.mr[at] <- (dat$epi$dep.HIV[at]/dat$epi$i.num[at])*52
 
-  # Care continuum stats (secondary)
-  dat$epi$cc.test.int[at] <- mean(at - last.neg.test[diag.status == 0], na.rm = TRUE)
-  dat$epi$cc.test.int.B[at] <- mean(at - last.neg.test[diag.status == 0 & race == 1], na.rm = TRUE)
-  dat$epi$cc.test.int.H[at] <- mean(at - last.neg.test[diag.status == 0 & race == 2], na.rm = TRUE)
-  dat$epi$cc.test.int.W[at] <- mean(at - last.neg.test[diag.status == 0 & race == 3], na.rm = TRUE)
-
-  dat$epi$cc.dx.delay[at] <- mean(diag.time[diag.time >= 2] - inf.time[diag.time >= 2], na.rm = TRUE)
-  dat$epi$cc.dx.delay.B[at] <- mean(diag.time[diag.time >= 2 & race == 1] -
-                                      inf.time[diag.time >= 2 & race == 1], na.rm = TRUE)
-  dat$epi$cc.dx.delay.H[at] <- mean(diag.time[diag.time >= 2 & race == 2] -
-                                      inf.time[diag.time >= 2 & race == 2], na.rm = TRUE)
-  dat$epi$cc.dx.delay.W[at] <- mean(diag.time[diag.time >= 2 & race == 3] -
-                                      inf.time[diag.time >= 2 & race == 3], na.rm = TRUE)
-
-  dat$epi$cc.dx.delay.int[at] <- mean(diag.time[diag.time >= 3380] - inf.time[diag.time >= 3380], na.rm = TRUE)
-  dat$epi$cc.dx.delay.int.B[at] <- mean(diag.time[diag.time >= 3380 & race == 1] -
-                                      inf.time[diag.time >= 3380 & race == 1], na.rm = TRUE)
-  dat$epi$cc.dx.delay.int.H[at] <- mean(diag.time[diag.time >= 3380 & race == 2] -
-                                      inf.time[diag.time >= 3380 & race == 2], na.rm = TRUE)
-  dat$epi$cc.dx.delay.int.W[at] <- mean(diag.time[diag.time >= 3380 & race == 3] -
-                                      inf.time[diag.time >= 3380 & race == 3], na.rm = TRUE)
-
-  # same as above, but with medians
-  dat$epi$cc.dx.delay.med[at] <- median(diag.time[diag.time >= 2] - inf.time[diag.time >= 2], na.rm = TRUE)
-  dat$epi$cc.dx.delay.B.med[at] <- median(diag.time[diag.time >= 2 & race == 1] -
-                                      inf.time[diag.time >= 2 & race == 1], na.rm = TRUE)
-  dat$epi$cc.dx.delay.H.med[at] <- median(diag.time[diag.time >= 2 & race == 2] -
-                                      inf.time[diag.time >= 2 & race == 2], na.rm = TRUE)
-  dat$epi$cc.dx.delay.W.med[at] <- median(diag.time[diag.time >= 2 & race == 3] -
-                                      inf.time[diag.time >= 2 & race == 3], na.rm = TRUE)
-
-  dat$epi$cc.dx.delay.int.med[at] <- median(diag.time[diag.time >= 3380] - inf.time[diag.time >= 3380], na.rm = TRUE)
-  dat$epi$cc.dx.delay.int.B.med[at] <- median(diag.time[diag.time >= 3380 & race == 1] -
-                                          inf.time[diag.time >= 3380 & race == 1], na.rm = TRUE)
-  dat$epi$cc.dx.delay.int.H.med[at] <- median(diag.time[diag.time >= 3380 & race == 2] -
-                                          inf.time[diag.time >= 3380 & race == 2], na.rm = TRUE)
-  dat$epi$cc.dx.delay.int.W.med[at] <- median(diag.time[diag.time >= 3380 & race == 3] -
-                                          inf.time[diag.time >= 3380 & race == 3], na.rm = TRUE)
-
-  # dat$epi$cc.tx.any1y[at] <- sum((at - dat$attr$tx.period.last <= 52), na.rm = TRUE) /
-  #   sum(dat$attr$diag.status == 1, na.rm = TRUE)
-  # dat$epi$cc.tx.any1y.B[at] <- sum((at - dat$attr$tx.period.last <= 52) & race == 1, na.rm = TRUE) /
-  #   sum(dat$attr$diag.status == 1 & race == 1, na.rm = TRUE)
-  # dat$epi$cc.tx.any1y.H[at] <- sum((at - dat$attr$tx.period.last <= 52) & race == 2, na.rm = TRUE) /
-  #   sum(dat$attr$diag.status == 1 & race == 2, na.rm = TRUE)
-  # dat$epi$cc.tx.any1y.W[at] <- sum((at - dat$attr$tx.period.last <= 52) & race == 3, na.rm = TRUE) /
-  #   sum(dat$attr$diag.status == 1 & race == 3, na.rm = TRUE)
-
-  # dat$epi$cc.dx.delay[at] <- mean(dat$attr$diag.time - dat$attr$inf.time, na.rm = TRUE)
-  # dat$epi$cc.testpy[at] <- 1-sum((at - dat$attr$last.neg.test) > 52 & status == 0,
-  #     is.na(dat$attr$last.neg.test) & status == 0, na.rm = TRUE) /
-  #   sum(status == 0)
-  # dat$epi$cc.linked[at] <- sum(dat$attr$cuml.time.on.tx > 0, na.rm = TRUE) /
-  #   sum(dat$attr$diag.status == 1, na.rm = TRUE)
-  # dat$epi$cc.tx[at] <- sum(dat$attr$tx.status == 1, na.rm = TRUE) /
-  #   sum(dat$attr$diag.status == 1, na.rm = TRUE)
-  # dat$epi$cc.tx.ret3m[at] <- sum((at - dat$attr$tx.period.last) <= 52 &
-  #       (dat$attr$tx.period.last - dat$attr$tx.period.first) > 13, na.rm = TRUE) /
-  #   sum(dat$attr$diag.status == 1, na.rm = TRUE)
-  # dat$epi$cc.vsupp.tt1[at] <- sum(dat$attr$vl <= log10(200) & dat$attr$tt.traj == 1, na.rm = TRUE) /
-  #   sum(dat$attr$diag.status == 1 & dat$attr$tt.traj == 1, na.rm = TRUE)
-  # dat$epi$cc.vsupp.tt2[at] <- sum(dat$attr$vl <= log10(200) & dat$attr$tt.traj == 2, na.rm = TRUE) /
-  #   sum(dat$attr$diag.status == 1 & dat$attr$tt.traj == 2, na.rm = TRUE)
-  # dat$epi$cc.vsupp.tt3[at] <- sum(dat$attr$vl <= log10(200) & dat$attr$tt.traj == 3, na.rm = TRUE) /
-  #   sum(dat$attr$diag.status == 1 & dat$attr$tt.traj == 3, na.rm = TRUE)
-
-
-  # HIV screening outcomes
-  dat$epi$mean.neg.tests[at] <- mean(dat$attr$num.neg.tests[diag.status == 0], na.rm = TRUE)
-  dat$epi$mean.neg.tests.B[at] <- mean(dat$attr$num.neg.tests[diag.status == 0 & race == 1], na.rm = TRUE)
-  dat$epi$mean.neg.tests.H[at] <- mean(dat$attr$num.neg.tests[diag.status == 0 & race == 2], na.rm = TRUE)
-  dat$epi$mean.neg.tests.W[at] <- mean(dat$attr$num.neg.tests[diag.status == 0 & race == 3], na.rm = TRUE)
-
-  dat$epi$test.past.year[at] <- sum(at - dat$attr$last.neg.test <= 52 & diag.status == 0, na.rm = TRUE) /
-    sum(diag.status == 0, na.rm = TRUE)
-  dat$epi$test.past.year.B[at] <- sum(at - dat$attr$last.neg.test <= 52 & diag.status == 0 & race == 1, na.rm = TRUE) /
-    sum(diag.status == 0 & race == 1, na.rm = TRUE)
-  dat$epi$test.past.year.H[at] <- sum(at - dat$attr$last.neg.test <= 52 & diag.status == 0 & race == 2, na.rm = TRUE) /
-    sum(diag.status == 0 & race == 2, na.rm = TRUE)
-  dat$epi$test.past.year.W[at] <- sum(at - dat$attr$last.neg.test <= 52 & diag.status == 0 & race == 3, na.rm = TRUE) /
-    sum(diag.status == 0 & race == 3, na.rm = TRUE)
-
   # HIV stage
   dat$epi$hstage.acute[at] <- sum(stage %in% 1:2 & diag.time >= 2, na.rm = TRUE) /
                               sum(status == 1 & diag.time >= 2, na.rm = TRUE)
@@ -269,6 +164,56 @@ prevalence_msm <- function(dat, at) {
   dat$epi$ir100.ct[at] <- ir100.rct + ir100.uct
   dat$epi$ir100.sti[at] <- dat$epi$ir100.gc[at] + dat$epi$ir100.ct[at]
 
+
+  # STI tx
+  dat$epi$gc.tx[at] <- sum(dat$attr$uGC.tx | dat$attr$rGC.tx, na.rm = TRUE)
+  dat$epi$ct.tx[at] <- sum(dat$attr$uCT.tx | dat$attr$rCT.tx, na.rm = TRUE)
+
+  dat$epi$gc.tx.B[at] <- sum(race == 1 & (dat$attr$uGC.tx | dat$attr$rGC.tx), na.rm = TRUE)
+  dat$epi$ct.tx.B[at] <- sum(race == 1 & (dat$attr$uCT.tx | dat$attr$rCT.tx), na.rm = TRUE)
+
+  dat$epi$gc.tx.H[at] <- sum(race == 2 & (dat$attr$uGC.tx | dat$attr$rGC.tx), na.rm = TRUE)
+  dat$epi$ct.tx.H[at] <- sum(race == 2 & (dat$attr$uCT.tx | dat$attr$rCT.tx), na.rm = TRUE)
+
+  dat$epi$gc.tx.W[at] <- sum(race == 3 & (dat$attr$uGC.tx | dat$attr$rGC.tx), na.rm = TRUE)
+  dat$epi$ct.tx.W[at] <- sum(race == 3 & (dat$attr$uCT.tx | dat$attr$rCT.tx), na.rm = TRUE)
+
+  gc <- rGC == 1 | uGC == 1
+  ct <- rCT == 1 | uCT == 1
+
+  dat$epi$gc[at] <- sum(gc, na.rm = TRUE)
+  dat$epi$ct[at] <- sum(ct, na.rm = TRUE)
+
+  dat$epi$gc.B[at] <- sum(race == 1 & gc, na.rm = TRUE)
+  dat$epi$ct.B[at] <- sum(race == 1 & ct, na.rm = TRUE)
+
+  dat$epi$gc.H[at] <- sum(race == 2 & gc, na.rm = TRUE)
+  dat$epi$ct.H[at] <- sum(race == 2 & ct, na.rm = TRUE)
+
+  dat$epi$gc.W[at] <- sum(race == 3 & gc, na.rm = TRUE)
+  dat$epi$ct.W[at] <- sum(race == 3 & ct, na.rm = TRUE)
+
+  # Mean degrees
+  main.deg <- get_degree(dat$el[[1]])
+  casl.deg <- get_degree(dat$el[[2]])
+  inst.deg <- get_degree(dat$el[[3]])
+
+  dat$epi$main.deg[at] <- mean(main.deg, na.rm = TRUE)
+  dat$epi$casl.deg[at] <- mean(casl.deg, na.rm = TRUE)
+  dat$epi$inst.deg[at] <- mean(inst.deg, na.rm = TRUE)
+
+  dat$epi$main.deg.B[at] <- mean(main.deg[which(race == 1)], na.rm = TRUE)
+  dat$epi$casl.deg.B[at] <- mean(casl.deg[which(race == 1)], na.rm = TRUE)
+  dat$epi$inst.deg.B[at] <- mean(inst.deg[which(race == 1)], na.rm = TRUE)
+
+  dat$epi$main.deg.H[at] <- mean(main.deg[which(race == 2)], na.rm = TRUE)
+  dat$epi$casl.deg.H[at] <- mean(casl.deg[which(race == 2)], na.rm = TRUE)
+  dat$epi$inst.deg.H[at] <- mean(inst.deg[which(race == 2)], na.rm = TRUE)
+
+  dat$epi$main.deg.W[at] <- mean(main.deg[which(race == 3)], na.rm = TRUE)
+  dat$epi$casl.deg.W[at] <- mean(casl.deg[which(race == 3)], na.rm = TRUE)
+  dat$epi$inst.deg.W[at] <- mean(inst.deg[which(race == 3)], na.rm = TRUE)
+
   return(dat)
 }
 
diff --git a/R/mod.simnet.R b/R/mod.simnet.R
index fff4009..95fe5e3 100644
--- a/R/mod.simnet.R
+++ b/R/mod.simnet.R
@@ -20,6 +20,11 @@ simnet_msm <- function(dat, at) {
   ## Main network
   nwparam.m <- EpiModel::get_nwparam(dat, network = 1)
 
+  nwparam.m$coef.form[1] <-
+    nwparam.m$coef.form[1] + log(dat$param$netresim.form.rr[1])
+  nwparam.m$coef.diss$coef.adj[1] <-
+    nwparam.m$coef.diss$coef.adj[1] - log(dat$param$netresim.disl.rr[1])
+
   dat$attr$deg.casl <- get_degree(dat$el[[2]])
   dat <- tergmLite::updateModelTermInputs(dat, network = 1)
 
@@ -31,10 +36,14 @@ simnet_msm <- function(dat, at) {
 
   plist1 <- update_plist(dat, at, ptype = 1)
 
-
   ## Casual network
   nwparam.p <- EpiModel::get_nwparam(dat, network = 2)
 
+  nwparam.p$coef.form[1] <-
+    nwparam.p$coef.form[1] + log(dat$param$netresim.form.rr[2])
+  nwparam.p$coef.diss$coef.adj[1] <-
+    nwparam.p$coef.diss$coef.adj[1] - log(dat$param$netresim.disl.rr[2])
+
   dat$attr$deg.main <- get_degree(dat$el[[1]])
   dat <- tergmLite::updateModelTermInputs(dat, network = 2)
 
@@ -54,6 +63,8 @@ simnet_msm <- function(dat, at) {
 
   ## One-off network
   nwparam.i <- EpiModel::get_nwparam(dat, network = 3)
+  nwparam.i$coef.form[1] <-
+    nwparam.i$coef.form[1] + log(dat$param$netresim.form.rr[3])
 
   dat$attr$deg.tot <- pmin(dat$attr$deg.main + get_degree(dat$el[[2]]), 3)
   dat <- tergmLite::updateModelTermInputs(dat, network = 3)
@@ -86,7 +97,9 @@ update_plist <- function(dat, at, ptype) {
 
   # look up new formations, row bind them
   news_uid_start <- news_uid[news_uid[, 3] == 1, , drop = FALSE]
-  plist1 <- rbind(plist1, cbind(news_uid_start[, 1:2, drop = FALSE], ptype, at, NA))
+  if (nrow(news_uid_start) > 0) {
+    plist1 <- rbind(plist1, cbind(news_uid_start[, 1:2, drop = FALSE], ptype, at, NA))
+  }
 
   return(plist1)
 }
diff --git a/R/mod.stitrans.R b/R/mod.stitrans.R
index 414e436..1c2803c 100644
--- a/R/mod.stitrans.R
+++ b/R/mod.stitrans.R
@@ -300,6 +300,20 @@ stitrans_msm <- function(dat, at) {
   dat$epi$incid.ct.B[at] <- length(intersect(union(idsInf_rct, idsInf_uct), which(race == 2)))
   dat$epi$incid.ct.W[at] <- length(intersect(union(idsInf_rct, idsInf_uct), which(race == 3)))
 
+  dat$epi$incid.rct.B[at] <- length(intersect(idsInf_rct, which(race == 1)))
+  dat$epi$incid.rct.H[at] <- length(intersect(idsInf_rct, which(race == 2)))
+  dat$epi$incid.rct.W[at] <- length(intersect(idsInf_rct, which(race == 3)))
+  dat$epi$incid.uct.B[at] <- length(intersect(idsInf_uct, which(race == 1)))
+  dat$epi$incid.uct.H[at] <- length(intersect(idsInf_uct, which(race == 2)))
+  dat$epi$incid.uct.W[at] <- length(intersect(idsInf_uct, which(race == 3)))
+
+  dat$epi$incid.rgc.B[at] <- length(intersect(idsInf_rgc, which(race == 1)))
+  dat$epi$incid.rgc.H[at] <- length(intersect(idsInf_rgc, which(race == 2)))
+  dat$epi$incid.rgc.W[at] <- length(intersect(idsInf_rgc, which(race == 3)))
+  dat$epi$incid.ugc.B[at] <- length(intersect(idsInf_ugc, which(race == 1)))
+  dat$epi$incid.ugc.H[at] <- length(intersect(idsInf_ugc, which(race == 2)))
+  dat$epi$incid.ugc.W[at] <- length(intersect(idsInf_ugc, which(race == 3)))
+
   # Check all infected have all STI attributes
   stopifnot(all(!is.na(rGC.infTime[rGC == 1])),
             all(!is.na(rGC.sympt[rGC == 1])),
diff --git a/R/params.R b/R/params.R
index 2a4301f..66e63ca 100644
--- a/R/params.R
+++ b/R/params.R
@@ -93,7 +93,8 @@
 #'
 #' @param epistats GLMs for epidemiological parameter from the standard ARTnet workflow.
 #' @param acts.aids.vl Viral load level after which sexual act rate goes to zero.
-#' @param acts.scale Scalar for main/casual act rate for model calibration.
+#' @param acts.scale.main Scalar for main act rate for model calibration.
+#' @param acts.scale.casl Scalar for casual act rate for model calibration.
 #' @param cond.scale Scalar for condom use probability for model calibration.
 #'
 #' @param riskh.start Time step at which behavioral risk history assessment occurs.
@@ -221,7 +222,8 @@ param_msm <- function(netstats,
                       # Behavioral
                       epistats,
                       acts.aids.vl = 5.75,
-                      acts.scale = 1,
+                      acts.scale.main = 1,
+                      acts.scale.casl = 1,
                       cond.scale = 1,
 
                       # STI epi
@@ -369,6 +371,7 @@ control_msm <- function(simno = 1,
                         nsteps = 100,
                         start = 1,
                         initialize.FUN = initialize_msm,
+                        param_updater.FUN = param_updater,
                         aging.FUN = aging_msm,
                         departure.FUN = departure_msm,
                         arrival.FUN = arrival_msm,
@@ -410,6 +413,7 @@ control_msm <- function(simno = 1,
 
   p$save.network <- FALSE
   p$verbose.int <- 1
+  p$tergmLite <- TRUE
 
   class(p) <- "control.net"
   return(p)
diff --git a/inst/msm-test-script.R b/inst/msm-test-script.R
index ae7fcfb..93c2fd0 100644
--- a/inst/msm-test-script.R
+++ b/inst/msm-test-script.R
@@ -3,40 +3,97 @@ rm(list = ls())
 suppressMessages(library("EpiModelHIV"))
 devtools::load_all("~/Dropbox/Dev/EpiModelHIV/EpiModelHIV-p")
 
-scr.dir <- "~/Dropbox/Dev/ARTnet/"
-netstats <- readRDS(file.path(scr.dir, "data/artnet.NetStats.Atlanta.rda"))
-epistats <- readRDS(file.path(scr.dir, "data/artnet.EpiStats.Atlanta.rda"))
-est <- readRDS(file.path(scr.dir, "data/artnet.NetEst.Atlanta.rda"))
-
-param <- param_msm(netstats = netstats,
-                   epistats = epistats,
-                   hiv.test.int = c(43, 43, 45),
-                   a.rate = 0.00055,
-                   riskh.start = 2,
-                   prep.start = 30,
-                   prep.start.prob = 0.10,
-                   tt.part.supp = c(0.20, 0.20, 0.20),
-                   tt.full.supp = c(0.40, 0.40, 0.40),
-                   tt.dur.supp = c(0.40, 0.40, 0.40),
-                   tx.halt.full.rr = 0.8,
-                   tx.halt.dur.rr = 0.1,
-                   tx.reinit.full.rr = 2.0,
-                   tx.reinit.dur.rr = 5.0,
-                   hiv.rgc.rr = 2.5,
-                   hiv.ugc.rr = 1.5,
-                   hiv.rct.rr = 2.5,
-                   hiv.uct.rr = 1.5,
-                   hiv.dual.rr = 0.0,
-                   rgc.tprob = 0.35,
-                   ugc.tprob = 0.25,
-                   rct.tprob = 0.20,
-                   uct.tprob = 0.16,
-                   rgc.ntx.int = 16.8,
-                   ugc.ntx.int = 16.8,
-                   rct.ntx.int = 32,
-                   uct.ntx.int = 32,
-                   acts.aids.vl = 5.75)
-init <- init_msm()
+scr.dir <- "~/Dropbox/Projects/SexualDistancing/"
+netstats <- readRDS(file.path(scr.dir, "out/est/netstats.rds"))
+epistats <- readRDS(file.path(scr.dir, "out/est/epistats.rds"))
+est <- readRDS(file.path(scr.dir, "out/est/netest.rds"))
+
+full_tx_eff <- rep(1, 3)
+
+param <- param_msm(
+  netstats = netstats,
+  epistats = epistats,
+  hiv.test.rate = c(0.00385, 0.00385, 0.0069),
+  hiv.test.late.prob = rep(0, 3),
+  tx.init.prob = c(0.1775, 0.19, 0.2521),
+  tt.part.supp = 1 - full_tx_eff,
+  tt.full.supp = full_tx_eff,
+  tt.dur.supp = rep(0, 3),
+  tx.halt.part.prob = c(0.0065, 0.0053, 0.003),
+  tx.halt.full.rr = rep(0.45, 3),
+  tx.halt.dur.rr = rep(0.45, 3),
+  tx.reinit.part.prob = rep(0.00255, 3),
+  tx.reinit.full.rr = rep(1, 3),
+  tx.reinit.dur.rr = rep(1, 3),
+  max.time.off.tx.full.int = 52 * 15,
+  max.time.on.tx.part.int = 52 * 10,
+  max.time.off.tx.part.int = 52 * 10,
+  aids.mr = 1 / 250,
+  trans.scale =  c(2.7, 0.35, 0.243), #c(2.21, 0.405, 0.255),
+  acts.scale.main = 1.00,
+  acts.scale.casl = 0.10,
+  acts.aids.vl = 5.75,
+  circ.prob = c(0.874, 0.874, 0.918),
+  a.rate = 0.00052,
+  prep.start = (52 * 60) + 1,
+  riskh.start = 52 * 59,
+  prep.adhr.dist = c(0.089, 0.127, 0.784),
+  prep.adhr.hr = c(0.69, 0.19, 0.01),
+  prep.start.prob =  0.71, # 0.00896,
+  prep.discont.rate = 0.02138792, # 1 - (2^(-1/(224.4237/7)))
+  ## prep.tst.int = 90/7,         # do I need that?
+  ## prep.risk.int = 182/7,       # do I need that?
+  ## prep.sti.screen.int = 182/7,
+  ## prep.sti.prob.tx = 1,
+  prep.risk.reassess.method = "year",
+  prep.require.lnt = TRUE, # FALSE -> start with random PrEP initiation
+
+  ## STI PARAMS (default: from combprev2, make it gaps)
+  ## Using values in prep-race: scripts/burnin/sim.burn.R
+  ## If not mentionned -> default from prep disparities
+  ## for H : mean(c(B, W))
+  #ok
+  rgc.tprob = 0.2267303, #logistic(logit(0.19) + log(1.25)) #0.357,  # gaps appendix 9.4
+  ugc.tprob = 0.19,# 0.248,  # gaps appendix 9.4
+  rct.tprob = 0.2038369, #logistic(logit(0.17) + log(1.25)) #0.3216, # gaps appendix 9.3
+  uct.tprob = 0.17,#0.213,  # gaps appendix 9.3
+  rgc.sympt.prob = 0.1,#0.077, # gaps appendix 10.3
+  ugc.sympt.prob = 0.9333333,#0.824, # gaps appendix 10.3
+  rct.sympt.prob = 0.1,#0.1035,# gaps appendix 10.2
+  uct.sympt.prob = 0.95,#0.885, # gaps appendix 10.2
+  rgc.ntx.int = 26,#35.11851, # gaps appendix 11.2
+  ugc.ntx.int = 26,#35.11851, # gaps appendix 11.2
+  gc.tx.int   = 2, # gaps appendix 11.2 - mentionned, not explicit
+  rct.ntx.int = 32,#44.24538, # gaps appendix 11.1
+  uct.ntx.int = 32,#44.24538, # gaps appendix 11.1
+  ct.tx.int   = 2, # gaps appendix 11.1 - mentionned, not explicit
+
+  gc.sympt.prob.tx =  rep(0.9, 3),  #c(0.86, 0.91, 0.96),
+  ct.sympt.prob.tx =  rep(0.9, 3),  #c(0.72, 0.785, 0.85),
+  gc.asympt.prob.tx = rep(0.1, 3), #c(0.10, 0.145, 0.19),
+  ct.asympt.prob.tx = rep(0.1, 3), #c(0.05, 0.525, 0.10),
+  # gaps appendix 9.3 - 9.4 (not explained this way but similar result)
+  sti.cond.eff = 0.95,
+  sti.cond.fail = c(0.39, 0.3, 0.21),
+  # gaps appendix 9.2
+  hiv.rgc.rr = 2.78,
+  hiv.ugc.rr = 1.73,
+  hiv.rct.rr = 2.78,
+  hiv.uct.rr = 1.73,
+  # if both ct + gc -> log(RRgc) + 0.2 * log(RRct) | swap ct and gc if RRct > RRgc
+  hiv.dual.rr = 0.2, # not mentionned in appendix
+
+  netresim.form.rr = rep(1, 3),
+  netresim.disl.rr = rep(1, 2),
+
+)
+
+init <- init_msm(
+  prev.ugc = 0.05,
+  prev.rct = 0.05,
+  prev.rgc = 0.05,
+  prev.uct = 0.05
+)
 control <- control_msm(simno = 1,
                        nsteps = 52 * 2,
                        nsims = 1,
@@ -44,6 +101,7 @@ control <- control_msm(simno = 1,
                        save.nwstats = FALSE,
                        save.clin.hist = FALSE)
 
+debug(acts_msm)
 sim <- netsim(est, param, init, control)
 
 # Explore clinical history
diff --git a/man/control_het.Rd b/man/control_het.Rd
index 23515a4..e3a37a8 100644
--- a/man/control_het.Rd
+++ b/man/control_het.Rd
@@ -5,15 +5,32 @@
 \title{Control Settings for Stochastic Network Model of HIV-1 Infection in
        Sub-Saharan Africa}
 \usage{
-control_het(simno = 1, nsteps = 100, start = 1, nsims = 1,
-  ncores = 1, par.type = "single", initialize.FUN = initialize_het,
-  aging.FUN = aging_het, cd4.FUN = cd4_het, vl.FUN = vl_het,
-  dx.FUN = dx_het, tx.FUN = tx_het, deaths.FUN = deaths_het,
-  births.FUN = births_het, resim_nets.FUN = simnet_het,
-  trans.FUN = trans_het, prev.FUN = prevalence_het,
-  verbose.FUN = verbose.net, module.order = NULL,
-  save.nwstats = FALSE, save.other = c("el", "attr"), verbose = TRUE,
-  skip.check = TRUE, ...)
+control_het(
+  simno = 1,
+  nsteps = 100,
+  start = 1,
+  nsims = 1,
+  ncores = 1,
+  par.type = "single",
+  initialize.FUN = initialize_het,
+  aging.FUN = aging_het,
+  cd4.FUN = cd4_het,
+  vl.FUN = vl_het,
+  dx.FUN = dx_het,
+  tx.FUN = tx_het,
+  deaths.FUN = deaths_het,
+  births.FUN = births_het,
+  resim_nets.FUN = simnet_het,
+  trans.FUN = trans_het,
+  prev.FUN = prevalence_het,
+  verbose.FUN = verbose.net,
+  module.order = NULL,
+  save.nwstats = FALSE,
+  save.other = c("el", "attr"),
+  verbose = TRUE,
+  skip.check = TRUE,
+  ...
+)
 }
 \arguments{
 \item{simno}{Simulation ID number.}
diff --git a/man/control_msm.Rd b/man/control_msm.Rd
index 44b017c..4f40d80 100644
--- a/man/control_msm.Rd
+++ b/man/control_msm.Rd
@@ -4,18 +4,38 @@
 \alias{control_msm}
 \title{Epidemic Model Control Settings}
 \usage{
-control_msm(simno = 1, nsims = 1, ncores = 1, nsteps = 100,
-  start = 1, initialize.FUN = initialize_msm, aging.FUN = aging_msm,
-  departure.FUN = departure_msm, arrival.FUN = arrival_msm,
-  hivtest.FUN = hivtest_msm, hivtx.FUN = hivtx_msm,
-  hivprogress.FUN = hivprogress_msm, hivvl.FUN = hivvl_msm,
-  resim_nets.FUN = simnet_msm, acts.FUN = acts_msm,
-  condoms.FUN = condoms_msm, position.FUN = position_msm,
-  prep.FUN = prep_msm, hivtrans.FUN = hivtrans_msm,
-  stitrans.FUN = stitrans_msm, stirecov.FUN = stirecov_msm,
-  stitx.FUN = stitx_msm, prev.FUN = prevalence_msm,
-  verbose.FUN = verbose.net, save.nwstats = FALSE,
-  save.clin.hist = FALSE, truncate.plist = TRUE, verbose = TRUE, ...)
+control_msm(
+  simno = 1,
+  nsims = 1,
+  ncores = 1,
+  nsteps = 100,
+  start = 1,
+  initialize.FUN = initialize_msm,
+  param_updater.FUN = param_updater,
+  aging.FUN = aging_msm,
+  departure.FUN = departure_msm,
+  arrival.FUN = arrival_msm,
+  hivtest.FUN = hivtest_msm,
+  hivtx.FUN = hivtx_msm,
+  hivprogress.FUN = hivprogress_msm,
+  hivvl.FUN = hivvl_msm,
+  resim_nets.FUN = simnet_msm,
+  acts.FUN = acts_msm,
+  condoms.FUN = condoms_msm,
+  position.FUN = position_msm,
+  prep.FUN = prep_msm,
+  hivtrans.FUN = hivtrans_msm,
+  stitrans.FUN = stitrans_msm,
+  stirecov.FUN = stirecov_msm,
+  stitx.FUN = stitx_msm,
+  prev.FUN = prevalence_msm,
+  verbose.FUN = verbose.net,
+  save.nwstats = FALSE,
+  save.clin.hist = FALSE,
+  truncate.plist = TRUE,
+  verbose = TRUE,
+  ...
+)
 }
 \arguments{
 \item{simno}{Unique ID for the simulation run, used for file naming purposes
diff --git a/man/init_het.Rd b/man/init_het.Rd
index d1c38ba..71d1599 100644
--- a/man/init_het.Rd
+++ b/man/init_het.Rd
@@ -5,9 +5,15 @@
 \title{Initial Conditions for Stochastic Network Model of HIV-1 Infection in
        Sub-Saharan Africa}
 \usage{
-init_het(i.prev.male = 0.05, i.prev.feml = 0.05, ages.male = seq(18,
-  55, 7/365), ages.feml = seq(18, 55, 7/365),
-  inf.time.dist = "geometric", max.inf.time = 5 * 365, ...)
+init_het(
+  i.prev.male = 0.05,
+  i.prev.feml = 0.05,
+  ages.male = seq(18, 55, 7/365),
+  ages.feml = seq(18, 55, 7/365),
+  inf.time.dist = "geometric",
+  max.inf.time = 5 * 365,
+  ...
+)
 }
 \arguments{
 \item{i.prev.male}{Prevalence of initially infected males.}
diff --git a/man/init_msm.Rd b/man/init_msm.Rd
index e12cbad..8710d6c 100644
--- a/man/init_msm.Rd
+++ b/man/init_msm.Rd
@@ -4,8 +4,13 @@
 \alias{init_msm}
 \title{Epidemic Model Initial Conditions}
 \usage{
-init_msm(prev.ugc = 0.005, prev.rgc = 0.005, prev.uct = 0.013,
-  prev.rct = 0.013, ...)
+init_msm(
+  prev.ugc = 0.005,
+  prev.rgc = 0.005,
+  prev.uct = 0.013,
+  prev.rct = 0.013,
+  ...
+)
 }
 \arguments{
 \item{prev.ugc}{Initial prevalence of urethral gonorrhea.}
diff --git a/man/param_het.Rd b/man/param_het.Rd
index 12fbee9..d7e42a5 100644
--- a/man/param_het.Rd
+++ b/man/param_het.Rd
@@ -5,19 +5,45 @@
 \title{Parameters for Stochastic Network Model of HIV-1 Infection in
        Sub-Saharan Africa}
 \usage{
-param_het(time.unit = 7, acute.stage.mult = 5, aids.stage.mult = 1,
-  vl.acute.topeak = 14, vl.acute.toset = 107, vl.acute.peak = 6.7,
-  vl.setpoint = 4.5, vl.aidsmax = 7, cond.prob = 0.09,
-  cond.eff = 0.78, act.rate.early = 0.362, act.rate.late = 0.197,
-  act.rate.cd4 = 50, acts.rand = TRUE, circ.prob.birth = 0.9,
-  circ.eff = 0.53, tx.elig.cd4 = 350, tx.init.cd4.mean = 120,
-  tx.init.cd4.sd = 40, tx.adhere.full = 0.76, tx.adhere.part = 0.5,
-  tx.vlsupp.time = 365/3, tx.vlsupp.level = 1.5,
-  tx.cd4.recrat.feml = 11.6/30, tx.cd4.recrat.male = 9.75/30,
-  tx.cd4.decrat.feml = 11.6/30, tx.cd4.decrat.male = 9.75/30,
-  tx.coverage = 0.3, tx.prev.eff = 0.96, b.rate = 0.03/365,
-  b.rate.method = "totpop", b.propmale = NULL, ds.exit.age = 55,
-  ds.rate.mult = 1, di.cd4.aids = 50, di.cd4.rate = 2/365, ...)
+param_het(
+  time.unit = 7,
+  acute.stage.mult = 5,
+  aids.stage.mult = 1,
+  vl.acute.topeak = 14,
+  vl.acute.toset = 107,
+  vl.acute.peak = 6.7,
+  vl.setpoint = 4.5,
+  vl.aidsmax = 7,
+  cond.prob = 0.09,
+  cond.eff = 0.78,
+  act.rate.early = 0.362,
+  act.rate.late = 0.197,
+  act.rate.cd4 = 50,
+  acts.rand = TRUE,
+  circ.prob.birth = 0.9,
+  circ.eff = 0.53,
+  tx.elig.cd4 = 350,
+  tx.init.cd4.mean = 120,
+  tx.init.cd4.sd = 40,
+  tx.adhere.full = 0.76,
+  tx.adhere.part = 0.5,
+  tx.vlsupp.time = 365/3,
+  tx.vlsupp.level = 1.5,
+  tx.cd4.recrat.feml = 11.6/30,
+  tx.cd4.recrat.male = 9.75/30,
+  tx.cd4.decrat.feml = 11.6/30,
+  tx.cd4.decrat.male = 9.75/30,
+  tx.coverage = 0.3,
+  tx.prev.eff = 0.96,
+  b.rate = 0.03/365,
+  b.rate.method = "totpop",
+  b.propmale = NULL,
+  ds.exit.age = 55,
+  ds.rate.mult = 1,
+  di.cd4.aids = 50,
+  di.cd4.rate = 2/365,
+  ...
+)
 }
 \arguments{
 \item{time.unit}{Unit of time relative to one day.}
diff --git a/man/param_msm.Rd b/man/param_msm.Rd
index 557c271..200ca62 100644
--- a/man/param_msm.Rd
+++ b/man/param_msm.Rd
@@ -4,41 +4,91 @@
 \alias{param_msm}
 \title{Epidemic Model Parameters}
 \usage{
-param_msm(netstats, hiv.test.rate = c(0.01325, 0.0125, 0.0124),
-  hiv.test.late.prob = c(0.25, 0.25, 0.25), test.window.int = 21/7,
-  tt.part.supp = c(0.2, 0.2, 0.2), tt.full.supp = c(0.4, 0.4, 0.4),
-  tt.dur.supp = c(0.4, 0.4, 0.4), tx.init.prob = c(0.092, 0.092,
-  0.127), tx.halt.part.prob = c(0.0102, 0.0102, 0.0071),
-  tx.halt.full.rr = c(0.9, 0.9, 0.9), tx.halt.dur.rr = c(0.5, 0.5,
-  0.5), tx.reinit.part.prob = c(0.00066, 0.00066, 0.00291),
-  tx.reinit.full.rr = c(1, 1, 1), tx.reinit.dur.rr = c(1, 1, 1),
-  max.time.off.tx.full.int = 52 * 15, max.time.on.tx.part.int = 52 *
-  10, max.time.off.tx.part.int = 52 * 10, vl.acute.rise.int = 6.4,
-  vl.acute.peak = 6.886, vl.acute.fall.int = 6.4, vl.set.point = 4.5,
-  vl.aids.onset.int = 520, vl.aids.int = 104, vl.aids.peak = 7,
-  vl.full.supp = 1.5, vl.part.supp = 3.5, vl.tx.down.slope = 0.25,
-  vl.tx.aids.down.slope = 0.25, vl.tx.up.slope = 0.25,
-  aids.mr = 1/104, a.rate = 0.00052, arrival.age = 15,
-  URAI.prob = 0.008938, UIAI.prob = 0.003379, trans.scale = c(1, 1,
-  1), acute.rr = 6, circ.rr = 0.4, cond.eff = 0.95,
-  cond.fail = c(0.25, 0.25, 0.25), circ.prob = c(0.874, 0.874, 0.918),
-  epistats, acts.aids.vl = 5.75, acts.scale = 1, cond.scale = 1,
-  rgc.tprob = 0.35, ugc.tprob = 0.25, rct.tprob = 0.2,
-  uct.tprob = 0.16, rgc.sympt.prob = 0.16, ugc.sympt.prob = 0.8,
-  rct.sympt.prob = 0.14, uct.sympt.prob = 0.58, rgc.ntx.int = 16.8,
-  ugc.ntx.int = 16.8, gc.tx.int = 1.4, rct.ntx.int = 32,
-  uct.ntx.int = 32, ct.tx.int = 1.4, gc.sympt.prob.tx = c(0.95, 0.95,
-  0.95), ct.sympt.prob.tx = c(0.9, 0.9, 0.9),
-  gc.asympt.prob.tx = c(0.15, 0.15, 0.15), ct.asympt.prob.tx = c(0.15,
-  0.15, 0.15), sti.cond.eff = 0.9, sti.cond.fail = c(0.2, 0.2, 0.2),
-  hiv.rgc.rr = 2.78, hiv.ugc.rr = 1.73, hiv.rct.rr = 2.78,
-  hiv.uct.rr = 1.73, hiv.dual.rr = 0.2, riskh.start = Inf,
-  prep.start = Inf, prep.start.prob = 0.2, prep.adhr.dist = c(0.089,
-  0.127, 0.784), prep.adhr.hr = c(0.69, 0.19, 0.01),
-  prep.discont.rate = 1 - (2^(-1/(224.4237/7))), prep.tst.int = 90/7,
-  prep.risk.int = 182/7, prep.sti.screen.int = 182/7,
-  prep.sti.prob.tx = 1, prep.risk.reassess.method = "year",
-  prep.require.lnt = TRUE, ...)
+param_msm(
+  netstats,
+  hiv.test.rate = c(0.01325, 0.0125, 0.0124),
+  hiv.test.late.prob = c(0.25, 0.25, 0.25),
+  test.window.int = 21/7,
+  tt.part.supp = c(0.2, 0.2, 0.2),
+  tt.full.supp = c(0.4, 0.4, 0.4),
+  tt.dur.supp = c(0.4, 0.4, 0.4),
+  tx.init.prob = c(0.092, 0.092, 0.127),
+  tx.halt.part.prob = c(0.0102, 0.0102, 0.0071),
+  tx.halt.full.rr = c(0.9, 0.9, 0.9),
+  tx.halt.dur.rr = c(0.5, 0.5, 0.5),
+  tx.reinit.part.prob = c(0.00066, 0.00066, 0.00291),
+  tx.reinit.full.rr = c(1, 1, 1),
+  tx.reinit.dur.rr = c(1, 1, 1),
+  max.time.off.tx.full.int = 52 * 15,
+  max.time.on.tx.part.int = 52 * 10,
+  max.time.off.tx.part.int = 52 * 10,
+  vl.acute.rise.int = 6.4,
+  vl.acute.peak = 6.886,
+  vl.acute.fall.int = 6.4,
+  vl.set.point = 4.5,
+  vl.aids.onset.int = 520,
+  vl.aids.int = 104,
+  vl.aids.peak = 7,
+  vl.full.supp = 1.5,
+  vl.part.supp = 3.5,
+  vl.tx.down.slope = 0.25,
+  vl.tx.aids.down.slope = 0.25,
+  vl.tx.up.slope = 0.25,
+  aids.mr = 1/104,
+  a.rate = 0.00052,
+  arrival.age = 15,
+  URAI.prob = 0.008938,
+  UIAI.prob = 0.003379,
+  trans.scale = c(1, 1, 1),
+  acute.rr = 6,
+  circ.rr = 0.4,
+  cond.eff = 0.95,
+  cond.fail = c(0.25, 0.25, 0.25),
+  circ.prob = c(0.874, 0.874, 0.918),
+  epistats,
+  acts.aids.vl = 5.75,
+  acts.scale.main = 1,
+  acts.scale.casl = 1,
+  cond.scale = 1,
+  rgc.tprob = 0.35,
+  ugc.tprob = 0.25,
+  rct.tprob = 0.2,
+  uct.tprob = 0.16,
+  rgc.sympt.prob = 0.16,
+  ugc.sympt.prob = 0.8,
+  rct.sympt.prob = 0.14,
+  uct.sympt.prob = 0.58,
+  rgc.ntx.int = 16.8,
+  ugc.ntx.int = 16.8,
+  gc.tx.int = 1.4,
+  rct.ntx.int = 32,
+  uct.ntx.int = 32,
+  ct.tx.int = 1.4,
+  gc.sympt.prob.tx = c(0.95, 0.95, 0.95),
+  ct.sympt.prob.tx = c(0.9, 0.9, 0.9),
+  gc.asympt.prob.tx = c(0.15, 0.15, 0.15),
+  ct.asympt.prob.tx = c(0.15, 0.15, 0.15),
+  sti.cond.eff = 0.9,
+  sti.cond.fail = c(0.2, 0.2, 0.2),
+  hiv.rgc.rr = 2.78,
+  hiv.ugc.rr = 1.73,
+  hiv.rct.rr = 2.78,
+  hiv.uct.rr = 1.73,
+  hiv.dual.rr = 0.2,
+  riskh.start = Inf,
+  prep.start = Inf,
+  prep.start.prob = 0.2,
+  prep.adhr.dist = c(0.089, 0.127, 0.784),
+  prep.adhr.hr = c(0.69, 0.19, 0.01),
+  prep.discont.rate = 1 - (2^(-1/(224.4237/7))),
+  prep.tst.int = 90/7,
+  prep.risk.int = 182/7,
+  prep.sti.screen.int = 182/7,
+  prep.sti.prob.tx = 1,
+  prep.risk.reassess.method = "year",
+  prep.require.lnt = TRUE,
+  ...
+)
 }
 \arguments{
 \item{netstats}{Target statistics and related network initialization data from
@@ -164,7 +214,9 @@ will be circumcised (vector of length 3).}
 
 \item{acts.aids.vl}{Viral load level after which sexual act rate goes to zero.}
 
-\item{acts.scale}{Scalar for main/casual act rate for model calibration.}
+\item{acts.scale.main}{Scalar for main act rate for model calibration.}
+
+\item{acts.scale.casl}{Scalar for casual act rate for model calibration.}
 
 \item{cond.scale}{Scalar for condom use probability for model calibration.}
 
diff --git a/man/param_updater.Rd b/man/param_updater.Rd
new file mode 100644
index 0000000..f0d51d0
--- /dev/null
+++ b/man/param_updater.Rd
@@ -0,0 +1,49 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mod.dat_updater.R
+\name{param_updater}
+\alias{param_updater}
+\title{Module to modify the `param` list at specified time steps during the simulation}
+\usage{
+param_updater(dat, at)
+}
+\arguments{
+\item{dat}{Master data list object of class \code{dat} containing networks,
+individual-level attributes, and summary statistics.}
+
+\item{at}{Current time step.}
+}
+\description{
+Module to modify the `param` list at specified time steps during the simulation
+}
+\details{
+if a list `dat$param$param_updaters` is present, this module will update the
+`param` list with new values at given timesteps.
+`dat$control$param_updaters` is a list containing `updaters`. An updater is a
+list containing an `at` element telling when the changes will happend, an
+optional `verbose` boolean controlling whether to output a message when a
+change is made (default = TRUE) and a `param` list similar
+to the simulation's `dat$param` list.
+if the updater is a function but not the value to replace, the
+function will be applied to the current element (see example) .
+}
+\examples{
+ ## example of a `param_updaters` list
+ list(
+   list(
+     at = 4860,
+     param = list(
+       hiv.test.rate = rep(0.0128, 3),
+       trans.scale = c(1.61, 0.836, 0.622)
+     )
+   ),
+   list(
+     at = 5160,
+     verbose = FALSE,
+     param = list(
+       hiv.test.rate = function(x) x * 3,
+       trans.scale = function(x) x^2 / 3
+     )
+   )
+ )
+
+}
diff --git a/man/update_list.Rd b/man/update_list.Rd
new file mode 100644
index 0000000..678bae6
--- /dev/null
+++ b/man/update_list.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mod.dat_updater.R
+\name{update_list}
+\alias{update_list}
+\title{Update list `x` using the elements of list `new_x`}
+\usage{
+update_list(x, new_x)
+}
+\arguments{
+\item{x}{a list}
+
+\item{new_x}{a list}
+}
+\value{
+the full `x` list with the modifications added by `new_x`
+}
+\description{
+Update list `x` using the elements of list `new_x`
+}
+\details{
+This function updates list `x` by name. If `x` and `new_x` elements are not
+named, the function will not work properly.
+If a function is provided to replace an element that was originaly not a
+function, this function will be applied to the original value.
+}