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

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 used mathematical modelling to analyse cytokine data from two studies administering a 2 ng kg−1 dose of endotoxin, one as a bolus and the other as a continuous infusion over 4 h. Using our model, we simulated the dynamics of mean and subject‐specific cytokine responses as well as the response to long‐term endotoxin administration. Cytokine measurements revealed that the bolus injection led to significantly higher peaks for interleukin (IL)‐8, while IL‐10 reaches higher peaks during continuous administration. Moreover, the peak timing of all measured cytokines occurred later with continuous infusion. We identified three model parameters that significantly differed between the two administration methods. Monocyte activation of IL‐10 was greater during the continuous infusion, while tumour necrosis factor α and IL‐8 recovery rates were 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, the continuous infusion model exhibited prolonged inflammation with recurrent peaks resolving within 2 days during long‐term (20–32 h) endotoxin administration.

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-α, IL-1β, 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-α, 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-α, 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-α, 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-α 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.

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-α, IL-6, IL-8, and IL-10.
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.
We calculated the mean cytokine response after removing abnormal responses (outlying data points outside the 1.5×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 proand anti-inflammatory cytokines TNF-α 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-α, 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 differ-   and subject-specific continuous infusion and bolus data.

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.

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 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 where k M R (hr −1 ) denoted the regeneration rate and M ∞ (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-α at rate k M T N F (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 M A (hr −1 ).
In equations ( 2)-( 3), the upregulation or downregulation of state Y by state X was as where h represents the steepness of the curve and η Y X the half-maximum value.

Model summary:
The mathematical model described above was an ODE system of the form where X ∈ R 7 denoted the time-varying states X = {E, M R , M A , T NF, IL6, IL8, IL10} determining endotoxin (E), monocytes (resting M R and activated M A ), TNF-α, IL-6, IL-8, and IL-10 concentrations.The model parameters θ ∈ 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 where r k is the residual vector for each state k = {T NF, IL6, IL8, IL10} and 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 ); θ) denotes the model output for the state k at time t i , 1 ≤ i ≤ 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.

Nominal parameters
Nominal parameter values were taken from Brady (2017) except for cytokine baseline concentrations (w T N F , 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 = {T NF, IL6, IL8, IL10}, denoted X i , was scaled using the technique from Windoloski et al. (2023) where for the scaling factor α and desired peak concentration Xi .We substituted equation ( 11) into the ODE for state X i giving The scaling factor 1/α 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 Thus, half-saturation values were scaled by 1/α.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 T N F , 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 α < 0.9 or α > 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 χ was given by where y = g(t, X(t), θ) was the model output at time t, θ the nominal parameter set, and y data the mean continuous infusion data.We approximated the (i, j) entry in the submatrix χ k using forward differences, where For submatrix X k , elements χ ij were given by where φ = 10 −8 was the solver tolerance, h = √ φ the step size (Pope et al., 2009), and e j the basis vector in the jth direction.We ranked relative sensitivities by computing the twonorm of each column of χ, obtaining a single sensitivity per parameter.We repeated the sensitivity analysis by simulating 100 runs sampling parameters from a uniform distribution varying ±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 = χ T χ.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 where (i, j) refers to θ i and θ 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 χ = UΣV T where U and V contained the left and right singular vectors of χ, and Σ contained the singular values σ of χ.
The largest k singular values of χ were determined by σ(k) ≥ 10 √ φ where φ 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 T k 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 θ as θ = P T θ. (17) The first k parameters of θ 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 ±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 θ i where We denoted θi as the mean and σ i as the standard deviation of the estimated parameter θ i .Parameters in each subset with CoV ≥ 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 (R 2 ) (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 θi were computed as where N was the total number of data points, q was the number of parameters that were estimated, t α/2 N −q was the t-value from the student's t-distribution for confidence level 1 − α with N − q degrees of freedom, and the variance estimator matrix Σ was given by We defined χ similarly to equations ( 14) and (15) where and the (i, j) element of submatrix χ k ( θ) with k ∈ {T NF, IL6, IL8, IL10} was approximated using forward differences given by Here, ỹk i = g(t i , X k (t i ), θ) was the optimal model output with optimal parameter vector θ for cytokine state k at time t i for 1 ≤ i ≤ 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 where σ k was a diagonal matrix of size N k × N k with entries with y k data defined in equation ( 10).The asymptotic prediction interval for cytokine k at time t i was given by and the confidence interval by We defined ỹk i and N k as in equation ( 22), q k was the number of estimated parameters that impacted cytokine state k, and t α/2 N k −q k was the t-value for confidence level 1 − α with N k − q k degrees of freedom.χ k ( θ) was given by ( 21) and ( 22), but columns in χ T N F ( θ), χ IL6 ( θ), and χ IL10 ( θ) 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 singular unless removed.The matrix G T ik was defined as which was ith row of the submatrix χ k ( θ), and the variance estimator s 2 k was given by 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 (α = 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 (α = 0.05) was also performed on the set of optimized parameters from subjectspecific 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.

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-α (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-α, IL-6, and IL-8 (p < 0.0001) and IL-10 (p = 0.00695).

Sensitivity analysis and subset selection
Single and repeated sensitivity analysis (Figure 4) highlights the system's dependence on endotoxin, activated monocytes, TNF-α, and IL-10 dynamics.The system's most sensitive parameters are the growth or decay of these states, where k T N F M (growth rate of TNF-α by monocytes) and k M A (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-α and IL-10 are the main cytokines that upregulate and downregulate   other states.The least sensitive rate constant is k M R , 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.
The parameter k M R 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 The third subset is found by performing SVD-QR followed by the SCM on all 15 rate constants, which results in the subset 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 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 (R 2 , 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 R 2 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 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-α and IL-10 have larger peak concentrations.The immune resolution time dur-  We generated N = 61 data points for each cytokine to determine the mean data and model uncertainty using confidence level (1 − α) with α = 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  and 7b) is sufficiently robust to capture variation in data.This is evidenced by high R   22.8 (13.5) 22.9 (11.2)  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 9B Scaled Parameter Values

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-α 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-α, IL-6, and IL-8 with a bolus injection plus continuous infusion.Our study exhibited an increased TNF-α 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-α and IL-6 occurred significantly later with the continuous infusion compared to a bolus injection and (ii) TNF-α 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 antiinflammatory 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-α 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-α 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α, 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 36hour 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β and TGF-β, and signaling pathways such as the NF-κ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β was measured in the bolus study, its concentrations were not detectable.Our previous work (Dobreva et al., 2021)

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

Figure 1 :
Figure1: 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 fromDanz (2023).
ence ( d) between the mean bolus ( b) and continuous infusion (c) baseline dj 0 = bj 0 − cj 0 for cytokine j = {T NF, IL6, IL8, IL10} and then shifting the concentrations by bji (k) = b j i (k) − dj 0 ,where b j i (k) is the original cytokine concentration j at time i for the kth bolus participant (1 ≤ k ≤ 20).The adjusted cytokine concentration is denoted by bj i (k).

Figure 2 :
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).

Figure 4 :
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.

Figure 5 :
Figure 5: Optimized model fit to mean continuous infusion and bolus data estimating S F inal .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.

Figure 6 :
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.

Figure 7 :
Figure 7: Optimal model responses for (A) subject 1 from the continuous infusion study and (B) subject 16 from the bolus study.

Figure 8 :
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 (α = 0.05) include k T N F (p < 0.0001), k 8 (p < 0.0001), and k 10M (p = 0.0142).Parameters not statistically significant include k M A (p = 0.465), k T N F M (p = 0.106), and k 8M (p = 0.0615).This figure is generated using MATLAB code adapted from Danz (2023).
explored interactions of immune, cardiovascular, thermal, and pain responses during a bolus endotoxin challenge, andWindoloski 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.
enhance understanding of potential mechanisms impacting immune responses to endotoxin, we devised a physiologically-based mathematical model simulating mean and subjectspecific 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-α 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.

Table 1 :
Berg et al. (2012))haracteristics.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 fromJanum et al. (2016), and m = 9 for the continuous infusion data fromBerg et al. (2012).All concentrations were rounded to the nearest whole number.

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.

Table 5 :
2values Optimal 95% parameter confidence bounds for the mean continuous infusion and bolus model.

Table 6 :
Subject-specific parameter values as mean (SD) and estimated parameter p-values with significance level α = 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.
displays the model response for a 4-hour continuous endotoxin infusion of 2,