#  Characterization of differences in immune responses during bolus and
continuous infusion endotoxin challenges using mathematical modeling

Kristen A. Windoloski  Department of Mathematics, North Carolina State
University, Raleigh, North Carolina, USA Susanne Janum  Frederiksberg and
Bispebjerg Hospitals, Frederiksberg, Denmark  Centre for Physical Activity
Research, Rigshospitalet, Denmark  Ronan M.G. Berg  Department of Clinical
Physiology and Nuclear Medicine, Copenhagen University Hospital, Denmark
Centre for Physical Activity Research, Rigshospitalet, Denmark  Department of
Biomedical Sciences, University of Copenhagen, Denmark  Neurovascular Research
Laboratory, University of South Wales, Pontypridd, UK
Mette S. Olufsen  Department of Mathematics, North Carolina State University,
Raleigh, North Carolina, USA

###### Abstract

Endotoxin administration is commonly used to study the inflammatory response,
and though traditionally given as a bolus injection, it can be administered as
a continuous infusion over multiple hours. Several studies hypothesize that
the latter better represents the prolonged and pronounced inflammation
observed in conditions like sepsis. Yet, very few experimental studies have
administered endotoxin using both strategies, leaving significant gaps in
determining the underlying mechanisms responsible for their differing immune
responses. We use mathematical modeling to analyze cytokine data from two
studies administering a 2 ng/kg dose of endotoxin, one as a bolus and the
other as a continuous infusion over four hours. Using our model, we simulate
the dynamics of mean and subject-specific cytokine responses as well as the
response to long-term endotoxin administration. Cytokine measurements reveal
that the bolus injection leads to significantly higher peaks for IL-8, while
IL-10 reaches significantly higher peaks during continuous administration.
Moreover, the peak timing of all measured cytokines occurs later in the
continuous infusion. We identify three model parameters that significantly
differ between the two administration methods. Monocyte activation of IL-10 is
greater during the continuous infusion, while recovery rates of IL-8 is faster
for the bolus injection. This suggests that a continuous infusion elicits a
stronger, longer-lasting systemic reaction through increased stimulation of
monocyte anti-inflammatory mediator production and decreased recovery of pro-
inflammatory catalysts. Furthermore, our continuous infusion model exhibits
prolonged inflammation with recurrent peaks resolving within two days during
long-term (20-32 hours) endotoxin administration.

Key Points

  * •

This study uses mathematical modeling to compare the dynamic response to
continuous and bolus endotoxin administration examining the hypothesis that
continuous infusion better represents the inflammation seen in clinical
scenarios such as sepsis.

  * •

Our modeling study augments limited experimental studies by providing a better
understanding of the pathways and mechanisms impacting the change in immune
responses with long-term endotoxin infusion.

  * •

We introduce a mathematical model incorporating essential cellular and
cytokine pathways and calibrate it to mean and subject-specific data from two
endotoxin studies.

  * •

Statistical analysis of optimized model parameters suggests that the monocyte
activation rate for IL-10 is greater for the continuous infusion and
degradation rates of TNF-$\alpha$ and IL-8 are greater for the bolus
injection.

  * •

Continuous infusion model simulations show that prolonged LPS administration
(20-32 hours) significantly extends the cytokine response and results in
recurrent cytokine spikes. After 32 hours, the stimuli becomes too weak to
elicit a system response.

## Introduction

Endotoxin (lipopolysaccharide, LPS) derived from gram-negative bacteria’s
outward membrane (Heine et al., , 2001) is an immunostimulant administered to
healthy subjects as an experimental procedure to study the inflammatory
response (Suffredini and Noveck, , 2014). This type of experiment, referred to
as an endotoxin challenge, has allowed insight into mechanisms and treatments
of inflammation events such as rheumatoid arthritis (Miller et al., , 1979;
Merrill et al., , 2011; Lorenz et al., , 2013), systemic lupus erythematosus
(Zuckerman et al., , 1996), cancer (Easson et al., , 1998; Ho et al., , 2013;
Yassine, , 2016), Alzheimer’s disease (Sly et al., , 2001; Cunningham et al.,
, 2005; Akimoto et al., , 2007), and sepsis (Fitzal et al., , 2003;
Fredriksson et al., , 2009; Shinozaki et al., , 2010; Leijte et al., , 2019).

In an endotoxin challenge, LPS can be administered as a bolus (instantaneous)
injection (Copeland et al., , 2005; Clodi et al., , 2008; Janum et al., ,
2016), a continuous infusion over several hours (Berg et al., , 2012), or a
combination of the two (Kiers et al., , 2017) in both humans and animals
(Bahador and Cross, , 2007). The response is an increase in pro-
(TNF-$\alpha$, IL-1$\beta$, IL-6, IL-8) and anti- (IL-10, IL-1ra) inflammatory
cytokines, immune cells, body temperature, heart rate, blood pressure, and
hormone levels (Givalois et al., , 1994; Bahador and Cross, , 2007; Clodi et
al., , 2008; Janum et al., , 2016). The peak of each measured quantity and the
time it takes to reach the peak vary depending on a host of controllable
(administration method and total dose administered) and uncontrollable
(individual variation due to genetics, sex, and health status) factors.

Taudorf et al.,  (2007) performed an endotoxin challenge in healthy men
administering 0.3 ng/kg of LPS as a bolus and a continuous infusion over 4
hours. They found that the administration method significantly affects
TNF-$\alpha$, IL-6, and neutrophil production rates. These measured quantities
peaked earlier and had larger magnitudes during the bolus administration than
the continuous infusion. Kiers et al.,  (2017) compared immune responses to 1
and 2 ng/kg bolus doses of LPS in addition to a 1 ng/kg bolus followed by a 3
ng/kg continuous infusion over 3 hours. This study showed a significant
difference in mean cytokine concentrations, flu-like symptoms (headache,
nausea, shivering, pain), temperature, and heart rate increases between the
bolus-only and the bolus plus continuous infusion dose. Cytokines responses
(TNF-$\alpha$, IL-6, IL-8, IL-10) reached significantly higher peak levels,
and subjects exhibited prolonged elevated flu-like symptoms during the bolus
plus continuous infusion method. These results demonstrate that continuous
infusion initiates a more durable and occasionally more pronounced impact on
the immune response during the incitement of inflammation.

Although experimental studies are suitable for investigating the effects of
endotoxin administration method, they do not provide insight into why
differing dynamics are observed. This is where the power of mathematical
modeling of physiological systems can be applied. Simulations with
mathematical models can highlight the underlying mechanisms of disease, aid in
disease diagnosis, test and validate treatments and predict patient trajectory
and mortality. Numerous mathematical models of inflammation have been
developed over the last two decades. Kumar et al.,  (2004); Day et al.,
(2006); Reynolds et al.,  (2006) developed small but novel mathematical
models, highlighting their ability to reproduce inflammation scenarios of
clinical relevance and potential to predict treatment strategies. Several
models built upon this foundation by adding specific immune cells and
cytokines activated during the inflammatory response (Chow et al., , 2005; Roy
et al., , 2007; Foteinou et al., , 2009; Su et al., , 2009; Parker et al., ,
2016; Brady et al., , 2018; Torres et al., , 2019). Others created detailed
model incorporating feedback from other physiological entities such as the
cardiovascular system, nervous system, the hypothalamic-pituitary-adrenal
(HPA) axis, pain perception, and thermal responses (Scheff et al., , 2010;
Foteinou et al., , 2011; Malek et al., , 2015; Bangsgaard et al., , 2017;
Dobreva et al., , 2021; Windoloski et al., , 2023). Several models were
calibrated to experimental data or validated in specific patients. Some used
data from a bolus endotoxin challenge in animals (mice or rats) (Chow et al.,
, 2005; Day et al., , 2006; Roy et al., , 2007; Parker et al., , 2016; Torres
et al., , 2019) while others used bolus data from human subjects (Foteinou et
al., , 2009; Scheff et al., , 2010; Foteinou et al., , 2011; Malek et al., ,
2015; Bangsgaard et al., , 2017; Brady et al., , 2018; Dobreva et al., , 2021;
Windoloski et al., , 2023). These studies demonstrate the need for
computational inflammation models that (i) utilize experimental data from a
continuous infusion of endotoxin and (ii) investigate the mechanisms behind
response differences observed during variations in endotoxin administration
method.

Recent experimental studies (Kiers et al., , 2017; van Lier et al., , 2019)
propose that a continuous endotoxin infusion is more appropriate to study the
prolonged system response during systemic inflammation and sepsis. To provide
more insight into understanding what immune signaling components are impacted
during the switch from a bolus to continuous infusion, we study the
inflammatory response to continuous infusion of endotoxin through the lens of
a mathematical model. Doing so provides (i) newfound insight into the response
differences between a bolus and continuous administration of endotoxin, (ii) a
better mathematical representation to study the dynamics of sepsis, and (iii)
a better model to investigate treatments of inflammatory conditions since the
continuous infusion prolongs the exposure window for treatment testing.

We present a novel inflammatory response mathematical model predicting innate
cytokine responses (TNF-$\alpha$, IL-6, IL-8, IL-10) to a 2 ng/kg bolus and
continuous infusion over four hours of endotoxin from Janum et al.,  (2016)
and Berg et al.,  (2012). The model structure is rigorously explored through
sensitivity and identifiability analysis, and parameter estimation calibrates
the model to mean and subject-specific cytokine data. We compare each study’s
cytokine data to characterize larger endotoxin doses than compared in previous
literature and develop statistical uncertainty bounds for the optimal mean
model. Mechanisms responsible for varying immune dynamics observed in bolus
and continuous infusion experimental studies are hypothesized via statistical
analysis of optimized model parameters. We propose that the transition from a
bolus to continuous infusion impacts physiologically-relevant components
related to IL-10 activation by monocytes and TNF-$\alpha$ and IL-8 degradation
rates. Moreover, we use our continuous infusion model to investigate the
system response to perturbations in infusion duration and total endotoxin dose
administered. This illustrates its ability as a clinically-realistic in silico
model that can simulate prolonged and pronounced responses.

## Methods

### Ethical approval

The current work utilized experimental data from two published studies by Berg
et al.,  (2012) and Janum et al.,  (2016). The study by Berg et al.,  (2012)
was approved by the Scientific Ethical Committee of Copenhagen and
Frederiksberg Municipalities in Denmark. The data from Berg et al.,  (2012)
was made available by Berg (coauthor of this study). The study by Janum et
al.,  (2016) received approval for the experimental protocol by the Regional
Committee on Health Research Ethics and the Regional Monitoring Board, and the
study followed the protocols listed in the Declaration of Helsinki. Individual
data from Janum et al.,  (2016) was made available by Janum (coauthor of this
study) and Mehlsen (coauthor of Janum et al.,  (2016)). All participants from
both studies gave their written and oral consent.

### Experimental data

The experiments by Berg et al.,  (2012) and Janum et al.,  (2016) administered
the same total dose of endotoxin (2 ng/kg) to healthy study participants, one
as a continuous infusion and one as a bolus injection. Mean and subject-
specific cytokine data were measured and used to calibrate our mathematical
model.

The study from Berg et al.,  (2012) investigated the effects of an increase in
mean arterial pressure on cerebral autoregulation; it included nine healthy
male participants aged 21-25. All study participants were subject to physical
examination. Data were only included from subjects with normal blood work and
cardiovascular markers. Participants did not take medication, had a typical
medical history, were non-sedentary, and infection-free at least four weeks
before the study. The study by Janum et al.,  (2016) was designed to
investigate the connection between pain and the innate immune system reaction
in 20 male athletes aged 18-35. All study participants had a healthy weight,
were non-smokers, and had no signs of illness two weeks prior to the study
day. Pre-screening activities involved a review of each subject’s medical
history, a physical examination, and laboratory work.

In Berg et al.,  (2012), participants were subject to a 4-hour continuous
infusion of 2 ng/kg (0.5 ng/kg/hr) of endotoxin (Batch G2 B274, US
Pharmacopeial Convention, Rockville, MD, USA) administered via an antecubital
catheter. In contrast, the study participants in Janum et al.,  (2016)
received a 2 ng/kg bolus endotoxin dose (Lot EC-6, National Institutes of
Health, Bethesda, MD, USA) via a peripheral intravenous catheter following two
hours of baseline stabilization. In Berg et al.,  (2012), blood samples were
taken hourly for the first 4 hours following the start of endotoxin
administration and 2 hours after completed endotoxin administration.
Measurements from Janum et al.,  (2016) were taken hourly, starting two hours
before endotoxin administration and continuing for six hours following
administration. To capture peak response, this study analyzed an additional
blood sample taken 1.5 hours after LPS administration. Janum et al.,  (2016)
used ELISA (Meso Scale Discovery, Rockville, Maryland, USA) and Berg et al.,
(2012) used SECTOR Imager 2400 (Meso Scale Diagnostics, Gaithersburg, MD, USA)
to determine concentrations of TNF-$\alpha$, IL-6, IL-8, and IL-10.

Figure 1: Box and whisker plots of the continuous infusion ($m=9$) (Berg et
al., , 2012) and bolus ($n=20$) (Janum et al., , 2016) data. The horizontal
axis represents the number of hours after the start of endotoxin
administration. Black boxplots above the symbol ‘B’ represent bolus data, and
blue boxplots above the symbol ‘C’ represent continuous infusion data. The red
cross symbol denotes abnormal responses (outliers) from each study. This
figure is generated using MATLAB code adapted from Danz,  (2023).

Subjects 1, 2, 6, 8, and 9 from the continuous infusion study were missing one
cytokine measurement, and subject 4 was missing four measurements. Of these,
subjects 4, 8, and 9 were missing baseline concentration measurements of IL-6.
Subject 4 also did not have a baseline concentration of IL-10. Figure 1 shows
data from both studies, identifying outliers from both data sets. Because the
immune response exhibits significant variation in individual responses to
stimuli, we considered these outlying data points abnormal but not
unrealistic.

Table 1: Experimental data characteristics. Abnormal concentrations (outliers) were not used to compute the mean data. Individual subject concentrations were calculated from a sample size $n=20$ for the bolus data from Janum et al.,  (2016), and $m=9$ for the continuous infusion data from Berg et al.,  (2012). All concentrations were rounded to the nearest whole number. Measurement | Study | TNF-$\alpha$ | IL-6 | IL-8 | IL-10
---|---|---|---|---|---
Peak of Mean Concentration (pg/mL) | Bolus | 326 | 702 | 714 | 40
Continuous | 532 | 707 | 456 | 136
Peak of Subject Concentration (pg/mL) | Bolus | 60-1297 | 351-1856 | 522-1124 | 20-105
Continuous | 212-1293 | 303-1335 | 170-683 | 57-376
Peak Timing of Mean Concentration (pg/mL) | Bolus | 1.5 | 2 | 2 | 3
Continuous | 3 | 4 | 4 | 4
Peak Timing of Subject Concentration (pg/mL) | Bolus | 1.5-2 | 2-3 | 1.5-3 | 2-3
Continuous | 3-4 | 4 | 3-4 | 3-6

We calculated the mean cytokine response after removing abnormal responses
(outlying data points outside the 1.5$\times$IQR range (the 25th and 75th
percentile of data) for each endotoxin data set. Because of the small number
of study participants, we only removed abnormal (outlying) measurements
instead of that individual’s entire cytokine profile. We used the mean of the
bolus and continuous infusion data to calibrate our mathematical model. In the
remainder of this study, we refer to this as the mean bolus or continuous
endotoxin administration.

Mean and subject-specific cytokine characteristics are reported in Table 1.
Primary pro- and anti-inflammatory cytokines TNF-$\alpha$ and IL-10 had higher
mean concentrations during continuous infusion, while secondary cytokines IL-6
did not depend on the administration method, but IL-8 had a higher peak value
for bolus injection. With continuous infusion, peak concentrations were later
for all measured cytokines. Individual subject concentrations from both
studies displayed a considerable variation in cytokine responses;
TNF-$\alpha$, IL-6, and IL-8 had a higher variance with bolus injection.

### Data calibration

To compare the cytokine responses from the two studies, we adjusted the bolus
data so that both studies had the same baseline concentration. This was done
by determining the difference ($\bar{d}$) between the mean bolus ($\bar{b}$)
and continuous infusion ($\bar{c}$) baseline

$\displaystyle\bar{d}^{j}_{0}=\bar{b}^{j}_{0}-\bar{c}^{j}_{0}$

for cytokine $j=\\{TNF,IL6,IL8,IL10\\}$ and then shifting the concentrations
by

$\displaystyle\hat{b}^{j}_{i}(k)=b^{j}_{i}(k)-\bar{d}^{j}_{0},$

where $b^{j}_{i}(k)$ is the original cytokine concentration $j$ at time $i$
for the $k$th bolus participant ($1\leq k\leq 20$). The adjusted cytokine
concentration is denoted by $\hat{b}^{j}_{i}(k)$. Figure 2 displays the mean
and subject-specific continuous infusion and bolus data.

Figure 2: Mean and subject-specific continuous infusion and bolus data. Red
circles denote the continuous infusion data ($m=9$) connected by solid lines,
and black squares denote the bolus data ($n=20$) connected by dotted lines.
Thin lines represent subject-specific responses and thick vertical lines
denote mean (SD).

### Mathematical model

Our mathematical model (Figure 3) adapted from our previous studies (Brady et
al., , 2018; Dobreva et al., , 2021; Windoloski et al., , 2023) predicting
dynamics of the innate immune response to endotoxin included a system of seven
ordinary differential equations (ODEs) with 45 parameters. The equations
characterized the time-varying concentrations of endotoxin, resting and
activated monocytes, and pro- and anti-inflammatory cytokines. Below we
briefly describe the model and refer to Brady,  (2017) for a detailed
derivation of the equations.

Figure 3: Cytokine interactions. Endotoxin administration recruits (activates)
monocytes from a large pool of resting monocytes ($M_{R}$). The active
monocytes ($M_{A}$) upregulates production of pro- (TNF-$\alpha$, IL-6, and
IL-8) and anti- (IL-10) inflammatory cytokines. TNF-$\alpha$ generates
positive feedback on monocyte, IL-6, and IL-8 production. IL-6 exhibits anti-
inflammatory properties through self-regulation, downregulation of
TNF-$\alpha$, and upregulation of IL-10. IL-10 downregulates all cytokine
production and the activation of monocytes. Solid black lines represent
stimulation, and dotted black lines represent inhibition.

Endotoxin: The equation determining the endotoxin concentration ($E$, ng/kg)
was adapted from Brady et al.,  (2018); Dobreva et al.,  (2021) to account for
a continuous infusion. This formulation is similar to that in Windoloski et
al.,  (2023) and motivated by Day et al.,  (2006). The endotoxin rate of
change was given by

$\frac{dE}{dt}=\begin{cases}D_{h}-k_{E}E&t\leq D_{ad}\\\
-k_{E}E&t>D_{ad}\end{cases},$ (1)

where $D_{h}$ (ng/kg/hr) was the endotoxin dose administered per hour,
$D_{ad}$ (hr) the dosing administration duration, and $k_{E}$ (hr-1) the
endotoxin decay rate. The 2 ng/kg continuous infusion was administered over 4
hours so $D_{h}=0.5$, $D_{ad}=4$, and $E(0)=0$, whereas $D_{h}=0$, $D_{ad}=0$,
and $E(0)=2$ for the bolus injection.

Monocytes: During the endotoxin challenge, resting monocytes circulating in
the blood are activated, upregulating cytokine production. This process
regulates inflammation via positive and negative feedback (Rossol et al., ,
2011). The resting ($M_{R}$) and activated ($M_{A}$) monocytes (number of
cells, noc) were found from

$\frac{dM_{R}}{dt}=k_{MR}M_{R}\left(1-\frac{M_{R}}{M_{\infty}}\right)-H_{M}^{U}(E)\left(k_{M}+k_{MTNF}H_{M}^{U}(TNF)\right)H_{M}^{D}(IL10)M_{R}\\\
$ (2)
$\frac{dM_{A}}{dt}=H_{M}^{U}(E)\left(k_{M}+k_{MTNF}H_{M}^{U}(TNF)\right)H_{M}^{D}(IL10)M_{R}-k_{MA}M_{A},$
(3)

where $k_{MR}$ (hr-1) denoted the regeneration rate and $M_{\infty}$ (noc) the
carrying capacity for the resting monocytes. The activated monocytes were
upregulated by endotoxin (Rossol et al., , 2011) at rate $k_{M}$ (hr-1) and
inflammatory cytokine TNF-$\alpha$ at rate $k_{MTNF}$ (hr-1). They were also
downregulated by anti-inflammatory cytokine IL-10 (Kucharzik et al., , 1998).
This process was activated relative to the resting monocytes. The increase in
activated monocytes caused an identical decrease in the resting monocytes.
Finally, the activated monocytes decayed at rate $k_{MA}$ (hr-1).

In equations (2)-(3), the upregulation or downregulation of state $Y$ by state
$X$ was as

$H^{U}_{Y}(X)=\frac{X^{h}}{\eta_{YX}^{h}+X^{h}},\ \ \
H^{D}_{Y}(X)=\frac{\eta_{YX}^{h}}{\eta_{YX}^{h}+X^{h}},$

where $h$ represents the steepness of the curve and $\eta_{YX}$ the half-
maximum value.

Inflammatory mediators: Activated monocytes upregulate cytokines, signaling
proteins that promote or suppress inflammation (Murphy, , 2012). Cytokines
that stimulate inflammation, called pro-inflammatory, include TNF-$\alpha$,
IL-6, and IL-8 (Johnston and Webster, , 2009). TNF-$\alpha$ is an early pro-
inflammatory mediator responsible for the induction of fever (Murphy, , 2012)
and recruitment of other pro-inflammatory cytokines (Johnston and Webster, ,
2009). IL-6 is a secondary pro-inflammatory mediator primarily involved in the
induction of the liver acute phase response (Johnston and Webster, , 2009),
although it can exhibit anti-inflammatory properties as well (Tilg et al., ,
1994). IL-8 is a late pro-inflammatory mediator mainly responsible for
recruiting neutrophils to the target site (Bickel, , 1993). Monocytes also
release anti-inflammatory cytokines, particularly IL-10, to counteract pro-
inflammatory responses and provide a balanced immune response (Johnston and
Webster, , 2009). These four cytokines are essential components of the innate
immune response. Their interactions can be predicted by

$\frac{dTNF}{dt}=k_{TNFM}H_{TNF}^{D}(IL6)H_{TNF}^{D}(IL10)M_{A}-k_{TNF}(TNF-
w_{TNF})\ $ (4)
$\frac{dIL6}{dt}=\left(k_{6M}+k_{6TNF}H_{IL6}^{U}(TNF)\right)H_{IL6}^{D}(IL6)H_{IL6}^{D}(IL10)M_{A}-k_{6}(IL6-w_{6})$
(5)
$\frac{dIL8}{dt}=\left(k_{8M}+k_{8TNF}H_{IL8}^{U}(TNF)\right)H_{IL8}^{D}(IL10)M_{A}-k_{8}(IL8-w_{8})$
(6)
$\frac{dIL10}{dt}=\left(k_{10M}+k_{106}H_{IL10}^{U}(IL6)\right)M_{A}-k_{10}(IL10-w_{10}).$
(7)

In equation (4), TNF-$\alpha$ was activated by monocytes (Johnston and
Webster, , 2009) at rate $k_{TNFM}$ (pg (mL hr noc)-1) and downregulated by
IL-6 and IL-10 (Tilg et al., , 1994). In equation (5), IL-6 was activated by
monocytes and TNF-$\alpha$ (Johnston and Webster, , 2009) at rates $k_{6M}$
(pg (mL hr noc)-1) and $k_{6TNF}$ (pg (mL hr noc)-1), and downregulated by
itself (Verboogen et al., , 2019) and IL-10 (Johnston and Webster, , 2009).
Similarly in equation (6), IL-8 was activated by monocytes and TNF-$\alpha$ at
rates $k_{8M}$ (pg (mL hr noc)-1) and $k_{8TNF}$ (pg (mL hr noc)-1), and
downregulated by IL-10 (Johnston and Webster, , 2009). In equation (7), IL-10
was activated by monocytes and IL-6 at rates $k_{10M}$ (pg (mL hr noc)-1) and
$k_{106}$ (pg (mL hr noc)-1), (Murphy, , 2012; Jin et al., , 2013). Cytokines
decayed to their baseline concentrations $w_{i}$ (pg (mL)-1) at rate $k_{i}$
(hr-1) for $i=\\{TNF,6,8,10\\}$.

Model summary: The mathematical model described above was an ODE system of the
form

$\frac{dX}{dt}=f(t,X;\theta),$ (8)

where $X\in\mathbb{R}^{7}$ denoted the time-varying states
$X=\\{E,M_{R},M_{A},TNF,IL6,IL8,IL10\\}$ determining endotoxin ($E$),
monocytes (resting $M_{R}$ and activated $M_{A}$), TNF-$\alpha$, IL-6, IL-8,
and IL-10 concentrations. The model parameters $\theta\in\mathbb{R}^{45}$ are
listed in Table 2 with subsections indicating what state the parameters belong
to.

To fit the mathematical model to the experimental data, we minimized the least
squares cost function, $J$, given by

$J=r^{T}r,\ \ \ r=[r_{TNF}\hskip 14.22636ptr_{IL6}\hskip
14.22636ptr_{IL8}\hskip 14.22636ptr_{IL10}],$ (9)

where $r_{k}$ is the residual vector for each state $k=\\{TNF,IL6,IL8,IL10\\}$
and

$r_{k}=\frac{1}{\sqrt{N}}\left(\frac{\left[y^{k}_{1}\ldots
y^{k}_{N}\right]-y^{k}_{data}}{\max\left(y^{k}_{data}\right)}\right).$ (10)

In equation (10), $N$ refers to the number of data points for each state $k$,
$y^{k}_{i}=g(t_{i},X_{k}(t_{i});\theta)$ denotes the model output for the
state $k$ at time $t_{i}$, $1\leq i\leq N$, and $y^{k}_{data}$ is the
associated data. The least squares cost $J$ was minimized using the nonlinear
optimization solver, fmincon, from MATLAB (MathWorks Inc., Natick, MA, USA).
Upper and lower parameter bounds were set by multiplying and dividing the
parameters’ nominal value by a factor of four.

Parameter | Description | Unit |  | Nominal
---
Value (C)
| Nominal
---
Value (B)
Endotoxin ($E$)
$D_{h}$ |  $E$ administered pr hour |  ng (kg hr)-1 | 0.5 | 0
$D_{ad}$ |  $E$ administered duration | hr | 4 | 0
$k_{E}$ |  $E$ decay rate |  hr-1 | 1.01 | 1.01
Monocytes ($M$)
$k_{MR}$ | Regeneration rate |  hr-1 | 0.006 | 0.006
$k_{MA}$ | Activated decay rate |  hr-1 | 2.51 | 2.51
$k_{MTNF}$ |  Activation rate by $TNF$ |  hr-1 | 9.000 | 8.650
$k_{M}$ |  Activation rate by $E$ |  hr-1 | 0.041 | 0.041
$\eta_{ME}$ |  Upregulation half max of $E$ |  ng (kg)-1 | 2 | 3.3
$\eta_{M10}$ |  Downregulation half max of $IL10$ |  pg (mL)-1 | 4.394 | 3.884
$\eta_{MTNF}$ |  Upregulation half max of $TNF$ |  pg (mL)-1 | 222.222 | 140.845
$h_{ME}$ |  Upregulation exp of $E$ on $M$ | non dim | 1 | 1
$h_{M10}$ |  Upregulation exp of $IL10$ on $M$ | non dim | 0.3 | 0.3
$h_{MTNF}$ |  Upregulation exp of TNF on $M$ | non dim | 3.16 | 3.16
$M_{\infty}$ | Carrying capacity | noc | 30000 | 30000
TNF-$\alpha$ ($TNF$)
$k_{TNF}$ | Decay rate |  hr-1 | 0.6 | 1
$k_{TNFM}$ |  Activation rate by $M$ |  pg (mL hr noc)-1 | 1.333 | 0.845
$\eta_{TNF10}$ |  Downregulation half max of $IL10$ |  pg (mL)-1 | 17.576 | 15.536
$\eta_{TNF6}$ |  Downregulation half max of $IL6$ |  pg (mL)-1 | 560 | 560
$h_{TNF10}$ |  Downregulation exp of $IL10$ on $TNF$ | non dim | 3 | 3
$h_{TNF6}$ |  Upregulation exp of $IL6$ on $TNF$ | non dim | 2 | 2
$w_{TNF}$ | Baseline concentration |  pg (mL)-1 | 13.873 | 8.793
IL-6 ($IL6$)
$k_{6}$ | Decay rate |  hr-1 | 0.66 | 0.66
$k_{6M}$ |  Activation rate by $M$ |  pg (mL hr noc)-1 | 0.81 | 0.81
$k_{6TNF}$ |  Activation rate by $TNF$ |  pg (mL hr noc)-1 | 0.81 | 0.81
$\eta_{610}$ |  Downregulation half max of $IL10$ |  pg (mL)-1 | 35.152 | 31.071
$\eta_{66}$ |  Downregulation half max of $IL6$ |  pg (mL)-1 | 560 | 560
$\eta_{6TNF}$ |  Upregulation half max of $TNF$ |  pg (mL)-1 | 411.11 | 260.563
$h_{610}$ |  Downregulation exp of $IL10$ on $IL6$ | non dim | 1 | 4
$h_{66}$ |  Downregulation exp of $IL6$ on $IL6$ | non dim | 1 | 1
$h_{6TNF}$ |  Upregulation exp of $TNF$ on $IL6$ | non dim | 2 | 2
$w_{6}$ | Baseline concentration |  pg (mL)-1 | 0.610 | 0.610
IL-8 ($IL8$)
$k_{8}$ | Decay rate |  hr-1 | 0.66 | 0.66
$k_{8M}$ |  Activation rate by $M$ |  pg (mL hr noc)-1 | 0.509 | 0.789
$k_{8TNF}$ |  Activation rate by $TNF$ |  pg (mL hr noc)-1 | 0.509 | 0.789
$\eta_{810}$ |  Downregulation half max of $IL10$ |  pg (mL)-1 | 17.576 | 15.536
$\eta_{8TNF}$ |  Upregulation half max of $TNF$ |  pg (mL)-1 | 411.11 | 260.563
$h_{810}$ |  Downregulation exp of $IL10$ on $IL8$ | non dim | 1.5 | 1.5
$h_{8TNF}$ |  Upregulation exp of $TNF$ on $IL8$ | non dim | 3 | 3
$w_{8}$ | Baseline concentration |  pg (mL)-1 | 2.695 | 4.175
IL-10 ($IL10$)
$k_{10}$ | Decay rate |  hr-1 | 0.4 | 0.8
$k_{10M}$ |  Activation rate by $M$ |  pg (mL hr noc)-1 | 0.019 | 0.017
$k_{106}$ |  Activation rate by $IL6$ |  pg (mL hr noc)-1 | 0.019 | 0.017
$\eta_{106}$ |  Upregulation half max of $IL6$ |  pg (mL)-1 | 560 | 560
$h_{106}$ |  Upregulation exp of $IL6$ on $IL10$ | non dim | 3.68 | 3.68
$w_{10}$ | Baseline concentration |  pg (mL)-1 | 4.239 | 3.747
Table 2: Parameter descriptions, units, and nominal values for the continuous
infusion (C) and bolus injection (B) endotoxin models.

### Nominal parameters

Nominal parameter values were taken from Brady,  (2017) except for cytokine
baseline concentrations ($w_{TNF},w_{6},w_{8}$ and $w_{10}$), which were set
to the mean continuous infusion and bolus data. We manually adjusted nominal
parameters affecting peak timing to account for the observation (Figure 1)
that the timing of cytokine activation depends on the administration method.
To improve the nominal model fit to the peak magnitudes of cytokine profiles,
the peak concentration of state $i=\\{TNF,IL6,IL8,IL10\\}$, denoted $X_{i}$,
was scaled using the technique from Windoloski et al.,  (2023) where

$X_{i}=\alpha\tilde{X_{i}},$ (11)

for the scaling factor $\alpha$ and desired peak concentration
$\tilde{X_{i}}$. We substituted equation (11) into the ODE for state $X_{i}$
giving

$\frac{dX_{i}}{dt}=f\left(t,\theta,X_{i},X_{j}\right)\hskip
14.22636pt\implies\hskip
14.22636pt\frac{d(\alpha\tilde{X_{i}})}{dt}=f\left(t,\theta,\alpha\tilde{X_{i}},X_{j}\right)$

for states $j\neq i$. Therefore,

$\frac{d\tilde{X_{i}}}{dt}=\frac{1}{\alpha}f\left(t,\theta,\alpha\tilde{X_{i}},X_{j}\right).$
(12)

The scaling factor $1/\alpha$ was distributed to each term on the right side
of the ODE, scaling the associated parameters in each term. State $X_{i}$ was
also scaled when it was upregulating another state variable, $Y$, as

$H^{U}_{Y}(X_{i})=H^{U}_{Y}(\alpha\tilde{X_{i}})=\frac{(\alpha\tilde{X_{i}})^{h}}{\eta_{YX_{i}}^{h}+(\alpha\tilde{X_{i}})^{h}}\hskip
14.22636pt\implies\hskip
14.22636ptH^{U}_{Y}(\tilde{X_{i}})=\frac{\tilde{X_{i}}^{h}}{\left(\frac{\eta_{YX_{i}}}{\alpha}\right)^{h}+\tilde{X_{i}}^{h}}.$

Thus, half-saturation values were scaled by $1/\alpha$. A similar approach was
applied for downregulation functions. Baseline cytokine parameters $(w_{i})$
were also scaled. Table 2 lists the nominal parameters for the continuous
infusion and bolus mean model.

Most nominal parameter values used for model calibration to the individual
data were set to the mean optimal values except for the initial cytokine
concentrations ($w_{TNF},w_{6},w_{8}$ and $w_{10}$), which were set to the
individual’s cytokine value at baseline. For subjects missing measurements at
time zero, we scaled their concentration after one hour based on values from
subjects where data were available; IL-6 and IL-10 baseline values were set at
$50\%$ and $75\%$ of their concentrations at hour one. The scaling analysis
was also applied to the individual subjects because of the significant
individual variation in cytokine responses between study participants. To
reduce the number of scaled parameters in the subject-specific optimizations,
we only scaled cytokines with scaling factor $\alpha<0.9$ or $\alpha>1.1$.

### Sensitivity analysis and subset selection

The highly nonlinear mathematical model had seven states and 45 parameters.
Because of its structure (Brady, , 2017; Brady et al., , 2018) and the
quantity of data, we selected a parameter subset from the rate constants to
estimate. We first conducted a local relative sensitivity analysis as
described in Olufsen and Ottesen,  (2013) on the mean continuous infusion
response using the residual vector in equation (9). The sensitivity matrix
$\chi$ was given by

$\chi=\frac{\partial r}{\partial\log(\theta)}=\frac{\partial
y}{\partial\theta}\cdot\frac{\theta}{\max(y_{data})},$ (13)

where $y=g(t,X(t),\theta)$ was the model output at time $t$, $\theta$ the
nominal parameter set, and $y_{data}$ the mean continuous infusion data. We
approximated the $(i,j)$ entry in the submatrix $\chi_{k}$ using forward
differences, where

$\chi=\left[\begin{array}[]{c}\chi_{TNF}\hskip 8.53581pt\chi_{IL6}\hskip
8.53581pt\chi_{IL8}\hskip 8.53581pt\chi_{IL10}\end{array}\right]^{T}.$ (14)

For submatrix $X_{k}$, elements $\chi_{ij}$ were given by

$\chi_{ij}=\frac{g\left(t_{i},X_{k}(t_{i}),\theta+he_{j}\right)-g\left(t_{i},X_{k}(t_{i}),\theta\right)}{h}\frac{\theta}{\max\left(y^{k}_{data}\right)},$
(15)

where $\phi=10^{-8}$ was the solver tolerance, $h=\sqrt{\phi}$ the step size
(Pope et al., , 2009), and $e_{j}$ the basis vector in the $j$th direction. We
ranked relative sensitivities by computing the two-norm of each column of
$\chi$, obtaining a single sensitivity per parameter. We repeated the
sensitivity analysis by simulating 100 runs sampling parameters from a uniform
distribution varying $\pm 30\%$ around the parameter’s nominal value to study
effects due to perturbations in parameter values.

Sensitive rate constants were used to select an identifiable parameter subset
that can be estimated. We utilized two practical identifiability techniques,
the structured correlation method (SCM) and the SVD-QR method (Miao et al., ,
2011; Olufsen and Ottesen, , 2013). The SCM used the Fisher-information matrix
$F=\chi^{T}\chi$. We checked the condition number to ensure that $F$ had an
inverse and calculated $G=F^{-1}$. The matrix $G$ was used to determine the
pairwise parameter covariance $C_{ij}$ by

$C_{ij}=\frac{G_{ij}}{\sqrt{G_{ii}\ G_{jj}}},$ (16)

where $(i,j)$ refers to $\theta_{i}$ and $\theta_{j}$. Parameter pairs for
which $|C_{ij}|>0.9$ were considered correlated. The parameter set with the
largest correlation was selected. The parameter within that set with the
smallest relative sensitivity was removed from the parameter set, and the
process was repeated until there were no correlated parameters.

The SVD-QR method used singular value decomposition (SVD) to determine
identifiable parameters. This method decomposed the sensitivity matrix
$\chi=U\Sigma V^{T}$ where $U$ and $V$ contained the left and right singular
vectors of $\chi$, and $\Sigma$ contained the singular values $\sigma$ of
$\chi$. The largest $k$ singular values of $\chi$ were determined by
$\sigma(k)\geq 10\sqrt{\phi}$ where $\phi$ was the ODE solver tolerance. The
first $k$ columns of the right singular vectors, $V_{k}$, were extracted from
$V$ and used to find a permutation matrix $P$ such that $V_{k}^{T}P=QR$, where
$Q$ was an orthogonal matrix and $R$ was an upper triangular matrix. The
permutation matrix $P$ was then used to reorder the parameter vector $\theta$
as

$\hat{\theta}=P^{T}\theta.$ (17)

The first $k$ parameters of $\hat{\theta}$ were considered identifiable.

The subset selection methods determine identifiable parameters near the
nominal values. To ensure that optimal values were also identifiable, we
investigated the parameter convergence for each subset. We conducted 20
optimizations for each subset with nominal parameters drawn from a uniform
distribution of $\pm 10\%$ of each estimated parameter’s nominal value. All
other parameters were fixed. For each of the 20 runs, we calculated the
coefficient of variation ($CoV$) for each estimated parameter $\theta_{i}$
where

$CoV\left(\theta_{i}\right)=\frac{\bar{\theta}_{i}}{\sigma_{i}}.$ (18)

We denoted $\bar{\theta}_{i}$ as the mean and $\sigma_{i}$ as the standard
deviation of the estimated parameter $\theta_{i}$. Parameters in each subset
with $CoV\geq 0.1$ were identified. The least sensitive parameter was removed
from the set, and this process was repeated until all estimated parameters in
each subset had a $CoV<0.1$. The resulting parameter subsets were considered
sensitive and identifiable and were used for parameter estimation.

### Statistical methods

For each parameter subset, the goodness of fit was computed using the
coefficient of determination (R2) (Dodge, , 2008), the corrected Akaike
information criterion (AICc) (Burnham and Anderson, , 2002), and the Bayesian
information criterion (BIC) (Schwarz, , 1978). Details of these measurements
are provided in Section 2 of the Supporting Information.

Additionally, we constructed parameter and model confidence and prediction
intervals using the frequentist approach detailed in Seber and Wild,  (2003);
Banks et al.,  (2009); Smith,  (2013). Parameter confidence intervals for
optimized parameter $\tilde{\theta}_{i}$ were computed as

$\tilde{\theta}_{i}\pm t_{N-q}^{\alpha/2}\sqrt{\Sigma_{ii}},$ (19)

where $N$ was the total number of data points, $q$ was the number of
parameters that were estimated, $t_{N-q}^{\alpha/2}$ was the t-value from the
student’s t-distribution for confidence level $1-\alpha$ with $N-q$ degrees of
freedom, and the variance estimator matrix $\Sigma$ was given by

$\Sigma=(\chi^{T}(\tilde{\theta})V^{-1}\chi(\tilde{\theta}))^{-1}.$ (20)

We defined $\chi$ similarly to equations (14) and (15) where

$\chi(\tilde{\theta})=\left[\begin{array}[]{c}\chi_{TNF}(\tilde{\theta})\hskip
8.53581pt\chi_{IL6}(\tilde{\theta})\hskip
8.53581pt\chi_{IL8}(\tilde{\theta})\hskip
8.53581pt\chi_{IL10}(\tilde{\theta})\end{array}\right]^{T}$ (21)

and the $(i,j)$ element of submatrix $\chi_{k}(\tilde{\theta})$ with
$k\in\\{TNF,IL6,IL8,IL10\\}$ was approximated using forward differences given
by

$\chi_{ij}(\tilde{\theta})=\frac{g\left(t_{i},X_{k}(t_{i}),\tilde{\theta}+he_{j}\right)-g\left(t_{i},X_{k}(t_{i}),\tilde{\theta}\right)}{h}.$
(22)

Here, $\tilde{y}^{k}_{i}=g(t_{i},X_{k}(t_{i}),\tilde{\theta})$ was the optimal
model output with optimal parameter vector $\tilde{\theta}$ for cytokine state
$k$ at time $t_{i}$ for $1\leq i\leq N_{k}$ where $N_{k}$ was the number data
points for cytokine $k$, and $h$ and $e_{j}$ were defined as in equation (15).
The diagonal variance matrix $V$ was given by

$V=\text{diag}\left(\sigma_{TNF},\sigma_{IL6},\sigma_{IL8},\sigma_{IL10}\right),$
(23)

where $\sigma_{k}$ was a diagonal matrix of size $N_{k}\times N_{k}$ with
entries

$\frac{1}{N_{k}-q}\left(r_{k}^{T}r_{k}\right),\ \ \
r_{k}=\left[\tilde{y}^{k}_{1}\ldots\tilde{y}^{k}_{N_{k}}\right]-y^{k}_{data}$
(24)

with $y^{k}_{data}$ defined in equation (10). The asymptotic prediction
interval for cytokine $k$ at time $t_{i}$ was given by

$PI^{k}_{i}=\tilde{y}^{k}_{i}\pm
t_{N_{k}-q_{k}}^{\alpha/2}s_{k}\sqrt{1+\textbf{G}_{ik}^{T}\left(\chi_{k}^{T}(\tilde{\theta})\chi_{k}(\tilde{\theta})\right)^{-1}\textbf{G}_{ik}}$
(25)

and the confidence interval by

$CI^{k}_{i}=\tilde{y}^{k}_{i}\pm
t_{N_{k}-q_{k}}^{\alpha/2}s_{k}\sqrt{\textbf{G}_{ik}^{T}\left(\chi_{k}^{T}(\tilde{\theta})\chi_{k}(\tilde{\theta})\right)^{-1}\textbf{G}_{ik}}.$
(26)

We defined $\tilde{y}^{k}_{i}$ and $N_{k}$ as in equation (22), $q_{k}$ was
the number of estimated parameters that impacted cytokine state $k$, and
$t_{N_{k}-q_{k}}^{\alpha/2}$ was the t-value for confidence level $1-\alpha$
with $N_{k}-q_{k}$ degrees of freedom. $\chi_{k}(\tilde{\theta})$ was given by
(21) and (22), but columns in $\chi_{TNF}(\tilde{\theta})$,
$\chi_{IL6}(\tilde{\theta})$, and $\chi_{IL10}(\tilde{\theta})$ corresponding
to IL-8 parameters were eliminated since they did not impact those state
variables (see Figure 3). Entries of these columns were approximately zero and
made $F_{k}=\chi_{k}^{T}(\tilde{\theta})\chi_{k}(\tilde{\theta})$ singular
unless removed. The matrix $\textbf{G}_{ik}^{T}$ was defined as

$\textbf{G}_{ik}^{T}=\left(\frac{\partial\tilde{y}^{k}_{i}}{\partial\theta_{1}}\ldots\frac{\partial\tilde{y}^{k}_{i}}{\partial\theta_{q_{k}}}\right),$
(27)

which was $i$th row of the submatrix $\chi_{k}(\tilde{\theta})$, and the
variance estimator $s^{2}_{k}$ was given by

$s^{2}_{k}=\frac{1}{N_{k}-q_{k}}r_{k}^{T}r_{k},$ (28)

with $r_{k}$ in equation (24). Due to the small number of data points per
cytokine, we generated pseudodata to compute uncertainty intervals. Data
points at eight and twelve hours were set by quartering the cytokine
concentration at six hours and returning the cytokine to the baseline value.
Then, a piecewise cubic spline interpolation was performed from $t=0-12$
hours.

Statistical data analysis included a two-sample unequal variances t-test
($\alpha=0.05$) on the continuous infusion and bolus data before data
calibration to compare their maximal concentrations and peak timing
statistically. Abnormal cytokine responses (outliers in Figure 1) were omitted
from the data sampled to conduct the hypothesis test. A two-sample unequal
variances t-test ($\alpha=0.05$) was also performed on the set of optimized
parameters from subject-specific optimizations to determine statistically
significant differences in parameter values between the two administration
methods. Parameter values that were outliers within their data set were not
included in the data sampled to conduct the hypothesis test.

## Results

###  Data

Statistical comparison (Table 3) of the continuous infusion and bolus
injection data show a significantly smaller concentration of IL-8
($p=0.00147$) and larger concentration of IL-10 ($p=0.00200$) with continuous
infusion. The peak concentration for TNF-$\alpha$ ($p=0.0809$) and IL-6
($p=0.702$) did not statistically differ significantly between the two
studies, but the time to peak cytokine concentration was significantly longer
for all cytokines during the continuous infusion study: TNF- $\alpha$, IL-6,
and IL-8 ($p<0.0001$) and IL-10 ($p=0.00695$).

Cytokine |  | Maximal
---
Concentration
| Time to Maximal
---
Concentration
TNF-$\alpha$ |  $p=0.0809$ ($m=9,n=19$) |  $p<0.0001$ ($m=9,n=19$)
IL-6 |  $p=0.702$ ($m=9,n=18$) |  $p<0.0001$ ($m=9,n=18$)
IL-8 |  $p=0.00147$ ($m=9,n=19$) |  $p<0.0001$ ($m=9,n=19$)
IL-10 |  $p=0.00200$ ($m=8,n=17$) |  $p=0.00695$ ($m=8,n=18$)
Table 3: Statistical significance $(\alpha=0.05$) of data attributes between
the continuous infusion ($m$ subjects) and bolus ($n$ subjects) studies in
Berg et al.,  (2012) and Janum et al.,  (2016). Subjects with abnormal
responses (outliers) for either quantity were omitted from the sample.

### Sensitivity analysis and subset selection

Single and repeated sensitivity analysis (Figure 4) highlights the system’s
dependence on endotoxin, activated monocytes, TNF-$\alpha$, and IL-10
dynamics. The system’s most sensitive parameters are the growth or decay of
these states, where $k_{TNFM}$ (growth rate of TNF-$\alpha$ by monocytes) and
$k_{MA}$ (activated monocyte decay rate) have the most significant impact on
the model. This can be explained by endotoxin and activated monocytes
promoting the activation of cytokines, where TNF-$\alpha$ and IL-10 are the
main cytokines that upregulate and downregulate other states. The least
sensitive rate constant is $k_{MR}$, the regeneration rate for resting
monocytes. Given that our study administers a finite dose of endotoxin that
does not deplete the resting monocytes before the system can recover, it is
reasonable that this parameter has a minute effect on the system. The single
and repeated sensitivity analysis results exhibit similar behavior with minor
variations in the order of sensitivity. This observation and careful scaling
of nominal parameter values provides a good foundation for choosing
identifiable subsets among the sensitive parameters.

Figure 4: Ranked sensitivities. Local sensitivities, scaled by the maximum
sensitivity, are denoted by yellow crosses. Boxplots of scaled relative
sensitivities are generated from $n=100$ local sensitivity analysis
simulations. Values are scaled by the maximum average sensitivity. Black dots
denote outliers.

The parameter $k_{MR}$ is insensitive, removed from the subset, and fixed at
its nominal value. Identifiability analysis using the SCM and SVD-QR method is
performed on the remaining 15 sensitive rate constants. Because different
identifiability analysis methods are not guaranteed to produce the same
results (Brady, , 2017), we generate three parameter subsets, two using only
the SCM and one using the SVD-QR followed by the SCM. Results show that all
rate constants for monocytes, IL-6, IL-8, and IL-10 cannot be uniquely
estimated. Therefore, the sensitive rate constants are split into two subsets,
one that includes monocyte-activated growth rates and the other that contains
cytokine-activated growth rates. The SCM performed on each of these subsets
results in the two subsets

$\displaystyle S_{1}$
$\displaystyle=\\{k_{E},k_{MA},k_{MTNF},k_{TNF},k_{TNFM},k_{6},k_{6TNF},k_{8},k_{8TNF}\\}$
$\displaystyle S_{2}$
$\displaystyle=\\{k_{MA},k_{M},k_{TNF},k_{TNFM},k_{6},k_{6M},k_{8},k_{10},k_{10M}\\}.$

The third subset is found by performing SVD-QR followed by the SCM on all 15
rate constants, which results in the subset

$\displaystyle S_{3}$
$\displaystyle=\\{k_{MA},k_{TNF},k_{TNFM},k_{6},k_{6M},k_{8},k_{8M},k_{10M}\\}.$

The identifiability and convergence of the above parameter subsets are checked
numerically using the coefficient of variation method, enabling us to reduce
the subsets further, obtaining three sensitive and identifiable parameter
subsets

$\displaystyle s_{1}$
$\displaystyle=\\{k_{MA},k_{MTNF},k_{TNF},k_{TNFM},k_{6},k_{8},k_{8TNF}\\}$
$\displaystyle s_{2}$
$\displaystyle=\\{k_{MA},k_{M},k_{TNF},k_{TNFM},k_{6},k_{6M},k_{8},k_{10},k_{10M}\\}$
$\displaystyle s_{3}$
$\displaystyle=\\{k_{MA},k_{TNF},k_{TNFM},k_{8},k_{8M},k_{10M}\\}.$

Note that $S_{2}=s_{2}$, indicating that the subset $S_{2}$ was identifiable.

### Parameter estimation and uncertainty quantification

Model fit for the mean continuous infusion data (R2, AICc, BIC, and least
squares cost $J$) for subsets $s_{1},s_{2},$ and $s_{3}$ are reported in Table
4. Subset $s_{3}$ has the lowest AICc and BIC values, but the R2 value and
least squares cost did not differ significantly between the three subsets.
Given the significance of the AICc and BIC values, we conduct the remaining
simulations using $S_{Final}=s_{3}$ including

$\displaystyle S_{Final}=\\{k_{MA},k_{TNF},k_{TNFM},k_{8},k_{8M},k_{10M}\\}.$

| Subset
---
Estimated
| Number of
---
Parameters
Average R2 | AICc | BIC | $J$
$s_{1}$ | 7 | 0.915 | 15.4 | 16.6 | 0.0602
$s_{2}$ | 9 | 0.923 | 24.7 | 22.5 | 0.0466
$s_{3}$ | 6 | 0.913 | 11.1 | 13.3 | 0.0547
Table 4: Goodness of fit measurements for the optimized subsets $s_{1},s_{2},$
and $s_{3}$. The coefficient of determination is denoted as $R^{2}$, AICc
represents the corrected Akaike information criterion, BIC the Bayesian
information criterion, and $J$ the least squares cost. Figure 5: Optimized
model fit to mean continuous infusion and bolus data estimating $S_{Final}$.
The mean continuous infusion fit is marked by red solid lines, and the mean
(SD) of the data by red circles and error bars. The mean bolus fit has black
dotted lines, and black triangles and error bars denote the mean (SD) of the
data.

The mean continuous infusion model exhibits later activation of monocytes and
cytokines compared to the bolus injection model (Figure 5). As a result, the
main pro- and anti-inflammatory cytokines TNF-$\alpha$ and IL-10 have larger
peak concentrations. The immune resolution time during the continuous infusion
model is approximately ten to twelve hours, whereas the mean bolus model is
only six to eight hours. Comparison of model fits by the coefficient of
determination (R2) for each cytokine reveal that TNF-$\alpha$ and IL-6 are
fitted better by the bolus model, while the continuous infusion model better
predicts IL-8 and IL-10. Differences are minor, though, specifically for IL-8
and IL-10.

We generated $N=61$ data points for each cytokine to determine the mean data
and model uncertainty using confidence level $(1-\alpha)$ with $\alpha=0.05$.
Confidence bounds on the optimal parameters from the continuous infusion and
bolus mean model responses are given in Table 5. The upper and lower bounds
remain within the physiological values except for $k_{10M}$, which has a
negative lower bound. Prediction and confidence intervals on the optimal mean
model are shown in Figure 6. Both prediction and confidence intervals for the
bolus (Figure 6b) are tighter than those for the continuous infusion model
(Figure 6a), indicating the variability of mean measurements and model output
is larger in the continuous infusion data. This is plausible, given the sample
sizes of the two studies. The lower bound for the prediction intervals of both
dose types extends into negative cytokine values, which is mathematically but
not physiologically appropriate.

Figure 6: $95\%$ prediction and confidence intervals for the mean (A) continuous infusion and (B) bolus responses. Mean and SD data points are marked by circles and error bars. Prediction intervals are the yellow dashed-dotted lines, and confidence intervals are the blue dotted lines. Parameter |  | Continuous Infusion
---
(Optimal Value $\pm$ Bound)
| Bolus
---
(Optimal Value $\pm$ Bound)
$k_{MA}$ | $3.49\pm 0.0994$ | $2.67\pm 0.0787$
$k_{TNF}$ | $0.423\pm 0.132$ | $1.40\pm 0.118$
$k_{TNFM}$ | $1.39\pm 0.0696$ | $0.998\pm 0.0398$
$k_{8}$ | $0.386\pm 0.119$ | $0.686\pm 0.100$
$k_{8M}$ | $0.613\pm 0.193$ | $0.746\pm 0.163$
$k_{10M}$ | $0.0365\pm 0.127$ | $0.0150\pm 0.124$
Table 5: Optimal $95\%$ parameter confidence bounds for the mean continuous
infusion and bolus model. Figure 7: Optimal model responses for (A) subject 1
from the continuous infusion study and (B) subject 16 from the bolus study.

We fit the model to the subject-specific cytokine profiles from the continuous
infusion ($m=9$) and bolus injection ($n=20$) studies by estimating the
parameters in $S_{Final}$. Results for continuous infusion subject 1 and bolus
injection subject 16 are shown in Figure 7, and dynamics for the remaining
subjects are presented in Figures S1-S29 in the Supporting Information.
Results show that our model captures varying cytokine responses to the same
total dose of endotoxin for both administration methods. While individual peak
cytokine concentrations and peak timing differ from that in the mean response,
the model (shown in Figures 7a and 7b) is sufficiently robust to capture
variation in data. This is evidenced by high R2 values for all subjects.

The mean and standard deviation for subject-specific optimal and scaled
parameters are listed in Table 6, and a boxplot of the optimized subject-
specific parameter values are shown in Figures 8A and B. We observe similar
median parameter values for the continuous infusion and bolus subject-specific
optimizations for parameters $k_{MA},k_{TNFM},$ and $k_{8M}$. For parameters
$k_{MA},k_{TNFM}$, and $k_{10M}$, there is a larger variance in the continuous
infusion than the bolus injection. Optimized parameter values denoted as
outliers in Figures 8A and B correspond to subjects 3, 5, and 9 from the
continuous infusion study and subjects 3, 9, 13, 14, and 20 from the bolus
study. These subjects all had abnormal endotoxin responses (at least one
outlying data point in Figure 1). Boxplots of all scaled subject-specific
parameters are shown in Figure 8C and all subject-specific parameter values
are listed in Table 6.

Statistical comparison of the continuous infusion and bolus optimized
parameters show that $k_{TNF}$ and $k_{8}$ ($p<0.0001$) are significantly
larger during the bolus injection, indicating the TNF-$\alpha$ and IL-8 decay
faster during the bolus dose. Additionally, $k_{10M}$ ($p=0.0142$) was
significantly larger during the continuous infusion, implying that monocyte
activation of IL-10 is more pronounced during a continuous infusion of
endotoxin. As a result, the continuous infusion had a significantly larger
activation response of IL-10 by monocytes and substantially smaller
TNF-$\alpha$ and IL-8 degradation rates. Parameters $k_{MA}$ ($p=0.465$),
$k_{TNFM}$ ($p=0.106$), and $k_{8M}$ ($p=0.0615$) were not significantly
different between the two administration methods, as reported in Table 6.
Abnormal responses denoted as outliers in Figure 8 were not included in the
sample from each study.

Parameter |  | Continuous Infusion
---
Mean (SD)
| Bolus
---
Mean (SD)
P-Value
$k_{MA}$ | 3.88 (1.06) | 3.59 (0.587) | $p=0.465(m=9,n=20)$
$\eta_{M10}$ | 5.69 (3.38) | 5.72 (2.79) |
$\eta_{MTNF}$ | 249 (152) | 152 (114) |
$k_{TNF}$ | 0.642 (0.274) | 1.78 (0.406) | $p<0.0001(m=9,n=19)$
$k_{TNFM}$ | 2.25 (1.56) | 1.44 (0.948) | $p=0.106(m=9,n=19)$
$\eta_{TNF10}$ | 22.8 (13.5) | 22.9 (11.2) |
$\eta_{TNF6}$ | 674 (356) | 686 (346) |
$w_{TNF}$ | 7.22 (5.09) | 6.86 (5.22) |
$k_{6M}$ | 0.975 (0.516) | 0.992 (0.500) |
$k_{6TNF}$ | 0.975 (0.516) | 0.992 (0.500) |
$\eta_{610}$ | 45.5 (27.1) | 45.8 (22.3) |
$\eta_{66}$ | 674 (356) | 686 (346) |
$\eta_{6TNF}$ | 461 (280) | 280 (211) |
$w_{6}$ | 0.658 (0.484) | 0.831 (0.690) |
$k_{8}$ | 0.439 (0.127) | 0.719 (0.174) | $p<0.0001(m=8,n=19)$
$k_{8M}$ | 0.733 (0.331) | 0.941 (0.386) | $p=0.0615(m=8,n=19)$
$k_{8TNF}$ | 0.542 (0.202) | 0.884 (0.154) |
$\eta_{810}$ | 22.8 (13.5) | 22.9 (11.2) |
$\eta_{8TNF}$ | 461 (280) | 280 (211) |
$w_{8}$ | 5.03 (5.73) | 3.33 (1.07) |
$k_{10M}$ | 0.0518 (0.0339) | 0.0238 (0.012) | $p=0.0142(m=7,n=19)$
$k_{106}$ | 0.0250 (0.0148) | 0.0251 (0.012) |
$\eta_{106}$ | 674 (356) | 686 (346) |
$w_{10}$ | 5.82 (5.04) | 6.280 (3.24) |
Table 6: Subject-specific parameter values as mean (SD) and estimated
parameter p-values with significance level $\alpha=0.05$. Estimated parameters
are marked in bold and remaining parameters were scaled from their nominal
values. Mean (SD) values are computed using subject parameter values,
including abnormal responses. P-values were calculated by removing the
abnormal responses prior to hypothesis testing, with $m$ subjects from the
continuous infusion and $n$ subjects from the bolus study being included.
Figure 8: (A) Boxplots of subject-specific optimized parameters from the bolus
‘B’ (black, $n=20$ subjects) and continuous ‘C’ (blue, $m=9$ subjects)
administration models. Associated p-values are listed next to each parameter,
(B) zoomed in boxplot of optimized parameter $k_{10M}$ from (A), and (C)
boxplots of subject-specific scaled parameters. On all plots outliers are
denoted by the red cross. Parameters considered statistically significant
($\alpha=0.05$) include $k_{TNF}$ ($p<0.0001$), $k_{8}$ ($p<0.0001$), and
$k_{10M}$ ($p=0.0142$). Parameters not statistically significant include
$k_{MA}$ ($p=0.465$), $k_{TNFM}$ ($p=0.106$), and $k_{8M}$ ($p=0.0615$). This
figure is generated using MATLAB code adapted from Danz,  (2023).

### Infusion perturbations

We use the optimal mean continuous infusion model to study the response to a
longer duration of inflammation and enhanced immune stimulation. Figure 9A
shows the model response when 2 ng/kg of endotoxin is given continuously over
4, 8, 12, and 24 hours. The infusion duration impacts peak cytokine
concentrations and the response’s resolution time. Peak concentrations
declined and occurred later as the infusion duration increased. Cytokine
concentrations returned to baseline approximately 12, 16, 20, and 36 hours
following the infusion start for the 4, 8, 12, and 24-hour continuous
infusions. The system exhibits oscillatory behavior when the infusion is
extended to 24 hours. The increase of anti-inflammatory cytokine IL-10 around
10 hours combats the initial pro-inflammatory response of TNF-$\alpha$ to
decline around 12 hours. However, because the endotoxin is still being
administered, it rebounds a second time once IL-10 levels begin to decline.
Following the termination of endotoxin administration, the monocytes are no
longer activated, and as a result, the inflammatory markers return to
baseline. This recurrent inflammatory behavior transpires when endotoxin is
administered for 20 to 32 hours, after which the stimulation from the
endotoxin is not strong enough to induce a pronounced response (Figure S35 in
the Supporting Information). This simulation also shows that the system takes
approximately 21-23 days to recover (Figure S36 in the Supporting Information)
relative to the resting monocyte population returning to the baseline value.

Figure 9: Continuous infusion mean model response when (A) 2 ng/kg of
endotoxin is administered as a 4 (red solid lines), 8 (yellow dashed-dotted
lines), 12 (light blue dotted lines), and 24-hour (dark blue dashed lines)
continuous infusion, and (B) 2 (red solid lines), 4 (yellow dashed-dotted
lines), 8 (light blue dotted lines), and 16 ng/kg (dark blue dashed lines) of
endotoxin is administered as a 4-hour continuous infusion.

Figure 9B displays the model response for a 4-hour continuous endotoxin
infusion of 2, 4, 8, and 16 ng/kg. The total endotoxin dose impacts peak
cytokine concentrations and the immune resolution time. Larger doses of
endotoxin result in earlier, greater peak cytokine concentrations, which occur
approximately 1.5-2 hours before peak cytokine concentrations for smaller
endotoxin doses. Additional simulations increasing both the duration of the
continuous infusion and the total dose of endotoxin are shown in Figures
S30-S34 the Supporting Information.

## Discussion

This study uses mathematical modeling to compare bolus and continuous
administration of LPS. The model is calibrated to data from Berg et al.,
(2012) and Janum et al.,  (2016). Data analysis reveals that IL-10 has a
significantly higher peak for the continuous dose, while the peak IL-8
concentration is higher with the bolus injection. For the continuous dose, the
peaks appear significantly later for all cytokines, a trend that continues
when the dose is given over a longer time. Model parameter analysis provide
insight into what processes may change with administration methods. Our
results suggest that the continuous infusion of endotoxin increases the
monocyte production rate of anti-inflammatory cytokine IL-10 and decreases the
clearance rates of significant pro-inflammatory markers TNF-$\alpha$ and IL-8.
Our continuous infusion model is crucial as it can replicate characteristics
of clinical inflammation associated with prolonged elevation of immune markers
when endotoxin infusion is extended and pronounced cytokine responses when the
endotoxin dosage is increased. Interestingly, administration over 20 and 32
hours produces double cytokine peaks indicating inflammation recurrence.
Finally, we found that it takes over 20 days before the resting monocytes have
reached the same level as before the stimulus.

Our findings agree with observations in Kiers et al.,  (2017) comparing a 2
ng/kg bolus response to a 1 ng/kg bolus followed by a 3 ng/kg continuous
infusion. A bolus injection followed by a continuous infusion showed higher
IL-10 production and prolonged symptoms due to an extended elevation of
cytokines. They also reported a significantly higher production of
TNF-$\alpha$, IL-6, and IL-8 with a bolus injection plus continuous infusion.
Our study exhibited an increased TNF-$\alpha$ production, though results were
not statistically significant, likely due to the small number of subjects and
high variance between subjects. IL-8 had significantly lower peaks during the
continuous infusion (Figure 5). Our findings also agree with the partial
conclusion of Taudorf et al.,  (2007), who reported that (i) the release of
cytokines TNF-$\alpha$ and IL-6 occurred significantly later with the
continuous infusion compared to a bolus injection and (ii) TNF-$\alpha$ and
IL-6 concentrations were significantly larger for the bolus dose. For our
study, IL-6 was higher for the bolus dose but again, results were not
significant. Taudorf et al.,  (2007) also reported larger neutrophil
concentrations that peaked earlier with the bolus dose. Our study did not
account for neutrophil dynamics, a component that could be added in future
studies. Differences in our findings could be a result of low endotoxin dosage
in Taudorf et al.,  (2007) and unequal total endotoxin dosing in Kiers et al.,
(2017). Overall, our findings implicate that continuous stimulation of the
system over hours could promote a more significant anti-inflammatory response
to counteract prolonged levels of pro-inflammatory cytokines, leading to lower
maximal concentrations of secondary cytokines such as IL-8.

A significant contribution of this study is the statistical analysis of
optimal parameter distributions between the two administration methods, which
has not been examined in earlier works. Results show that the activation rate
of IL-10 by monocytes was significantly larger, and the TNF-$\alpha$ and IL-8
decay rates were substantially lower in the continuous infusion versus the
bolus injection. These results suggest that continual endotoxin infusion
amplifies the monocyte production of IL-10 and dulls the resolution of pro-
inflammatory cytokines TNF-$\alpha$ and IL-8. Kiers et al.,  (2017) indicates
that a continuous infusion of endotoxin is a more probable model of prolonged
inflammation in conditions like sepsis, where a hyperinflammatory state (often
referred to as a cytokine storm) is accompanied by an immunosuppressive phase
with elevated anti-inflammation levels (Nedeva et al., , 2019; Torres et al.,
, 2022). Thus, stimulation by a continuous infusion of endotoxin may exhibit
mild but clear signs of a prolonged pro-inflammatory response and a
hyperactive anti-inflammatory response, similar to dynamics observed in sepsis
and supporting the hypothesis of Kiers et al.,  (2017).

Model analysis demonstrated the reliance of dynamics on the endotoxin,
monocyte, TNF-$\alpha$, and IL-10 states. These constituents encompass primary
elements of the inflammatory response - immune cells that respond to stimuli
and the main pro- and anti-inflammatory cytokines that modulate the response
strength (Johnston and Webster, , 2009). Thus, it is plausible that these
components strongly dictate immune dynamics. It is also reasonable that the
least influential component is the monocyte regeneration rate since the
challenges analyzed here are from short-lived low-dose endotoxin exposure and
the monocyte pool is not depleted prior to endotoxin clearance. However, if
simulating a pathogenic insult, we suspect that the influence of this
parameter on system dynamics would significantly increase to clear an
infection of much greater magnitude than is safely observed in an endotoxin
challenge.

In Figure 8, several parameter values were marked as outliers for continuous
infusion and bolus subject-specific model fits. These parameters correspond to
subjects from both studies that exhibited abnormal cytokine responses for at
least one of the measured cytokines. Given this, we hypothesize that these
subjects may experience a more severe response or even enhanced complications
to a clinical inflammation event. While this requires further investigation,
it could be explored in silico by mathematical modeling.

We also explore variations of endotoxin infusion duration and total dose in
Figure 9. A comparable simulation was conducted in Windoloski et al.,  (2023)
on a bolus endotoxin model where the total dose was increased, representing
the administration of more potent immune stimuli that cannot safely be given
to humans and the stimuli strength of clinical infection. Both our study and
Windoloski et al.,  (2023) observed enhanced cytokine production. Our
simulations extending the continuous infusion duration correspond to the
clinical scenario of continual systemic aggravation by inflammatory stimuli.
In this case, the model produces up to approximately 36 hours of elevated
immune markers depending on the length of the endotoxin infusion. It also
displays attributes similar to endotoxin tolerance, a clinical phenomenon
related to a reduced response to endotoxin after initial exposure (West and
Heagy, , 2002), through the appearance of multiple decreasing cytokine peaks
when continual endotoxin administration is given across 20 to 32 hours.
Oscillations in cytokine concentrations also arise for 24 and 36-hour
infusions when the total endotoxin dose is increased from 2 ng/kg to 4, 8, and
16 ng/kg (Figures S33-S34 in Supporting Information), showing that the system
can produce fluctuating behavior for prolonged periods of inflammation if the
stimuli are large enough. These cytokine oscillations that occur could also be
clinically-relevant with reference to recurrent infections where the system
returns close to baseline before peaking again. In a clinical setting,
however, the inflammatory stimulus is a live pathogen whose concentration
would also fluctuate, compared to an endotoxin challenge where the
administration of endotoxin is constant until cessation of infusion. This
model of prolonged inflammation can be used to study inflammatory dynamics
over longer periods and test or validate treatments for inflammatory
conditions given the longer endotoxin exposure window.

###  Limitations

A major limitation is that our model is calibrated to two data sets from Berg
et al.,  (2012) and Janum et al.,  (2016). Although both administered a total
dose of 2 ng/kg of endotoxin and had similar experimental protocols, the
endotoxin was sourced from different vendors. There is widespread individual
variation in the human immune response, evidenced by the individual subject
data used in this study (Figures 1 and 2). However, this is not uncommon. It
is well-known that immune responses vary between individuals due in part to
uncontrollable factors such as genetics, age, sex, seasonal and circadian
influences, and environmental effects (Brodin and Davis, , 2017). Therefore,
utilizing additional endotoxin challenge data or, more ideally, comparing the
immune responses during both endotoxin administration strategies on the same
subjects using the same endotoxin batch would yield the best results. Only a
few studies administer large endotoxin doses (such as 2 ng/kg as used here) as
both a bolus and a continuous infusion. While Taudorf et al.,  (2007)
administers 0.3 ng/kg of endotoxin as a bolus and a continuous infusion, the
cytokine concentrations are notably lower than those from a larger dose
(Krabbe et al., , 2001; Janum et al., , 2016) and near or below reported
concentrations in septic patients (Casey et al., , 1993; Wu et al., , 2009;
Berg et al., , 2012). Torres et al.,  (2022) suggests that most patients are
likely in the immunosuppressive stage of sepsis upon hospital admittance, so
we suspect initial inflammation levels could be higher than reported in sepsis
studies. Therefore data from a 2 ng/kg endotoxin challenge likely yield more
realistic cytokine concentrations as observed in sepsis and should be used in
our study.

Another limitation is that we do not have enough data to validate our
endotoxin perturbation results on the continuous infusion model. Experimental
data for a continuous infusion of larger endotoxin doses is not seen in
literature except in Kiers et al.,  (2017) (who administers a total of 4 ng/kg
of endotoxin) since larger doses of endotoxin are considered unsafe (Bahador
and Cross, , 2007). Safety may also play a role in the lack of experimental
studies administering endotoxin for an extensive time beyond 4 hours. Another
limitation is that, although our mathematical model is highly nonlinear and
complex, there are direct elements of the immune response (cells such as
macrophages and neutrophils, cytokines such as IL-1$\beta$ and TGF-$\beta$,
and signaling pathways such as the NF-$\kappa$B pathway) and other sources of
immune regulation (cardiovascular, nerve, hormonal, metabolic) that are not
included here. Although Janum et al.,  (2016) reports that IL-1$\beta$ was
measured in the bolus study, its concentrations were not detectable. Our
previous work (Dobreva et al., , 2021) explored interactions of immune,
cardiovascular, thermal, and pain responses during a bolus endotoxin
challenge, and Windoloski et al.,  (2023) expanded on that model to include
hormonal regulation. Including these additional factors in the model dynamics
could provide clearer insight into processes that activate at different speeds
or strengths when the endotoxin challenge administration method is varied
between a bolus and continuous infusion. A deeper understanding of continuous
infusion dynamics could provide a better translational mathematical model of
systemic inflammation such as sepsis, encompassing multi-organ dynamics.

### Future work

Further investigation of this work involves expanding our study of continuous
infusion dynamics to include immune interactions with other systems, such as
the cardiovascular and neuroendocrine systems, thermal, pain, and metabolic
regulation, building upon the work in Windoloski et al.,  (2023). These
components are well-known to impact immune response and regulation (Miller et
al., , 2010; Hjemdahl et al., , 2011; Kenney and Ganta, , 2014; Janum et al.,
, 2016; Varela et al., , 2018; Dobreva et al., , 2021), and studying how a
continuous infusion affects these elements can enhance understanding of
clinically prolonged inflammation events. Furthermore, while an endotoxin
challenge attempts to mimic the dynamics of a clinical-level immune insult,
its duration is finite. It cannot simulate the extensive effects of an actual
infection or trauma. Therefore, mathematical modeling can extrapolate dynamics
from controlled environments to clinical relevance by looking at the impact of
age, smoking, diabetes, or cancer on immune responses. Additional future
directions of this study focus on transitioning our endotoxin immune response
model to a model of sepsis, a life-threatening condition involving hyperactive
immune responses and subsequent organ failure that is still not fully
understood (Nedeva et al., , 2019). Much research has been focused on
identifying a universal biomarker and treatment of sepsis, but one has yet to
be accepted within the scientific community (Cecconi et al., , 2018). However,
recent progress has proposed several candidates, including administering
vitamin C (Kashiouris et al., , 2020; Wald et al., , 2022) to sepsis patients.
Adapting our current model to a model of sepsis could help improve our
understanding of the mechanisms of sepsis and provide insight into the
efficacy of potential sepsis treatments.

## Conclusion

To enhance understanding of potential mechanisms impacting immune responses to
endotoxin, we devised a physiologically-based mathematical model simulating
mean and subject-specific dynamics in human volunteers exposed to continuous
and bolus endotoxin administration. Comparison of subject-specific optimized
parameter values revealed significant differences in the monocyte activation
rate of IL-10 and recovery rates of pro-inflammatory cytokines TNF-$\alpha$
and IL-8. This suggests that increased IL-10 activation by monocytes and
slower recovery rates of pro-inflammatory cytokines could play a role in the
more pronounced anti-inflammatory response and smaller secondary cytokine
response seen in the continuous infusion data. Additionally, these factors
likely influence the system’s elongated and more gradual response to the
endotoxin, as seen by the statistically significant later peak concentration
times of all cytokines during the continuous infusion. Individuals with
abnormal cytokine responses also reported statistically outlying optimal
parameter values, suggesting their responses to a clinical infection could
result in enhanced (outlying) reactions and complications. Simulations of the
continuous infusion for a longer duration or increased dose amount display the
model’s capability to predict immune responses to prolonged inflammation or
more potent inflammatory stimuli. Future directions of this work focus on
including a whole-body response model to study the immune mechanisms occurring
during a continuous infusion and translating this model to study clinically
observed inflammation in sepsis patients.

## References

  * Akimoto et al., (2007) Akimoto, H., Yamada, A., & Takikawa, O. (2007).  Up-regulation of the brain indoleamine 2,3-dioxygenase activity in a mouse model of alzheimer’s disease by systemic endotoxin challenge.  Int Congress Series, 1304,357–361.
  * Bahador and Cross, (2007) Bahador, M. & Cross, A. (2007).  Review: From therapy to experimental model: A hundred years of endotoxin administration to human subjects.  J Endotoxin Res, 13,251–279.
  * Bangsgaard et al., (2017) Bangsgaard, E., Hjorth, P., Olufsen, M., Mehlsen, J., & Ottesen, J. (2017).  Integrated inflammatory stress (itis) model.  Bull Math Biol, 79,1487–1509.
  * Banks et al., (2009) Banks, H., Davidian, M., Samuels, J., & Sutton, K. (2009).  An inverse problem statistical methodology summary.  In Mathematical and Statistical Estimation Approaches in Epidemiology, pages 249–302. Springer, Dordrecht, Netherlands.
  * Berg et al., (2012) Berg, R., Plovsing, R., Ronit, A., Bailey, D., Holstein-Rathlou, N., & Møller, K. (2012).  Disassociation of static and dynamic cerebral autoregulatory performance in healthy volunteers after lipopolysaccharide infusion and in patients with sepsis.  Am J Physiol, 303,R1127–R1135.
  * Bickel, (1993) Bickel, M. (1993).  The role of interleukin-8 in inflammation and mechanisms of regulation.  J Periodontol, 64 (5 Suppl),456–460.
  * Brady, (2017) Brady, R. (2017).  Mathematical modeling of the acute inflammatory response & cardiovascular dynamics in young men.  PhD thesis, North Carolina State University.
  * Brady et al., (2018) Brady, R., Frank-Ito, D., Tran, H., Janum, S., Møller, K., Brix, S., Ottesen, J., Mehlsen, J., & Olufsen, M. (2018).  Personalized mathematical model of endotoxin-induced inflammatory responses in young men and associated changes in heart rate variability.  Math Model Nat Phenom, 13,42.
  * Brodin and Davis, (2017) Brodin, P. & Davis, M. (2017).  Human immune system variation.  Nat Rev Immunol, 17,21–29.
  * Burnham and Anderson, (2002) Burnham, K. & Anderson, D. (2002).  Model selection and multi-model inference: a practical information-theoretic approach.  Springer-Verlag.
  * Casey et al., (1993) Casey, L. C., Balk, R. A., & Bone, R. C. (1993).  Plasma cytokine and endotoxin levels correlate with survival in patients with the sepsis syndrome.  Ann Intern Med, 119(8),771–778.
  * Cecconi et al., (2018) Cecconi, M., Evans, L., Levy, M., & Rhodes, A. (2018).  Sepsis and septic shock.  The Lancet, 392,75–87.
  * Chow et al., (2005) Chow, C., Clermont, G., Kumar, R., Lagoa, C., Tawadrous, Z., Gallo, D., Betten, B., Bartels, J., Constantine, G., Fink, M., Billiar, T., & Vodovotz, Y. (2005).  The acute inflammatory response in diverse shock states.  Shock, 24,74–84.
  * Clodi et al., (2008) Clodi, M., Vila, G., Geyeregger, R., Riedl, M., Stulnig, T., Struck, J., Luger, T., & Luger, A. (2008).  Oxytocin alleviates the neuroendocrine and cytokine response to bacterial endotoxin in healthy men.  Am J Physiol, 295,E686–E691.
  * Copeland et al., (2005) Copeland, S., Warren, H., Lowry, S., Calvano, S., & Remick, D. (2005).  Inflammation and the host response to injury investigators. acute inflammatory response to endotoxin in mice and humans.  Clin Diagn Lab Immunol, 12(1),60–67.
  * Cunningham et al., (2005) Cunningham, C., Wilcockson, D., Campion, S., Lunnon, K., & Perry, V. (2005).  Central and systemic endotoxin challenges exacerbate the local inflammatory response and increase neuronal death during chronic neurodegeneration.  J Neurosci, 25,9275–9284.
  * Danz, (2023) Danz, A. (2023).  boxplotgroup.  MATLAB Central File Exchange.
  * Day et al., (2006) Day, J., Rubin, J., Vodovotz, Y., Chow, C., Reynolds, A., & Clermont, G. (2006).  A reduced mathematical model of the acute inflammatory response ii. capturing scenarios of repeated endotoxin administration.  J Theor Biol, 242,237–256.
  * Dobreva et al., (2021) Dobreva, A., Brady-Nicholls, R., Larripa, K., Puelz, C., Mehlsen, J., & Olufsen, M. (2021).  A physiological model of the inflammatory-thermal-pain-cardiovascular interactions during an endotoxin challenge.  J Physiol, 599,1459–1485.
  * Dodge, (2008) Dodge, Y. (2008).  The concise encyclopedia of statistics.  Springer Science & Business Media.
  * Easson et al., (1998) Easson, A., Bode, B., Fischer, C., & Souba, W. (1998).  Effects of endotoxin challenge on hepatic amino acid transport during cancer.  J Surg Res, 77,29–34.
  * Fitzal et al., (2003) Fitzal, F., DeLano, F., Young, C., Rosario, H., Junger, W. G., & Schmid-Schönbein, G. (2003).  Pancreatic enzymes sustain systemic inflammation after an initial endotoxin challenge.  Surgery, 134,446–456.
  * Foteinou et al., (2009) Foteinou, P., Calvano, S., Lowry, S., & Androulakis, I. (2009).  Modeling endotoxin-induced systemic inflammation using an indirect response approach.  Math Biosci, 217,27–42.
  * Foteinou et al., (2011) Foteinou, P., Calvano, S., Lowry, S., & Androulakis, I. (2011).  A physiological model for autonomic heart rate regulation in human endotoxemia.  Shock, 35,229–239.
  * Fredriksson et al., (2009) Fredriksson, K., Fläring, U., Guillet, C., Wernerman, J., & Rooyackers, O. (2009).  Muscle mitochondrial activity increases rapidly after an endotoxin challenge in human volunteers.  Acta Anaesthesiol Scand, 53,299–304.
  * Givalois et al., (1994) Givalois, L., Dornand, J., Mekaouche, M., Solier, M., Bristow, A., Ixart, G., Siaud, P., Assenmacher, I., & Barbanel, G. (1994).  Temporal cascade of plasma level surges in acth, corticosterone, and cytokines in endotoxin-challenged rats.  Am J Physiol, 267,R164–R170.
  * Heine et al., (2001) Heine, H., Rietschel, E. T., & Ulmer, A. (2001).  The biology of endotoxin.  Mol Biotechnol, 19,279–296.
  * Hjemdahl et al., (2011) Hjemdahl, P., Rosengren, A., & Steptoe, A. (2011).  Stress and cardiovascular disease.  Springer, London, UK.
  * Ho et al., (2013) Ho, T., Clermont, G., & Parker, R. (2013).  A model of neutrophil dynamics in response to inflammatory and cancer chemotherapy challenges.  Comput Chem Eng, 51,187–196.
  * Janum et al., (2016) Janum, S., Nielsen, S., Werner, M., Mehlsen, J., Kehlet, H., & Møller, K. (2016).  Pain perception in healthy volunteers: Effect of repeated exposure to experimental systemic inflammation.  Innate Immun, 22,546–556.
  * Jin et al., (2013) Jin, J., Han, X., & Yu, Q. (2013).  Interleukin-6 induces the generation of il-10-producing tr1 cells and suppresses autoimmune tissue inflammation.  J Autoimmun, 40,28–44.
  * Johnston and Webster, (2009) Johnston, G. & Webster, N. (2009).  Cytokines and the immunomodulatory function of the vagus nerve.  Br J Anaesth, 102,453–462.
  * Kashiouris et al., (2020) Kashiouris, M., L’heureux, M., Cable, C., Fisher, B., Leichtle, S., & Fowler, A. (2020).  The emerging role of vitamin c as a treatment for sepsis.  Nutrients, 12,1–16.
  * Kenney and Ganta, (2014) Kenney, M. & Ganta, C. (2014).  Autonomic nervous system and immune system interactions.  Compr Physiol, 4,1177.
  * Kiers et al., (2017) Kiers, D., Koch, R., Hamers, L., Gerretsen, J., Thijs, E., Van Ede, L., Riksen, N., Kox, M., & Pickkers, P. (2017).  aracterization of a model of systemic inflammation in humans in vivo elicited by continuous infusion of endotoxin.  Sci Rep, 7,1–10.
  * Krabbe et al., (2001) Krabbe, K., Bruunsgaard, H., Hansen, C., Møller, K., Fonsmark, L., Qvist, J., Madsen, P., Kronborg, G., Andersen, H., Skinhøj, P., & Pedersen, B. (2001).  Ageing is associated with a prolonged fever response in human endotoxemia.  Clin Diagn Lab Immunol, 8,333–338.
  * Kucharzik et al., (1998) Kucharzik, T., Lügering, N., Pauels, H., Domschke, W., & Stoll, R. (1998).  Il-4, il-10 and il-13 down-regulate monocyte-chemoattracting protein-1 (mcp-1) production in activated intestinal epithelial cells.  Clin Exp Immunol, 111,152–157.
  * Kumar et al., (2004) Kumar, R., Clermont, G., Vodovotz, Y., & Chow, C. (2004).  The dynamics of acute inflammation.  J Theor Biol, 230,145–155.
  * Leijte et al., (2019) Leijte, G., Kiers, D., Van Der Heijden, W., Jansen, A., Gerretsen, J., Boerrigter, V., Netea, M., Kox, M., & Pickkers, P. (2019).  Treatment with acetylsalicylic acid reverses endotoxin tolerance in humans in vivo: A randomized placebo-controlled study.  Crit Care Med, 47,508–516.
  * Lorenz et al., (2013) Lorenz, W., Buhrmann, C., Mobasheri, A., Lueders, C., & Shakibaei, M. (2013).  Bacterial lipopolysaccharides form procollagen-endotoxin complexes that trigger cartilage inflammation and degeneration: Implications for the development of rheumatoid arthritis.  Arthritis Res Ther, 15,R111.
  * Malek et al., (2015) Malek, H., Ebadzadeh, M., Safabakhsh, R., Razavi, A., & Zaringhalam, J. (2015).  Dynamics of the hpa axis and inflammatory cytokines: Insights from mathematical modeling.  Comput Biol Med, 67,1–12.
  * Merrill et al., (2011) Merrill, J., You, J., Constable, C., Leeman, S., & Amar, S. (2011).  Whole-body deletion of lps-induced tnf-$\alpha$ factor (litaf) markedly improves experimental endotoxic shock and inflammatory arthritis.  PNAS, 108,21247–21252.
  * Miao et al., (2011) Miao, H., Xia, X., Perelson, A. S., & Wu, H. (2011).  On identifiability of nonlinear ode models and applications in viral dynamics.  SIAM Rev, 53,3–39.
  * Miller et al., (2010) Miller, A., Pearce, B., Ruzek, M., & Biron, C. (2010).  Interactions between the hypothalamic-pituitary-adrenal axis and immune system during viral infection: pathways for environmental effects on disease expression.  In Handbook of Physiology, The Endocrine System, Coping with the Environment: Neural and Endocrine Mechanisms, volume Supplement 23, pages 425–450. Comprehensive Physiology.
  * Miller et al., (1979) Miller, M., Samuelson, C., Hiramoto, R., & Ward, J. (1979).  Endotoxin toxicity in rats with 6-sulfanilamidoindazole arthritis.  Infect Immun, 25,337–344.
  * Murphy, (2012) Murphy, K. (2012).  Janeway’s Immunobiology.  Garland Science, London, UK, 10th edition.
  * Nedeva et al., (2019) Nedeva, C., Menassa, J., & Puthalakath, H. (2019).  Sepsis: Inflammation is a necessary evil.  Front Cell Dev Biol, 7,1–12.
  * Olufsen and Ottesen, (2013) Olufsen, M. & Ottesen, T. (2013).  A practical approach to parameter estimation applied to model predicting heart rate regulation.  J Math Biol, 67,39–68.
  * Parker et al., (2016) Parker, R., Hogg, J., Roy, A., Kellum, J., Rimmelé, T., Daun-Gruhn, S., Fedorchak, M., Valenti, I., Federspiel, W., Rubin, J., Vodovotz, Y., Lagoa, C., & Clermont, G. (2016).  Modeling and hemofiltration treatment of acute inflammation.  Processes, 4,1–33.
  * Pope et al., (2009) Pope, S. R., Ellwein, L. M., Zapata, C. L., Novak, V., Kelley, C. T., & Olufsen, M. S. (2009).  Estimation and identification of parameters in a lumped cerebrovascular model.  Math Biosci Eng, 6(1),93–115.
  * Reynolds et al., (2006) Reynolds, A., Rubin, J., Clermont, G., Day, J., Vodovotz, Y., & Ermentrout, B. (2006).  A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation.  J Theor Biol, 242,220–236.
  * Rossol et al., (2011) Rossol, M., Heine, H., Meusch, U., Quandt, D., Klein, C., Sweet, M., & Hauschildt, S. (2011).  Lps-induced cytokine production in human monocytes and macrophages.  Crit Rev Immunol, 31,379–446.
  * Roy et al., (2007) Roy, A., Clermont, G., Daun, S., & Parker, R. (2007).  A mathematical model of acute inflammatory response to endotoxin challenge.  In AIChE Annual Meeting,, 538g, Salt Lake City, UT.
  * Scheff et al., (2010) Scheff, J., Calvano, S., Lowry, S., & Androulakis, I. (2010).  Modeling the influence of circadian rhythms on the acute inflammatory response.  J Theor Biol, 264,1068–1076.
  * Schwarz, (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  Ann Stat, pages 461–464.
  * Seber and Wild, (2003) Seber, G. & Wild, C. (2003).  Nonlinear regression.  Wiley-Interscience Paperback Series. Wiley, Hoboken, NJ.
  * Shinozaki et al., (2010) Shinozaki, S., Inoue, Y., Yang, W., Fukaya, M., Carter, E., Ming-Yu, Y., Fischman, A., Tompkins, R., & Kaneki, M. (2010).  Farnesyltransferase inhibitor improved survival following endotoxin challenge in mice.  Biochem Biophys Res Commun, 391,1459–1464.
  * Sly et al., (2001) Sly, L., Krzesicki, R., Brashler, J., Buhl, A., McKinley, D., Carter, D., & Chin, J. (2001).  Endogenous brain cytokine mrna and inflammatory responses to lipopolysaccharide are elevated in the tg2576 transgenic mouse model of alzheimer’s disease.  Brain Res Bull, 56,581–588.
  * Smith, (2013) Smith, R. (2013).  Uncertainty quantification: theory, implementation, and applications.  SIAM, Philadelphia, PA.
  * Su et al., (2009) Su, B., Zhou, W., Dorman, K., & Jones, D. (2009).  Mathematical modelling of immune response in tissues.  Comput Math Methods Med, 10,9–38.
  * Suffredini and Noveck, (2014) Suffredini, A. & Noveck, R. (2014).  Human endotoxin administration as an experimental model in drug development.  Clin Pharmacol Ther, 96,418–422.
  * Taudorf et al., (2007) Taudorf, S., Krabbe, K., Berg, R., Pedersen, B., & Møller, K. (2007).  Human models of low-grade inflammation: Bolus versus continuous infusion of endotoxin.  Clin Vaccine Immunol, 14,250–255.
  * Tilg et al., (1994) Tilg, H., Trehu, E., Atkins, M., Dinarello, C., & Mier, J. (1994).  Interleukin-6 (il-6) as an anti-inflammatory cytokine: induction of circulating il-1 receptor antagonist and soluble tumor necrosis factor receptor p55.  Blood, 83,113–118.
  * Torres et al., (2022) Torres, L., Pickkers, P., & van der Poll, T. (2022).  Sepsis-induced immunosuppression.  Annu Rev Physiol, 84,157–181.
  * Torres et al., (2019) Torres, M., Wang, J., Yannie, P., Ghosh, S., Segal, R., & Reynolds, A. (2019).  Identifying important parameters in the inflammatory process with a mathematical model of immune cell influx and macrophage polarization.  PLoS Comput Biol, 15,1–27.
  * van Lier et al., (2019) van Lier, D., Geven, C., Leijte, G., & Pickkers, P. (2019).  Experimental human endotoxemia as a model of systemic inflammation.  Biochimie, 159,99–106.
  * Varela et al., (2018) Varela, M., Mogildea, M., Moreno, I., & Lopes, A. (2018).  Acute inflammation and metabolism.  Inflammation, 41,1115–1127.
  * Verboogen et al., (2019) Verboogen, D. R., Revelo, N. H., Ter Beest, M., & van den Bogaart, G. (2019).  Interleukin-6 secretion is limited by self-signaling in endosomes.  J Mol Cell Biol, 11,144–157.
  * Wald et al., (2022) Wald, E., Badke, C., Hintz, L., Spewak, M., & Sanchez-Pinto, L. (2022).  Vitamin therapy in sepsis.  Pediatr Res, 91,328–336.
  * West and Heagy, (2002) West, M. A. & Heagy, W. (2002).  Endotoxin tolerance: a review.  Crit Care Med, 30(1),S64–S73.
  * Windoloski et al., (2023) Windoloski, K., Bangsgaard, E., Dobreva, A., Ottesen, J., & Olufsen, M. (2023).  A unified computational model for the human response to lipopolysaccharide-induced inflammation.  In Multiplicity of time scales in complex systems, Mathematics Online First Collections. Springer Cham, Switzerland.
  * Wu et al., (2009) Wu, H.-P., Chen, C.-K., Chung, K., Tseng, J.-C., Hua, C.-C., Liu, Y.-C., Chuang, D.-Y., & Yang, C.-H. (2009).  Serial cytokine levels in patients with severe sepsis.  Inflamm Res, 58,385–393.
  * Yassine, (2016) Yassine, F. (2016).  Effect of endotoxin challenge on normal, tumor initiated, and invasive human breast cells.  PhD thesis, Am University of Beirut, Libanon.
  * Zuckerman et al., (1996) Zuckerman, S., Ahmari, S., Bryan-Poole, N., Evans, G., Short, L., & Glasebrook, A. (1996).  Estriol: A potent regulator of tnf and il-6 expression in a murine model of endotoxemia.  Inflammation, 20,581–597.

## Additional Information

### Data availability

This study utilizes cytokine data from two previously published studies (Berg
et al., , 2012; Janum et al., , 2016). Time series data and mathematical
modeling computer simulations for the mean and each subject in both studies
are shown in the manuscript, and computer code for this study is available on
the GitHub repository: https://github.com/mjcolebank/CDG_NCSU.

### Competing interests

The authors declare that they have no competing interests.

### Author contributions

KAW, RMGB, and MSO were responsible for the design of the study and the
acquisition, analysis, or interpretation of data. SJ was responsible for the
acquisition of data. KAW and MSO were responsible for the mathematical model
calibration, simulation, and interpretations. KAW, SJ, RMGB, and MSO were
responsible for drafting or critically editing the manuscript. All authors
were responsible for approval of the final version of the manuscript. All
authors agree to be accountable for all aspects of the work. All persons
designated as authors qualify for authorship and all those who qualify for
authorship are listed.

### Funding

Ronan M.G. Berg and Susanne Janum were supported in part via funding to the
Centre for Physical Activity Research (CFAS) receiving suppot from TrygFonden
(grants ID 101390 and ID 20045).

### Acknowledgements

We would like to thank Jesper Mehlsen (Rigshospitalet, Denmark) for his
contribution to the acquisition of experimental data and Mitchel Colebank
(University of California, Irvine) for his guidance on uncertainty
quantification methods.

### Keywords

Inflammation, cytokines, mathematical modeling, data analysis, continuous
infusion, endotoxin challenge, administration method