m <- mod

if (!is.null(init)) m <- update(m, init = init)

out <- m %>% param(...) %>%

mrgsim(end = hours, delta = hours, recsort = 3) %>% as.data.frame()

tail(out, 1)

cmts <- c("S", "T4", "T3", "rT3", "Ggut", "Ipl", "Ithy", "NIS", "TSH",

"NEFA", "TBG", "SNS", "Rb", "HR", "CR", "Tc", "Vol", "CNSx",

"Bili", "Cort", "GIx", "AFb", "Hz", "BWi", "Prec", "PGE",

"Aptu", "Cptu", "Ammi", "Cmmi", "Apro", "Cpro", "Cesm", "Cgc",

"Achol", "Ciop", "Casa", "Capap")

as.list(row[1, cmts])

st <- ts_euthyroid()

ladder <- unique(c(seq(1.5, min(TRAb, 4.0), by = 0.5), TRAb))

for (tr in ladder) {

st <- ts_state(ts_settle(ts_mod, 4000, init = st, TRAb = tr))

}

st$Hz <- 0; st$BWi <- 0 # reset the outcome accumulators

st

`1 Untreated` = list(ev = ts_none(), p = list()),

`2 PTU alone` = list(ev = ev_ptu(), p = list()),

`3 Methimazole alone` = list(ev = ev_mmi(), p = list()),

`4 Iodide (1 h after PTU)` = list(ev = c(ev_ptu(), ev_iod(1)),

p = list()),

`5 Propranolol alone` = list(ev = ev_pro(), p = list()),

`6 Esmolol alone` = list(ev = ts_none(),

p = list(ESMRATE = 1.60)),

`7 Hydrocortisone alone` = list(ev = ev_hc(), p = list()),

`8 Cooling + fluids alone` = list(ev = ts_none(),

p = list(COOL = 1.2,

FLUID = 0.0045)),

`9 FULL BUNDLE` = list(

ev = c(ev_ptu(), ev_iod(1), ev_pro(), ev_hc()),

p = list(COOL = 1.2, FLUID = 0.0045)),

`10 BUNDLE with methimazole` = list(

ev = c(ev_mmi(), ev_iod(1), ev_pro(), ev_hc()),

p = list(COOL = 1.2, FLUID = 0.0045)),

`11 BUNDLE, iodide BEFORE thionamide` = list(

ev = c(ev_ptu(), ev_iod(0), ev_pro(), ev_hc()),

p = list(COOL = 1.2, FLUID = 0.0045)),

`12 BUNDLE minus beta-blockade` = list(

ev = c(ev_ptu(), ev_iod(1), ev_hc()),

p = list(COOL = 1.2, FLUID = 0.0045)),

`13 BUNDLE minus cooling/fluids` = list(

ev = c(ev_ptu(), ev_iod(1), ev_pro(), ev_hc()), p = list()),

`14 BUNDLE + aspirin` = list(

ev = c(ev_ptu(), ev_iod(1), ev_pro(), ev_hc(), ev_asa()),

p = list(COOL = 1.2, FLUID = 0.0045)),

`15 BUNDLE + acetaminophen` = list(

ev = c(ev_ptu(), ev_iod(1), ev_pro(), ev_hc(), ev_apap()),

p = list(COOL = 1.2, FLUID = 0.0045)),

`16 BUNDLE + iopanoic + cholestyramine` = list(

ev = c(ev_ptu(), ev_iod(1), ev_pro(), ev_hc(), ev_iop(), ev_chol()),

p = list(COOL = 1.2, FLUID = 0.0045)),

`17 BUNDLE + plasma exchange at 12 h` = list(

ev = c(ev_ptu(), ev_iod(1), ev_pro(), ev_hc()),

p = list(COOL = 1.2, FLUID = 0.0045, TPE_ON = 1)),

`18 BUNDLE with esmolol` = list(

ev = c(ev_ptu(), ev_iod(1), ev_hc()),

p = list(COOL = 1.2, FLUID = 0.0045, ESMRATE = 1.60))

prec = TS_PREC, pge = TS_PGE) {

if (is.null(base)) base <- ts_thyrotoxic(TRAb)

init <- base

init$Prec <- prec

init$PGE <- pge

m <- ts_mod %>% update(init = init) %>% param(TRAb = TRAb)

if (length(spec$p)) m <- do.call(param, c(list(m), spec$p))

m %>% mrgsim(events = spec$ev, end = hours, delta = delta, recsort = 3) %>%

as.data.frame()

base <- ts_thyrotoxic(TRAb)

sc <- ts_scenarios()

out <- lapply(names(sc), function(nm) {

d <- ts_run(sc[[nm]], base = base, TRAb = TRAb, hours = hours)

d$scenario <- nm

d

})

bind_rows(out)

res %>% group_by(scenario) %>%

summarise(

Tmax = max(Tc),

Tc_24h = Tc[which.min(abs(time - 24))],

BWPS_24h = BWPS[which.min(abs(time - 24))],

BWPS_72h = BWPS[which.min(abs(time - 72))],

totT3_24h = T3[which.min(abs(time - 24))],

freeT3_24 = fT3[which.min(abs(time - 24))],

HR_24h = HR[which.min(abs(time - 24))],

runaway = max(Tc) >= 40,

mort_pct = 100 * max(mort),

.groups = "drop") %>%

arrange(scenario) %>% as.data.frame()

hours = 72) {

st <- ts_euthyroid()

rows <- list()

for (tr in trabs) {

st <- ts_state(ts_settle(ts_mod, 4000, init = st, TRAb = tr))

s2 <- st; s2$Hz <- 0; s2$BWi <- 0

d <- ts_mod %>% update(init = s2) %>% param(TRAb = tr) %>%

mrgsim(end = hours, delta = 0.25) %>% as.data.frame()

rows[[length(rows) + 1]] <- data.frame(

TRAb = tr, totT4 = tail(d$T4, 1), totT3 = tail(d$T3, 1),

freeT3 = tail(d$fT3, 1), Fd = tail(d$Fd, 1), Sig = tail(d$Sig, 1),

store_d = tail(d$storeD, 1), TSH = tail(d$TSH, 1),

Tc = tail(d$Tc, 1), HR = tail(d$HR, 1), NEFA = tail(d$NEFA, 1),

BWPS = tail(d$BWPS, 1), Lambda = tail(d$Lambda, 1),

Tmax72 = max(d$Tc), storm = max(d$Tc) >= 40)

}

bind_rows(rows)

1.3, 1.5, 2.0),

TRAb = 4.19, hours = 168) {

base <- ts_thyrotoxic(TRAb)

bind_rows(lapply(precs, function(pr) {

d <- ts_run(list(ev = ts_none(), p = list()), base = base, TRAb = TRAb,

hours = hours, prec = pr, pge = 0.6 * pr)

data.frame(Prec = pr, Tmax = max(d$Tc),

Tc_24h = d$Tc[which.min(abs(d$time - 24))],

HRmax = max(d$HR), BWPS_24h = d$BWPS[which.min(abs(d$time - 24))],

BWPSmax = max(d$BWPS), Bilimax = max(d$Bili),

mort_pct = 100 * max(d$mort),

storm = max(d$Tc) >= 40)

}))

TRAb = 4.19, hours = 168) {

base <- ts_thyrotoxic(TRAb)

arms <- list(

propranolol = list(ev = c(ev_ptu(), ev_iod(1), ev_pro(), ev_hc()),

p = list(COOL = 1.2, FLUID = 0.0045)),

esmolol = list(ev = c(ev_ptu(), ev_iod(1), ev_hc()),

p = list(COOL = 1.2, FLUID = 0.0045, ESMRATE = 1.60)))

bind_rows(lapply(CRs, function(cr) {

b <- base; b$CR <- cr

bind_rows(lapply(names(arms), function(a) {

d <- ts_run(arms[[a]], base = b, TRAb = TRAb, hours = hours)

data.frame(CR0 = cr, arm = a, CO_nadir = min(d$perf),

shock_max = max(d$shock),

BWPS_24h = d$BWPS[which.min(abs(d$time - 24))],

Bilimax = max(d$Bili), mort_pct = 100 * max(d$mort))

}))

}))

base <- ts_thyrotoxic(TRAb)

d <- ts_run(list(ev = ev_iod(1), p = list()), base = base, TRAb = TRAb,

hours = 24 * days, delta = 1, prec = 0, pge = 0)

d %>% filter(time %% 24 == 0) %>%

transmute(day = time / 24, Ipl, NIS, Ithy, relinh, A_TPO = 1 - inhTPO,

Synth, S, totT4 = T4, totT3 = T3)

nefas = c(.4, .6, .8, 1.0, 1.5, 2.0, 2.5, 3.0)) {

base <- ts_thyrotoxic(TRAb)

bind_rows(lapply(nefas, function(nf) {

b <- base; b$NEFA <- nf

d <- ts_mod %>% update(init = b) %>% param(TRAb = TRAb) %>%

mrgsim(end = 0, delta = 1) %>% as.data.frame()

data.frame(NEFA = nf, measured_totT3 = d$T3[1], Fd = d$Fd[1],

freeT3 = d$fT3[1], Sig = d$Sig[1],

equivalent_totT3 = d$T3[1] * d$Fd[1])

}))

keep = c("1 Untreated", "2 PTU alone",

"3 Methimazole alone", "5 Propranolol alone",

"8 Cooling + fluids alone", "9 FULL BUNDLE")) {

d <- res[res$scenario %in% keep, ]

op <- par(mfrow = c(2, 3), mar = c(4, 4, 2.5, 1))

vars <- list(c("Tc", "core temperature (degC)"),

c("HR", "heart rate (bpm)"),

c("fT3", "FREE T3 (pmol/L)"),

c("T3", "TOTAL T3 (nmol/L)"),

c("BWPS", "Burch-Wartofsky score"),

c("mort", "cumulative mortality"))

cols <- seq_along(keep)

for (v in vars) {

plot(NA, xlim = range(d$time), ylim = range(d[[v[1]]]), xlab = "time (h)",

ylab = v[2], main = v[2])

for (i in seq_along(keep)) {

s <- d[d$scenario == keep[i], ]

lines(s$time, s[[v[1]]], col = cols[i], lwd = 2)