flowchart LR
S["<b>S</b><br/>Susceptible"] -->|"infection<br/>(si.flow)"| I
I["<b>I</b><br/>Infectious"] -->|"natural clearance<br/>(is.flow.natural)"| S
I -->|"screen + treat index<br/>(n.index.tx)"| S
I -->|"partner notification<br/>(n.partner.cleared.pn)"| S
style S fill:#3498db,color:#fff
style I fill:#e74c3c,color:#fff
Partner Notification for an Endemic STI
SIS
STI
intervention
partner notification
cumulative edgelist
intermediate
An SIS endemic-STI model with partner notification. Compares no PN to Patient Referral and Expedited Partner Therapy, varying trace rate and lookback window. Showcases EpiModel’s cumulative-edgelist API.
Overview
This example builds a partner notification (PN) intervention for an endemic bacterial STI, modeled as an SIS process on a dynamic sexual contact network. Screening identifies infected indices, indices are treated, and a custom partner_services module looks up each index’s recent partners from the cumulative edgelist and routes them through one of two arms: Patient Referral (PR), where the partner returns for testing and treatment if positive, or Expedited Partner Therapy (EPT), where the partner is dispensed medication directly without testing. Five scenarios share the same network and the same disease parameters and differ only in the PN configuration.
The pedagogical core is the cumulative-edgelist API. EpiModel stores every partnership the engine has seen, optionally truncated to a lookback window, and exposes a small set of accessors for retrieving the past contacts of a given node: get_partners(), get_cumulative_edgelist(), and the unique-vs-positional-ID conversion helpers get_posit_ids() and get_unique_ids(). Once a custom module can ask “who were this person’s partners in the last 60 weeks?” any partner-based intervention becomes a matter of choosing what to do with the returned set.
TipDownload standalone scripts
- model.R - Main simulation script
- module-fx.R - Custom module functions
Model Structure
| Compartment | Label | Description |
|---|---|---|
| Susceptible | S | Not infected; at risk and reinfectable |
| Infectious | I | Infected and transmitting |
Three return paths from I to S: slow natural clearance, treatment of screened-positive indices, and treatment of notified partners (PR or EPT). Reinfection is allowed throughout, which motivates the SIS framing: the post-PN equilibrium prevalence is the quantity of interest, not the time to elimination.
The Cumulative Edgelist Pattern
The cumulative-edgelist machinery lives in three places.
Three flags on
control.net():control.net( ..., cumulative.edgelist = TRUE, # turn the engine on truncate.el.cuml = pn.lookback, # drop edges older than this save.cumulative.edgelist = TRUE # attach to returned sim )truncate.el.cumlis destructive: edges older than the window are dropped from storage and cannot be recovered. Set it to the maximum lookback across all scenarios so each scenario has the full edge history it needs.Inside the
partner_servicesmodule, ask for each index’s recent partners:idsIndex <- which(active == 1 & dx.this.step == 1) part_df <- get_partners(dat, idsIndex, truncate = pn.lookback, only.active.nodes = TRUE)The unique-vs-positional ID round-trip.
get_partners()returns unique IDs in thepartnercolumn, because past partners may have already departed the population and no longer have a positional ID. Per-node attributes (status,active, etc.) are indexed by positional ID. Convert before indexing:partner_pid <- get_posit_ids(dat, part_df$partner)only.active.nodes = TRUEguarantees the returned partners are still active, so the conversion will succeed, but the conversion is still required. AnyNAs that slip through (e.g., a partner who departs between the cumulative-edgelist update and the conversion call) should be dropped.
Treat the unique-vs-positional distinction as a hard rule of the cumulative-edgelist API: values returned by get_partners() are in unique-ID space and must be converted before any attribute indexing.
Setup
suppressMessages(library(EpiModel))
nsims <- 5
ncores <- 5
nsteps <- 600
n <- 1000
pn.start <- 300Custom Modules
Screening Module
Routine screening of infected actives at rate screen.rate per week. Positives set the per-step dx.this.step flag (the PN trigger) and stamp dx.time; diag.status is sticky and remains 1 even after later clearance, marking the node as previously diagnosed.
screen <- function(dat, at) {
active <- get_attr(dat, "active")
status <- get_attr(dat, "status")
1 if (at == 2) {
n <- length(active)
dat <- set_attr(dat, "diag.status", rep(0, n))
dat <- set_attr(dat, "dx.time", rep(NA_integer_, n))
dat <- set_attr(dat, "dx.this.step", rep(0, n))
dat <- set_attr(dat, "tx.this.step", rep(0, n))
dat <- set_attr(dat, "pn.notified", rep(NA_integer_, n))
dat <- set_attr(dat, "infections", ifelse(status == "i", 1, 0))
}
diag.status <- get_attr(dat, "diag.status")
dx.time <- get_attr(dat, "dx.time")
# Reset per-step flags before this step's events fire.
dx.this.step <- rep(0L, length(active))
tx.this.step <- rep(0L, length(active))
screen.rate <- get_param(dat, "screen.rate")
pn.start <- get_param(dat, "pn.start")
nDxNew <- 0
if (screen.rate > 0) {
idsElig <- which(active == 1 & status == "i")
if (length(idsElig) > 0) {
vec <- which(rbinom(length(idsElig), 1, screen.rate) == 1)
if (length(vec) > 0) {
idsDx <- idsElig[vec]
nDxNew <- length(idsDx)
2 dx.this.step[idsDx] <- 1
diag.status[idsDx] <- 1
dx.time[idsDx] <- at
}
}
}
dat <- set_attr(dat, "diag.status", diag.status)
dat <- set_attr(dat, "dx.time", dx.time)
dat <- set_attr(dat, "dx.this.step", dx.this.step)
dat <- set_attr(dat, "tx.this.step", tx.this.step)
dat <- set_epi(dat, "n.indices", at, nDxNew)
dat <- set_epi(dat, "pn.on", at, as.numeric(at >= pn.start))
dat <- set_epi(dat, "n.diag.ever", at,
sum(active == 1 & diag.status == 1, na.rm = TRUE))
return(dat)
}- 1
-
Lazy initialization at
at == 2. EpiModel reservesat == 1forinit.net()setup. - 2
-
dx.this.stepis the canonical PN trigger. It is reset to0at the top of the screening module each step and set to1on fresh-positive screens; downstream modules read it within the same step.
Partner Notification Module
The partner_services module is kept intentionally minimal. It identifies fresh-positive indices, queries the cumulative edgelist for their recent partners, converts unique IDs to positional IDs, applies a Bernoulli trace, and stamps the reached partners with pn.notified = at. The downstream treat module reads that stamp to determine how to handle each notified partner under the active PN arm.
partner_services <- function(dat, at) {
active <- get_attr(dat, "active")
status <- get_attr(dat, "status")
dx.this.step <- get_attr(dat, "dx.this.step")
pn.notified <- get_attr(dat, "pn.notified")
pn.lookback <- get_param(dat, "pn.lookback")
pn.trace.prob <- get_param(dat, "pn.trace.prob")
pn.start <- get_param(dat, "pn.start")
n.partners.elig <- 0
n.partners.reached <- 0
if (at >= pn.start && pn.trace.prob > 0) {
1 idsIndex <- which(active == 1 & dx.this.step == 1)
if (length(idsIndex) > 0) {
2 part_df <- get_partners(
dat, idsIndex,
truncate = pn.lookback,
only.active.nodes = TRUE
)
if (!is.null(part_df) && nrow(part_df) > 0) {
3 partner_pid <- get_posit_ids(dat, part_df$partner)
partner_pid <- partner_pid[!is.na(partner_pid)]
partner_pid <- unique(partner_pid)
4 partner_pid <- setdiff(partner_pid, idsIndex)
n.partners.elig <- length(partner_pid)
if (n.partners.elig > 0) {
reached <- partner_pid[
5 rbinom(n.partners.elig, 1, pn.trace.prob) == 1
]
n.partners.reached <- length(reached)
if (n.partners.reached > 0) {
pn.notified[reached] <- at
}
}
}
}
}
dat <- set_attr(dat, "pn.notified", pn.notified)
dat <- set_epi(dat, "n.partners.elig", at, n.partners.elig)
dat <- set_epi(dat, "n.partners.reached", at, n.partners.reached)
return(dat)
}- 1
-
Indices are this step’s fresh-positive screens. They are guaranteed positional IDs (we are reading them out of
active/dx.this.step), so we can pass them straight intoget_partners(). - 2
-
get_partners()returns one row per (index, partner) edge that falls inside the lookback window.only.active.nodes = TRUEdrops partners who have departed the population since the partnership. - 3
-
ID space conversion.
part_df$partneris in unique-ID space, while every per-node vector indatis in positional-ID space. Conversion must occur before indexing. Without this round-trip the next line would silently corrupt the state of the wrong nodes. - 4
-
A partner could also be a fresh index this step. Such nodes are handled by the index-treatment path in
treat(); double-counting them as notified partners would inflate the cascade statistics. - 5
- The Bernoulli trace probability summarizes the “successfully notified” rate. In practice this is the product of contact-information availability, contact-tracer effort, and partner cooperation; here it is a single tunable parameter.
Treatment Module
The treatment cascade has two pathways. Indices are always offered treatment at tx.efficacy. Notified partners are handled per pn.arm:
"none": nothing. (Settingpn.trace.prob = 0upstream means no notifications are emitted; this branch never fires.)"PR"(Patient Referral): test the partner first at sensitivitypn.test.prob, then treat if positive attx.efficacy. Uninfected partners are not treated. Effective coverage on infected partners ispn.test.prob * tx.efficacy."EPT"(Expedited Partner Therapy): dispense medication directly to the partner atept.efficacy. Treats both infected and uninfected partners; uninfected treatment has no epidemic effect but is counted as a “wasted dose” in the cascade analysis (the dispensed-but-not-needed cost).
treat <- function(dat, at) {
active <- get_attr(dat, "active")
status <- get_attr(dat, "status")
dx.this.step <- get_attr(dat, "dx.this.step")
tx.this.step <- get_attr(dat, "tx.this.step")
pn.notified <- get_attr(dat, "pn.notified")
diag.status <- get_attr(dat, "diag.status")
tx.efficacy <- get_param(dat, "tx.efficacy")
ept.efficacy <- get_param(dat, "ept.efficacy")
pn.test.prob <- get_param(dat, "pn.test.prob")
pn.arm <- get_param(dat, "pn.arm")
n.index.tx <- 0
n.partner.tx <- 0
n.partner.tx.wasted <- 0
n.partner.cleared.pn <- 0
# Index treatment
idsIndex <- which(active == 1 & dx.this.step == 1)
if (length(idsIndex) > 0) {
vec <- which(rbinom(length(idsIndex), 1, tx.efficacy) == 1)
if (length(vec) > 0) {
idsClear <- idsIndex[vec]
idsClear <- idsClear[status[idsClear] == "i"]
n.index.tx <- length(idsClear)
if (n.index.tx > 0) {
status[idsClear] <- "s"
tx.this.step[idsClear] <- 1
}
}
}
# Notified-partner treatment
idsPart <- which(active == 1 & pn.notified == at)
if (length(idsPart) > 0 && pn.arm %in% c("PR", "EPT")) {
1 if (pn.arm == "PR") {
idsPartI <- idsPart[status[idsPart] == "i"]
if (length(idsPartI) > 0) {
testPos <- which(rbinom(length(idsPartI), 1, pn.test.prob) == 1)
idsTestPos <- idsPartI[testPos]
if (length(idsTestPos) > 0) {
rxOk <- which(rbinom(length(idsTestPos), 1, tx.efficacy) == 1)
idsRx <- idsTestPos[rxOk]
n.partner.tx <- length(idsRx)
if (n.partner.tx > 0) {
status[idsRx] <- "s"
diag.status[idsRx] <- 1
tx.this.step[idsRx] <- 1
n.partner.cleared.pn <- n.partner.tx
}
}
}
2 } else if (pn.arm == "EPT") {
rxOk <- which(rbinom(length(idsPart), 1, ept.efficacy) == 1)
idsRx <- idsPart[rxOk]
n.partner.tx <- length(idsRx)
if (n.partner.tx > 0) {
tx.this.step[idsRx] <- 1
isI <- status[idsRx] == "i"
n.partner.tx.wasted <- sum(!isI)
idsClear <- idsRx[isI]
if (length(idsClear) > 0) {
status[idsClear] <- "s"
n.partner.cleared.pn <- length(idsClear)
}
}
}
}
dat <- set_attr(dat, "status", status)
dat <- set_attr(dat, "tx.this.step", tx.this.step)
dat <- set_attr(dat, "diag.status", diag.status)
dat <- set_epi(dat, "n.index.tx", at, n.index.tx)
dat <- set_epi(dat, "n.partner.tx", at, n.partner.tx)
dat <- set_epi(dat, "n.partner.tx.wasted", at, n.partner.tx.wasted)
dat <- set_epi(dat, "n.partner.cleared.pn", at, n.partner.cleared.pn)
return(dat)
}- 1
- Patient Referral: filter to infected partners, then test, then treat. Uninfected partners exit the cascade at the testing gate.
- 2
- Expedited Partner Therapy: no filter, no test. Every notified partner is dispensed. Wasted doses on uninfected partners are recorded but have no epidemic effect.
Recovery and Infection Modules
Recovery is a slow natural clearance at rec.rate. Infection is the standard discordant-edge transmission, augmented to increment a per-node infections counter on every fresh S to I (used downstream for reinfection analysis).
recov <- function(dat, at) {
active <- get_attr(dat, "active")
status <- get_attr(dat, "status")
rec.rate <- get_param(dat, "rec.rate")
nRec <- 0
idsElig <- which(active == 1 & status == "i")
if (length(idsElig) > 0) {
vec <- which(rbinom(length(idsElig), 1, rec.rate) == 1)
if (length(vec) > 0) {
idsRec <- idsElig[vec]
nRec <- length(idsRec)
status[idsRec] <- "s"
}
}
dat <- set_attr(dat, "status", status)
dat <- set_epi(dat, "is.flow.natural", at, nRec)
return(dat)
}infect <- function(dat, at) {
active <- get_attr(dat, "active")
status <- get_attr(dat, "status")
infTime <- get_attr(dat, "infTime")
infections <- get_attr(dat, "infections")
inf.prob <- get_param(dat, "inf.prob")
act.rate <- get_param(dat, "act.rate")
idsInf <- which(active == 1 & status == "i")
nActive <- sum(active == 1)
nElig <- length(idsInf)
nInf <- 0
if (nElig > 0 && nElig < nActive) {
del <- discord_edgelist(dat, at)
if (!is.null(del) && nrow(del) > 0) {
del$transProb <- inf.prob
del$actRate <- act.rate
del$finalProb <- 1 - (1 - del$transProb)^del$actRate
transmit <- rbinom(nrow(del), 1, del$finalProb)
del <- del[which(transmit == 1), ]
idsNewInf <- unique(del$sus)
nInf <- length(idsNewInf)
if (nInf > 0) {
status[idsNewInf] <- "i"
infTime[idsNewInf] <- at
infections[idsNewInf] <- infections[idsNewInf] + 1
}
}
}
if (nInf > 0) {
dat <- set_attr(dat, "status", status)
dat <- set_attr(dat, "infTime", infTime)
dat <- set_attr(dat, "infections", infections)
}
dat <- set_epi(dat, "si.flow", at, nInf)
return(dat)
}Network Model
A sparse sexual contact network with some concurrency. Mean degree 0.7 and roughly 7% concurrency are in the empirically observed range for sexual partnership networks. Each timestep is one week.
nw <- network_initialize(n)
formation <- ~edges + concurrent
target.stats <- c(round(0.7 * n / 2), round(0.07 * n))
coef.diss <- dissolution_coefs(~offset(edges), duration = 100)
est <- netest(nw, formation, target.stats, coef.diss)Network Diagnostics
dx <- netdx(est, nsims = nsims, ncores = ncores, nsteps = nsteps,
nwstats.formula = ~edges + concurrent + degree(0:3))
Network Diagnostics
-----------------------
- Simulating 5 networks
- Calculating formation statisticsprint(dx)EpiModel Network Diagnostics
=======================
Diagnostic Method: Dynamic
Simulations: 5
Time Steps per Sim: 600
Formation Diagnostics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 350 344.237 -1.647 2.275 -2.533 7.267 13.775
concurrent 70 67.821 -3.112 1.383 -1.575 3.927 7.987
degree0 NA 415.040 NA 2.318 NA 9.526 18.190
degree1 NA 517.139 NA 2.070 NA 6.336 16.420
degree2 NA 42.122 NA 0.641 NA 2.553 5.868
degree3 NA 17.958 NA 0.441 NA 1.528 4.015
Duration Diagnostics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 100 99.231 -0.769 0.805 -0.955 0.933 4.328
Dissolution Diagnostics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 0.01 0.01 0.463 0 0.488 0 0.005plot(dx)
Epidemic Simulation
Control Settings
The cumulative-edgelist machinery is enabled here. truncate.el.cuml is set to the maximum lookback across scenarios so the edge history covers every scenario’s window.
init <- init.net(i.num = round(0.10 * n))
max.lookback <- 120
control <- control.net(
type = NULL,
nsims = nsims,
ncores = ncores,
nsteps = nsteps,
infection.FUN = infect,
screen.FUN = screen,
partner_services.FUN = partner_services,
treat.FUN = treat,
recovery.FUN = recov,
1 cumulative.edgelist = TRUE,
2 truncate.el.cuml = max.lookback,
3 save.cumulative.edgelist = TRUE,
verbose = FALSE
)- 1
-
Turn on cumulative-edgelist tracking. Required for
get_partners()to return anything. - 2
-
Drop edges older than
max.lookbackweeks. This is destructive: there is no recovering edges that were dropped earlier. - 3
-
Save the final cumulative edgelist on the returned
netsimobject so it can be inspected post-hoc.
Base Parameters
# Base parameter set holds defaults for every parameter the modules
# read. PN is off by default (pn.arm = "none"); per-scenario overrides
# are applied via the scenarios API in the next section.
param_base <- param.net(
inf.prob = 0.35,
act.rate = 1,
rec.rate = 0.01,
screen.rate = 0.015,
tx.efficacy = 0.95,
ept.efficacy = 0.85,
pn.test.prob = 0.85,
pn.arm = "none",
pn.trace.prob = 0,
pn.lookback = 60,
pn.start = pn.start
)Scenarios
All five scenarios share the same burn-in (steps 1 to pn.start = 300, PN off). PN switches on at step 300 and runs to step 600.
scenarios.df <- data.frame(
.scenario.id = c("none", "pr", "ept", "ept_long", "ept_max"),
.at = 0,
pn.arm = c("none", "PR", "EPT", "EPT", "EPT"),
pn.trace.prob = c(0, 0.5, 0.5, 0.5, 0.8),
pn.lookback = c(60, 60, 60, 120, 120),
stringsAsFactors = FALSE
)
scenarios.list <- create_scenario_list(scenarios.df)
sims <- list()
for (scn in scenarios.list) {
sims[[scn$id]] <- netsim(est, use_scenario(param_base, scn),
init, control)
}
sim_none <- sims[["none"]]
sim_pr <- sims[["pr"]]
sim_ept <- sims[["ept"]]
sim_ept_long <- sims[["ept_long"]]
sim_ept_max <- sims[["ept_max"]]The progression across scenarios isolates each intervention component:
- Scenario 1 establishes endemic prevalence under screening alone (the counterfactual).
- Scenario 2 adds PR, the standard public-health intervention.
- Scenario 3 swaps PR for EPT at the same trace rate, isolating the gain from bypassing the test gate.
- Scenario 4 extends the lookback window from 60 to 120 weeks, illustrating diminishing returns at the tail of the edge-age distribution.
- Scenario 5 combines both levers (high trace and long lookback). The improvement over Scenario 4 should be visible but moderate.
Note that 60 weeks (~14 months) is much longer than the CDC chlamydia/gonorrhea recommendation of 60 days. The partnership durations and act rates here are tuned for clean endemic-equilibrium teaching, not for calibrating to U.S. chlamydia surveillance data; treat the numbers as illustrative.
Analysis
sims <- list(
none = sim_none,
pr = sim_pr,
ept = sim_ept,
ept_long = sim_ept_long,
ept_max = sim_ept_max
)
labels <- c(none = "Screening only",
pr = "PR, 50%, lookback 60",
ept = "EPT, 50%, lookback 60",
ept_long = "EPT, 50%, lookback 120",
ept_max = "EPT, 80%, lookback 120")
cols <- c(none = "gray40",
pr = "#8e44ad",
ept = "steelblue",
ept_long = "#16a085",
ept_max = "firebrick")
sims <- lapply(sims, function(s) mutate_epi(s, prev = i.num / num))Prevalence Trajectory
The shared burn-in (steps 1 to 300) shows all five scenarios converging to the same endemic prevalence under screening alone. At step 300 PN switches on (vertical dashed line) and the scenarios bifurcate as each settles into a new lower equilibrium.
par(mfrow = c(1, 1), mar = c(3, 3, 2, 1), mgp = c(2, 1, 0))
plot(sims$none, y = "prev",
main = "Prevalence by Partner Notification Scenario",
ylab = "Prevalence (I / N)", xlab = "Week",
mean.col = cols["none"], mean.lwd = 2, mean.smooth = TRUE,
qnts = FALSE, legend = FALSE,
ylim = c(0, 1))
for (k in c("pr", "ept", "ept_long", "ept_max")) {
plot(sims[[k]], y = "prev",
mean.col = cols[k], mean.lwd = 2, mean.smooth = TRUE,
qnts = FALSE, legend = FALSE, add = TRUE)
}
abline(v = pn.start, lty = 2, col = "gray40")
legend("topright", legend = labels, col = cols, lwd = 2, bty = "n",
cex = 0.75)
Reinfection Pressure
For each scenario we compute the total subsequent infections divided by the total indices detected in the post-PN window. This ratio captures the per-detected-case transmission burden remaining in the population. Under aggressive PN the absolute number of new infections drops, but so does the number of indices detected; the ratio can rise when index detection falls faster than incident transmission. The plot shows the per-simulation values and the across-simulation mean.
post_window <- function(df) df$time > pn.start
reinf_rate <- function(s) {
df <- as.data.frame(s)
df <- df[post_window(df), ]
n_subsequent <- tapply(df$si.flow, df$sim, sum, na.rm = TRUE)
n_indices <- tapply(df$n.indices, df$sim, sum, na.rm = TRUE)
ifelse(n_indices > 0, n_subsequent / n_indices, NA)
}
ri_by_scen <- lapply(sims, reinf_rate)
ri_mean <- sapply(ri_by_scen, mean, na.rm = TRUE)
xrange <- range(unlist(ri_by_scen), na.rm = TRUE)
par(mfrow = c(1, 1), mar = c(4, 12, 2, 1), mgp = c(2.5, 1, 0))
plot(NA, xlim = c(0, xrange[2] * 1.05), ylim = c(0.5, length(ri_by_scen) + 0.5),
yaxt = "n", xlab = "Subsequent infections per index (post-PN)",
ylab = "",
main = "Reinfection Pressure by Scenario")
axis(2, at = seq_along(ri_by_scen), labels = labels[names(ri_by_scen)],
las = 1, cex.axis = 0.85)
abline(h = seq_along(ri_by_scen), col = "gray90", lty = 1)
for (i in seq_along(ri_by_scen)) {
vals <- ri_by_scen[[i]]
points(vals, rep(i, length(vals)),
pch = 21, bg = adjustcolor(cols[names(ri_by_scen)[i]], 0.6),
col = cols[names(ri_by_scen)[i]], cex = 1.2)
segments(ri_mean[i], i - 0.25, ri_mean[i], i + 0.25,
lwd = 3, col = cols[names(ri_by_scen)[i]])
}
Tracing Cascade
Four bars per scenario, showing where the funnel narrows: (a) total indices detected, (b) partners notified, (c) partners actually treated, (d) partners cleared (infected and treated successfully) by PN.
cascade_summary <- function(s) {
df <- as.data.frame(s)
df <- df[post_window(df), ]
c(indices = mean(tapply(df$n.indices, df$sim, sum, na.rm = TRUE)),
notified = mean(tapply(df$n.partners.reached, df$sim, sum, na.rm = TRUE)),
treated = mean(tapply(df$n.partner.tx, df$sim, sum, na.rm = TRUE)),
cleared_pn = mean(tapply(df$n.partner.cleared.pn, df$sim, sum, na.rm = TRUE)),
wasted = mean(tapply(df$n.partner.tx.wasted, df$sim, sum, na.rm = TRUE)),
nat_clear = mean(tapply(df$is.flow.natural, df$sim, sum, na.rm = TRUE)),
index_clear = mean(tapply(df$n.index.tx, df$sim, sum, na.rm = TRUE)),
eq_prev = mean(df$prev[df$time > nsteps - 100], na.rm = TRUE))
}
casc <- sapply(sims, cascade_summary)
cascade_mat <- casc[c("indices", "notified", "treated", "cleared_pn"), ]
cascade_labels <- c(none = "Screening\nonly",
pr = "PR 50%\nlb 60",
ept = "EPT 50%\nlb 60",
ept_long = "EPT 50%\nlb 120",
ept_max = "EPT 80%\nlb 120")
colnames(cascade_mat) <- cascade_labels[colnames(cascade_mat)]
par(mfrow = c(1, 1), mar = c(4, 5, 4, 1), mgp = c(3.5, 1, 0))
bp <- barplot(cascade_mat, beside = TRUE,
col = c("#7f8c8d", "#3498db", "#f39c12", "#27ae60"),
names.arg = colnames(cascade_mat),
las = 1, cex.names = 0.85,
ylab = "Mean count over post-PN window",
main = "Partner-Notification Cascade",
ylim = c(0, max(cascade_mat, na.rm = TRUE) * 1.25))
legend("top",
legend = c("Indices", "Notified", "Treated", "Cleared"),
fill = c("#7f8c8d", "#3498db", "#f39c12", "#27ae60"),
bty = "n", cex = 0.85, horiz = TRUE, inset = c(0, -0.08),
xpd = TRUE)
Summary Table
totals <- as.data.frame(t(casc))
totals$Scenario <- labels
totals$Per_index_notified <- ifelse(totals$indices > 0,
totals$notified / totals$indices, NA)
totals$Per_notified_treated <- ifelse(totals$notified > 0,
totals$treated / totals$notified, NA)
totals$PN_share_of_clearance <- with(totals,
ifelse((cleared_pn + index_clear +
nat_clear) > 0,
cleared_pn /
(cleared_pn + index_clear +
nat_clear),
NA))
totals$Reinf_per_index <- sapply(ri_by_scen, function(x) mean(x, na.rm = TRUE))
knitr::kable(data.frame(
Scenario = totals$Scenario,
EqPrev = round(totals$eq_prev, 3),
Indices = round(totals$indices),
Notified_per_idx = round(totals$Per_index_notified, 2),
Treated_per_ntfd = round(totals$Per_notified_treated, 2),
PN_share_clear = round(totals$PN_share_of_clearance, 2),
Reinf_per_idx = round(totals$Reinf_per_index, 2),
row.names = NULL
))| Scenario | EqPrev | Indices | Notified_per_idx | Treated_per_ntfd | PN_share_clear | Reinf_per_idx |
|---|---|---|---|---|---|---|
| Screening only | 0.214 | 974 | 0.00 | NA | 0.00 | 1.63 |
| PR, 50%, lookback 60 | 0.116 | 625 | 0.86 | 0.65 | 0.26 | 2.03 |
| EPT, 50%, lookback 60 | 0.087 | 483 | 0.87 | 0.84 | 0.27 | 2.01 |
| EPT, 50%, lookback 120 | 0.114 | 666 | 1.16 | 0.85 | 0.30 | 2.12 |
| EPT, 80%, lookback 120 | 0.079 | 481 | 1.80 | 0.84 | 0.39 | 2.35 |
PN_share_clear is the fraction of all clearance events in the post-PN window attributable to partner-treatment. It compares the PN intervention’s contribution against natural clearance and index treatment combined: a value of 0.30 means roughly one in three clearance events was caused by partner notification.
Parameters
Transmission and Recovery
| Parameter | Description | Default |
|---|---|---|
inf.prob |
Per-act transmission probability | 0.35 |
act.rate |
Acts per partnership per week | 1 |
rec.rate |
Weekly natural clearance rate | 0.01 (mean ~100 wk) |
Screening and Treatment
| Parameter | Description | Default |
|---|---|---|
screen.rate |
Weekly probability of screening for an infected | 0.015 (~67 wk between screens) |
tx.efficacy |
Probability a treated index actually clears | 0.95 |
ept.efficacy |
Probability a partner takes EPT meds and clears | 0.85 |
pn.test.prob |
Sensitivity of the returning-partner test in PR | 0.85 |
Partner Notification
| Parameter | Description | Varied across scenarios |
|---|---|---|
pn.arm |
"none", "PR", or "EPT" |
Yes |
pn.trace.prob |
Bernoulli probability a partner is reached | 0, 0.5, 0.8 |
pn.lookback |
Lookback window (weeks) on the cumulative edgelist | 60 or 120 |
pn.start |
Step at which PN switches on (after burn-in) | 300 |
Network
| Parameter | Description | Default |
|---|---|---|
| Population size | Number of nodes | 1000 |
| Mean degree | Edges per node | 0.7 |
| Concurrency | Nodes with degree > 1 | ~7% |
| Partnership duration | Mean edge duration (weeks) | 100 |
Module Execution Order
resim_nets -> summary_nets -> infection -> screen -> partner_services ->
treat -> recovery -> nwupdate -> prevalencescreen runs first among the custom modules so that this step’s fresh indices (dx.this.step == 1) are visible to partner_services and treat. partner_services then asks the cumulative edgelist for those indices’ partners and stamps notifications. treat reads both dx.this.step (indices) and pn.notified == at (notified partners) and applies the arm-specific treatment cascade. recovery handles slow natural clearance for everyone still in I.
Caveats
pn.lookback = 60weeks is much longer than the real CDC guidance (60 days for chlamydia/gonorrhea). Partnership durations and act rates here are tuned for clean endemic-equilibrium teaching, not for calibration to chlamydia surveillance data. Treat all numbers as illustrative.- The natural clearance rate is geometric with mean ~50 weeks. Real chlamydia clears more slowly and with higher between-host variability.
- Partner notification is modeled as instantaneous within the timestep. Real systems have multi-week delays between index diagnosis, contact-tracer outreach, and partner uptake.
- EPT is permissive in this model: any notified partner accepts dispensed medication at
ept.efficacy. Real EPT acceptance varies by partnership type, region, and pharmacy availability.
Next Steps
- Add a clinic-visit delay between index diagnosis and partner outreach to study the cost of delayed PN.
- Stratify the trace rate by partnership type (main vs. casual) using a multilayer network and the
networksargument ofget_partners(). - Allow re-notification of partners who were notified earlier but did not get treated, by relaxing the
pn.notified == atfilter intreat()to a window check. - Replace EPT efficacy with a stochastic uptake decision that depends on partnership recency, modeled from the
startcolumn of the cumulative edgelist. - Pair with contact tracing for a respiratory pathogen example, sharing the cumulative-edgelist API but with different triggers and a different downstream intervention.