- 1Center for Research and Education in Aging, University of California, Berkeley, Berkeley, CA, United States
- 2Lawrence Berkeley National Laboratory, Berkeley, CA, United States
- 3Department of Integrative Biology, University of California, Berkeley, Berkeley, CA, United States
- 4Department of Computer Science, Federal University of Juiz de Fora, Juiz de Fora, Brazil
- 5Department of Mathematics, North Carolina State University, Raleigh, NC, United States
- 6College of Chemistry, University of California, Berkeley, Berkeley, CA, United States
- 7Mayo Clinic Alix School of Medicine, Scottsdale, AZ, United States
- 8Department of Molecular and Cell Biology, University of California, Berkeley, Berkeley, CA, United States
Cell-based mathematical models have previously been developed to simulate the immune system in response to pathogens. Mathematical modeling papers which study the human immune response to pathogens have predicted concentrations of a variety of cells, including activated and resting macrophages, plasma cells, and antibodies. This study aims to create a comprehensive mathematical model that can predict cytokine levels in response to a gram-positive bacterium, S. aureus by coupling previous models. To accomplish this, the cytokines Tumor Necrosis Factor Alpha (TNF-α), Interleukin 6 (IL-6), Interleukin 8 (IL-8), and Interleukin 10 (IL-10) are included to quantify the relationship between cytokine release from macrophages and the concentration of the pathogen, S. aureus, ex vivo. Partial differential equations (PDEs) are used to model cellular response and ordinary differential equations (ODEs) are used to model cytokine response, and interactions between both components produce a more robust and more complete systems-level understanding of immune activation. In the coupled cellular and cytokine model outlined in this paper, a low concentration of S. aureus is used to stimulate the measured cellular response and cytokine expression. Results show that our cellular activation and cytokine expression model characterizing septic conditions can predict ex vivo mechanisms in response to gram-negative and gram-positive bacteria. Our simulations provide new insights into how the human immune system responds to infections from different pathogens. Novel applications of these insights help in the development of more powerful tools and protocols in infection biology.
1 Introduction
The human immune system consists of the innate and adaptive immune response. Innate immunity comprises multiple lines of defense beginning with skin, saliva, and various secretions, and ending with non-specific leukocytes, while adaptive immunity refers to a long-term specific response initiated to eliminate a specific pathogen (Kim et al., 2007; Mirzaei et al., 2020). The immune system leukocytes are usually in the form of neutrophils, macrophages, eosinophils, or natural killer cells. Their primary function is to perform phagocytosis of pathogens and cell debris through engulfment and chemical degradation (Badwey and Karnovsky, 1980). If the innate immune system is insufficient in eliminating the pathogen, it activates the adaptive immune system, which is primarily composed of B cells and T cells, known for their specificity in function (Zhang and An, 2007). B cells produce and recruit antibodies that tag antigens on infected cells for T cells to subsequently destroy (Winer et al., 2011). In addition, both B and T cells play a role in the production of cytokines, which are small signaling proteins released by leukocytes that facilitate communication between immune cells.
One such signaling protein is TNF-α, which is instrumental to the acute phase reaction during an inflammatory response. In response to sepsis, TNF-α, a regulator of immune cells, functions by upregulating other cytokines, including IL-1 and IL-6 (Feuerstein et al., 1994). IL-6 is involved in inflammation and homeostatic processes, but though it is primarily a proinflammatory cytokine, it can also act as an anti-inflammatory cytokine through its inhibitory effects on TNF-α. During inflammation and the delay in apoptosis of T cells, IL-6 is critical for the recruitment of T cells and the production of B and T cells (Fournier and Philpott, 2005). Similarly, IL-8 is also involved in the recruitment of T cells, as well as basophils and neutrophils (Holmes et al., 1991). IL-8 has been shown to be induced by TNF-α and inhibited by IL-10 (Yao et al., 1996; Osawa et al., 2002) which is an anti-inflammatory cytokine critical in the regulation of immune responses (de Waal Malefyt et al., 1991). In particular, IL-10 limits the production of the proinflammatory cytokines such as IL-6, while downregulating the expression of TNF-α, T helper type 1 cytokines, and major histocompatibility complex class 2 molecules (Iyer and Cheng, 2012; Ip et al., 2017). Additionally, IL-10 acts as an immuno-regulator by maintaining homeostasis and preventing host damage during infection (Fournier and Philpott, 2005; Ip et al., 2017).
On a fundamental level, the innate immune system is able to recognize and respond to a wide range of triggers. To do so, it has evolved to include pattern recognition receptors (PRRs) that are able to recognize pathogen-associated molecular patterns (PAMPS) (Chandler and Ernst, 2017). PAMPS are features in microbes that are vital, common, and most importantly, conserved. Present in nearly all bacterial cell walls, peptidoglycan (PepG) is essential for maintaining cell structure and shape (Hergott et al., 2016). PepG protects cells from bursting due to turgor or in response to environmental stressors. The lipid components of cell membranes also contain distinguishing markers that are recognizable and stimulatory. Lipoteichoic acid (LTA), a hallmark of gram-positive bacteria, consists of a glycolipid covalently bound to a hydrophilic glycerophosphate polymer. The glycolipid component localizes in the lipid bilayer of cell membranes. Both PepG and LTA are recognized by the innate immune system and can trigger the systemic release of cytokines (Hergott et al., 2016). When injected into rats during an in vivo study, PepG and LTA acted in conjunction to trigger the release of TNF-α (De Kimpe et al., 1995).
A primary designation in the classification of gram-negative bacteria is the presence of lipopolysaccharide (LPS), instead of LTA, as part of the outer cell membrane. LPS, also known as endotoxin, is one of the most potent immunostimulants (Alexander and Rietschel, 2001) and consists of a glycan polymer, oligosaccharide core, and membrane-anchor lipid (Alexander and Rietschel, 2001; Cavaillon, 2018). Physiological recognition of the lipid component by the immune system causes pro-inflammatory cytokine activity (Alexander and Rietschel, 2001; Liu et al., 2018; Monguió-Tortajada et al., 2018; Tawfik et al., 2020). Endotoxin is the main driver of pro- and anti-inflammatory cytokine activation in gram-negative bacterial infections (Feezor et al., 2003; Tawfik et al., 2020). Endotoxin-induced inflammatory conditions cause similar levels of cytokine expression as those found in septic conditions due to a gram-positive bacteria such as S. aureus, one of the most common bacteria found on the surface of human skin (Foster, 2004; Tawfik et al., 2020). S. aureus, characterized by its thick, PepG layer and LTA layer, is known to stimulate the immune system through the release of toxins into the bloodstream (Feezor et al., 2003; McNicholas et al., 2014). LTA and PepG are the primary sources of activation of cytokines in response to gram-positive bacterial infections (Feuerstein et al., 1994; De Kimpe et al., 1995; Lowings et al., 2009; Cole et al., 2014) whose levels can be predicted using mathematical models.
Mathematical modeling of complex biological systems holds the potential for elucidating emergent properties of intricate biological pathways within the human body (DiLeo et al., 2009; Möbius and Laan, 2015; Morel et al., 2017; Caudill and Lynch, 2018; Meier-Schellersheim et al., 2019). In particular, immune system research can benefit from in silico simulations of drug and pathogen responses, which provide a deeper understanding of system dynamics that can be applied to design better diagnostic and treatment protocols (Chow et al., 2005; Eftimie et al., 2016). In general, in silico experimentation provides an alternative to tests that are difficult, impractical, expensive, or potentially unethical to perform in vivo (Winslow et al., 2012; Nijhout et al., 2015). Additionally, they provide a safer and more cost-effective platform for clinical drug testing, as potential drug candidates can be “pre-screened” (Winslow et al., 2012) Previous simulations successfully modeled activated and resting macrophages, plasma cells, antibodies, helper T cells, T-lymphocytes, and B-lymphocytes in response to the pathogen S. aureus (Quintela et al., 2014; Alvarez et al., 2019). However, they failed to account for the interactions between the cellular and cytokine response (Altan-Bonnet and Mukherjee, 2019; Torres et al., 2019; Du and Yuan, 2020).
The mathematical model presented herein, extends previous studies by taking into account the interconnectivity between cellular and cytokine responses and captures ex vivo and potential in vivo dynamics of the immune response (Quintela et al., 2014; Brady et al., 2016). This study focuses on examining the immune response to acute stimulation by a pathogen, which leads to a cascade of cellular signals recruiting leukocytes, or white blood cells, throughout the body (Hoebe et al., 2004; Anderson et al., 2019; Xue and Falcon, 2019). The cytokines of interest include TNF-α, IL-6, IL-8, and IL-10 as they are crucial in downstream signaling pathways that affect inflammatory and other responses within the immune system (Khoa et al., 2003). Additionally, several pharmaceutical agents can influence these cytokines, including Fosfomycin (FOM) (Michalopoulos et al., 2011), Clarithromycin (CAM) (PubMedHealth, 2018), Dexamethasone (DEX), and glucocorticoids (Johnson et al., 2021). FOM is an antibiotic that disrupts the cell walls of bacteria, specifically the gram-positive bacteria S. aureus, which blocks the production of IL-8 and amplifies the synthesis of IL-6 and IL-10 (Morikawa et al., 1996; Lan et al., 2018; Zheng et al., 2020). CAM is a macrolide antibiotic that diminishes the spread of bacteria while increasing the production of IL-10 (Morikawa et al., 1996). DEX is a steroid that hinders the immune system in the presence of inflammation (Johnson et al., 2021) by reducing the production of both IL-6 and IL-10 (Morikawa et al., 1996).
The mathematical model representing the basic pathway of the standard human immune response to S. aureus (Quintela et al., 2014) was coupled with a mathematical model predicting the cytokine response to LPS (Wang et al., 2000; Brady et al., 2016). Although the two models that were coupled represented responses to different pathogens (gram-positive vs gram-negative bacteria), in both cases the cytokine concentrations are calculated as functions of activated macrophages. Thus, the fundamental interactions between macrophages and cytokines are similar. Several literature sources exploring the immune response to PepG, LTA, and S. aureus in conjunction with LPS suggest that our assumption involving both antigens is reasonable (Fournier and Philpott, 2005; Zhang and An, 2007). S. aureus was chosen because it is an opportunistic pathogen, and notably is one of the leading causes of life-threatening infections including sepsis (Kwiecinski and Horswill, 2020).
The following mathematical model is formulated using a system of ordinary differential equations (ODEs) and partial differential equations (PDEs) emulating the intricate relationships among the pathogen, cytokines, and cells within the human immune system. As described in detail below, PDEs are used to represent the spatial features of the cellular model, while ODEs are used to model the cytokine activation dynamics. This model is used to simulate the cellular and cytokine activation induced by a low dose of S. aureus. Results show the proposed cellular-cytokine mathematical model can be used to predict ex vivo and in vivo experimental data induced by a given pathogen.
2 Methods
2.1 Experimental Data
Our coupled mathematical models are validated using data from two studies. The first measured cytokine concentrations from whole human blood induced by S. aureus PAMPs in healthy males and the second measured cytokine expression from human endothelial cells removed and induced by a low dose of S. aureus. Experiment descriptions are further referenced in Table 1.
2.1.1 Human Whole Blood Ex Vivo Response to S. aureus PAMPs
We first incorporated data from ex vivo studies predicting the response to S. aureus PAMPs. This study used whole blood data collected from healthy volunteers (Wang et al., 2000). The PepG was isolated from S. aureus, purified with hydrofluoric acid, and stored at -20°C. Venous blood from healthy subjects was anticoagulated with Na-citrate and analyzed after removal from incubation in the presence or absence of either 10 μg of PepG or 100 μg of LTA per ml of blood at the 1, 3, 6, 12, and 24-hour mark. The plasma removed during those time periods was centrifuged at 7,000 × g for 2 min and stored at -20°C for analyses using the enzyme immunoassay specific to TNF-α, IL-6, and IL-10 (Wang et al., 2000).
2.1.2 Human Endothelial Cells Ex Vivo Response to S. aureus
In the second study, endothelial cells (EC) from the human umbilical vein were removed and kept in a 5.5% CO2 tissue culture. The ECs were analyzed after removal from incubation in the presence of 108 CFU of S. aureus per ml at the 1, 3, 6, 12, and 24-hour mark. The ECs were removed from incubation and centrifuged for 30 min for analysis using enzyme-linked immunosorbent assay (ELISA) to measure IL-8 protein levels in the infected samples (Yao et al., 1996).
The experimental data were collected to formulate cellular-cytokine mathematical models. However, the ex vivo experimental data for each cytokine were measured in different units. The whole blood experiments for TNF-α, IL-6, and IL-10 were measured in relative concentrations whereas the endothelial cells from an umbilical vein for IL-8 were measured in ng/ml. The endotoxin cytokine-based model experiments for all aforementioned cytokines were measured in pg/ml. As a result, a conversion factor was implemented to compare the cytokines both to each other and to the simulation model results. This conversion factor used a relative unit metric to the peak value of the 24-hour data. A time-over-time evaluation of the minimum value was divided by the difference between the maximum and minimum values of that instance.
To ensure proper coupling, the cytokine mathematical expressions were initially derived from cytokine concentrations in healthy males induced with low doses of LPS.
2.1.3 Human Cytokine levels In Vivo Response to LPS
Cytokine data were extracted from previous studies, approved by the Regional Committee on Health Research Ethics (protocol ID H-3-2012-011) and the Regional Data Monitoring board (ID j-2007-58-0015, local 30-0766), and reported to clinicaltrials.gov (NCT01592526) (Brady et al., 2016). The participants were recruited by means of public advertising in Copenhagen, Denmark, and were required to meet specific safety protocols for the study. The study used data from 20 healthy adult males aged 18-35 years. The inclusion criteria were male, age 18-35, good general health, body mass index < 30 kg/m2, and written informed consent to participate in the study. The exclusion criteria were daily medicine intake, smoking, allergic reaction to nicotine, and splenectomy.
Each participant was under EKG signal supervision during the entire experiment. The experiment was initiated with an injection of a low dose (2 ng/kg body weight) of LPS. Blood samples were taken before the LPS injections, then at t=2, 3, 3.5, 4, 5, 6, 7, and 8 hours. The experimental protocol was repeated twice, and to reduce the risk of tolerance towards the endotoxin, the two days of experimentation were spaced 4 weeks apart. Following each experimentation day, samples were collected in EDTA tubes (Greiner bio-one, Germany), centrifuged at 4⁰C at 3500 rpm, then analyzed for cytokine concentrations using ELISA (Brady et al., 2016). The data from the participants were then analyzed to understand the effects of changes in cytokine expression. This in vivo data measuring the cytokine-cytokine response was used to test the differential equation models (Quintela et al., 2014; Brady et al., 2016).
2.2 Mathematical Model
A novel cellular-cytokine mathematical model (shown in Figure 1) combining the cellular model by Quintela et al. and the inflammatory cytokine model by Brady et al. is proposed (Quintela et al., 2014; Brady et al., 2016). The cellular model by Quintela et al. predicts the activation of the acquired immune response by activated macrophages acting as antigen presenting cells to the bacterium S. aureus, present in a portion of lung tissue. The prediction of resting and activated macrophages is used as an input to the novel cellular-cytokine model predicting the cytokine response to a pathogen. The novel cellular-cytokine model described by a system of partial (PDEs) and ordinary (ODEs) differential equations is solved over a 24-hour period. The cellular model is first solved in C++ and then used as input for the novel cellular-cytokine model, which ODEs are solved in MATLAB, and predicted states are compared to the experimental data for cellular and cytokine concentrations in response to bacterial infection with S. aureus.
Figure 1 A visual representation of the model components. (A) The simulation created by Quintela et al. representing the relationships between S. aureus and the different immune response cells were incorporated (Quintela et al., 2014). (B) A visual representation of the cytokine mathematical model as outlined by Brady et al. (2016). These relationships were used in the simulation via the relationships between the cytokines and the active macrophages rather than the pathogen itself (See Figure 2C). (C) The coupled model represented in this study. The green and red arrows indicate a positive (up-regulation) and negative (down-regulation) response, respectively. The various shapes indicate the parts of the simulation from Quintela et al. while the dark blue boxes indicate the parts of the simulation from the mathematical equations from Brady et al. The concentrations of the cytokines are solely dependent on the macrophage concentration and are not directly affected by the pathogen itself.
As noted in the introduction, numerous research studies have concluded that the host response to gram-positive and gram-negative bacteria provides similar activation of pro- and anti-inflammatory cytokines (Feezor et al., 2003; Fournier and Philpott, 2005). In an ex vivo study of 52 healthy patients, Feezor et al. concluded that activation of TNF-⍺, IL-8, and IL-10 induced by LPS and S. aureus show no significant difference in cellular activation and cytokine expression. This is further supported in studies by De Kimpe et al., 1995; Feezor et al., 2003; Fournier and Philpott, 2005, finding that LTA and PepG promoting activation of pro-inflammatory, anti-inflammatory, and chemoattractant properties have similar effects as a host response to LPS. As a result, the combination of LPS from gram-negative bacteria and PAMPs from gram-positive bacteria is reasonable due to their ability to induce similar inflammatory responses (Fan et al., 2007).
The widespread inflammation and septic shock caused by gram-positive bacteria, such as S. aureus, is primarily due to the function of LTA and PepG during an inflammatory response (Feezor et al., 2003; Moreillon and Majcherczyk, 2003; Fournier and Philpott, 2005). LTA and PepG work in conjunction to induce cytokine expression in the host’s innate and adaptive immune response to S. aureus. Mechanisms of action in the host including phagocytosis, neutrophil flux, and Sbi protein activation differ between the wall components and the pathogen itself, but they elicit similar inflammatory cytokine responses (Feezor et al., 2003; Moreillon and Majcherczyk, 2003; Fournier and Philpott, 2005). Results from these studies promote our assumption that the inflammatory response to the gram-positive bacteria S. aureus, the gram-positive bacteria cell wall constituents LTA and PepG, and the gram-negative bacteria wall constituent LPS, induce similar cytokine expression ex vivo (Feezor et al., 2003; Moreillon and Majcherczyk, 2003; Fournier and Philpott, 2005) but differ in their mechanism of action and host recognition.
These studies serve as a foundation for our coupled mathematical model combining the cellular PDEs and the cytokine ODEs via a shared variable, representing the concentration of activated macrophages averaged over the inflamed tissue.
2.2.1 Cellular Model
The cellular model by Quintela et al. predicts the activation of the acquired immune response to s. aureus as a function of the concentrations of bacteria (A), resting (M) and activated macrophages (M), and antibodies (F) spatially distributed in 1cm3 of lung tissue (x = (x,y,z)) as well as concentrations of lymphocytes (T,B), plasma cells (p), and antibodies (F) varying in time at the nearest lymph node. In this study, we assume that the the activated macrophages act as antigen presenting cells and migrate to the nearest lymph node where the specific response is triggered. Specific antibodies then travel to the infection site to opsonize the antigen. The activated macrophages are modelled both as spatially interacting with the antigen in the tissue (MA(x,t)) and at the lymph node interacting with the lymphocytes with concentrations varying only in time. The model is formulated using a system of partial and ordinary differential equations with the following dependent variables:
Spatial variables (pg/mm3):
- S. aureus bacteria (A(x,t))
- Resting Macrophages (MR(x,t))
- Activated Macrophages (MA(x,t))
- Specific Antibodies (F(x,t))
Temporal variables (pg/mm3):
Average Activated Macrophages
- T-Lymphocytes (T(t))
- B-Lymphocytes (B(t))
- Plasma cells (P(t))
- Antibodies (FL(t))
The details of the cellular model equations and the coupling from the tissue and nearest lymph node are available in Quintela et al. (2014). The novel cellular-cytokine model that is based both on the macrophage activation part of this cellular model and on the inflammatory cytokine model by Brady et al, is described below.
2.2.2 Novel Cellular-Cytokine Model
The cellular and cytokine dynamics are described using ODEs. We predict concentrations of TNF-α, IL-6, IL-8, and IL-10 as a function of the resting (MR(t)) and activated macrophages (MA(t)) considering the average tissue concentration obtained from the cellular model as initial condition. This model was originally developed to study the response of cytokines to LPS, but since LPS and PAMPs induce similar pro- and anti-inflammatory responses, we couple the cytokine model to the cell-based model described in Section 2.2.1. This is supported by findings from Fan et al. who discovered that Gi proteins present in both gram-negative and gram-positive bacteria contribute to the regulation of several cytokines and chemokines in response to bacterial stimuli (Fan et al., 2007). In the cytokine model, the dependent variables include:
Temporal variables:
- Bacteria: A(t)
- Resting Macrophages: MR(t)
- Activated Macrophages: MA(t)
- Tumor Necrosis Factor α (pro-inflammatory): TNF(t)
- Interleukin-6 (pro-inflammatory): : IL6(t)
- Interleukin-8 (pro-inflammatory): : IL8(t)
- Interleukin-10 (anti-inflammatory): IL10(t)
For each cytokine, up-and down-regulation is modeled using sigmoidal functions given by
where X represents the cytokine inducing up-regulation (superscript U) or down-regulation (superscript D) of cytokine Y. The half-maximum value is represented by η. These sigmoidal functions are incorporated within all the cytokine equations to describe the relationships between each of the cytokines. Specific parameter values are given in Table 2.
Table 2 Parameters, values, and units for the variables in the partial and ordinary differential equations found in the simulation (6, 7).
The cellular portion of the model proposed herein, predicts concentrations (pg/mm3) of bacteria (A), resting (MR) and activated macrophages (MA) varying over time.
S. aureus bacterium (A(x, t)) growth rate and rate of decline is modeled as:
where the first term represents the logistic growth of the bacteria, the constant kA represents the carrying capacity, and βA the replication rate. The second term gives the natural decay rate of the bacteria in the absence of any immune system processes through the natural decay coefficient, μA. The third and fourth terms describe the phagocytosis of S. aureus through activated and resting macrophages, with the constants λMA and λMR representing the rate of decline caused by activated and resting macrophages, respectively.
Resting MR(t) macrophages response to the pathogen S. aureus are modeled as:
where the first term in Eq (3) represents the constant influx rate of the resting macrophage (μMR up to MR MAX and the second term represents macrophage activation at the rate γMA in response to the bacteria and activation rate kMTNF considering the influence of the cytokines TNF and IL-10. The initial condition is given by the average of resting macrophages in the tissue as an outcome of the cellular model simulation over 24h.
Activated macrophages MA(t) are modeled as:
where the first term in Eq (4) represents macrophage activation at the rate γMA and kMTNF considering the influence of the cytokines TNF and IL-10, and the last term represents the decay rate of the activated (μMA) macrophages. The initial condition is given by the average of activated macrophages in the tissue as an outcome of the cellular model simulation over 24h. The initial concentration of both resting and activated macrophages are constant and at the boundary of the tissue neither the resting or activated macrophage concentration change.
In response to macrophage activation, the Tumor Necrosis Factor alpha (TNF-α) dynamics can be modeled as:
where the first term represents the down-regulating interactions that the cytokines IL-6 and IL-10 have with TNF-⍺ growth (at rate kTNFM) mediated by the average concentration of activated macrophages. The second term describes the rate in which TNF-⍺ naturally decays over time. As noted in the equation, the rate of change of TNF-α depends on the activated macrophages, which is predicted from the cellular model. In the cellular model, MA depends on both x and t. This response is integrated here and described further in Section 2.2.3.
Interleukin 6 (IL-6) activation is modeled as:
where the first term represents the interactions between TNF-⍺ (upregulating) and IL-10 (downregulating) affecting IL-6 production at a rate (kIL6TNF), and an IL-6, which also induce auto-negative feedback. The second term represents the natural decay (at rate k_IL6) of IL-6 towards a resting level (qIL6).
Interleukin 8 (IL-8) activation is modeled as:
where the first term represents the interactions between the opposing effects of TNF-⍺ (upregulating) at a rate (kIL8TNF) and IL-10 (downregulating) at a rate (kIL8M) stimulating the growth of IL-8 at a rate proportional to the average concentration of activated macrophages production, while the second term represents the natural decay rate of IL-8.
Interleukin 10 (IL-10) activation is modeled as:
where the first term describes the up-regulation of IL-10 due to IL-6 (at a rate kIL10IL6) and average concentration of activated macrophages (at a rate kIL10M), while the second term describes the natural decay rate of IL-10.
2.2.3 Coupled Model Numerical Solution
We first solve the cellular spatiotemporal model (over a 24-hour period) predicting the average concentrations of resting and activated macrophages in response to a low dose of S. aureus in the tissue. It is assumed that pathogen and macrophage movements can be represented as diffusion according to Fick’s Law (Crank, 1975). Therefore, each PDEs include a diffusion term with a specific diffusivity coefficient (Di, = MR, MA) estimated by Quintela et al. based on experimental observations for the pathogen and the cells that are included in this study. The diffusion term representing the rate of transfer of cells from one site to another is proportional to their concentration gradient. For simplicity, it is assumed that the medium is isotropic and has the same diffusion coefficient for every direction (Quintela et al., 2014), i.e., we model diffusion via the term, DiΔMj, i = MR, MA and j= R, A denote the diffusion of the macrophages in the tissue, again Δ refer to the second order derivative in space. At the onset of the simulation there aren’t any activated macrophages in the tissue (MA0 = 0), i.e., the resting macrophages that are equally distributed over the domain at concentration MR0. The cellular model is solved in C++ using the finite differences method.
Following the simulation of the cellular model, we calculate the average number of resting and activated macrophages by integrating the resulting concentrations of each over the discretized 1 cubic cm domain as:
The average concentrations for resting , and activated , macrophages are used as initial condition for the novel cellular-cytokine model. The ODEs ((1) - (8)) are then solved in MATLAB using the ode45 function.
2.3 Model Parameters
The model parameters include diffusion coefficients and replication, decay, activation, phagocytic, and migration rates of bacteria and cellular macrophages from the cellular model, and cytokine rate constants, source terms, and half-maximum values exponents from the cytokine model. To simplify the complexity of the human immune response, the cellular model parameters, cytokine half-maximum value and Hill function exponent parameters are held constant. They can be adjusted for model specifications, including representing in vivo septic conditions.
Diffusion coefficients: Di denotes the diffusion of the particular bacteria or cell into the tissue. Values from all cellular model parameters were obtained from literature research (Marchuk, 1997; Sarah and Richard, 2009; Pigozzo et al., 2013).
DA - Bacteria diffusion coefficient
DMR - Resting macrophage diffusion coefficient
DMA - Activated macrophage diffusion coefficient
Replication, decay, activation, and phagocytic rates: The following parameters denote the different rates of growth and decay for the bacteria and cells.
βA - Replication rate of the bacteria
kA - Carrying capacity of the bacteria
μA - Natural decay rate of the bacteria
μMR - Natural decay rate of the resting macrophages
μMA - Natural decay rate of the activated macrophages
γMA - Rate in which resting macrophages become active
λMR - Activation of the macrophages
λMA - Destruction rate of the bacteria by activate macrophages
λAF│MR- Destruction rate of opsonized bacteria by resting macrophages
λAF│MA - Destruction rate of opsonized bacteria by activated macrophages
Cellular rate of migration: The αi denotes the migration rate of the macrophages to the site of infection.
αMA - Migration rate of activated macrophages
αMR - Migration rate of resting macrophages
Rate constants: ki denotes cytokine activation or elimination rates, and kij determines the rate of change in the upregulation rate of a cytokine secreted from activated macrophages. These parameters have initial values based on predicted conditions of our model activated with a low dose of S. aureus.
kMTNF – Activation rate of resting macrophages influenced by TNF-α
kTNF - Activation rate (per hour) of TNF-α
kTNFM- Upregulation of TNF-α by the activated macrophages
k6- Activation rate (per hour) of IL-6
k6M- Upregulation of IL-6 by the activated macrophages
k6TNF- Upregulation of IL-6 by TNF-α
k8- Activation rate (per hour) of IL-8
k8M- Upregulation of IL-8 by the activated macrophages
k8TNF - Upregulation of IL-8 by TNF-α
k10 - Activation rate (per hour) of IL-10
k10M - Upregulation of IL-10 by the activated macrophages
k106 - Upregulation of IL-10 by IL-6
Source parameters: qi represents the base concentration of each cytokine in the absence of pathogen stimulation. qi are also used to set initial conditions for each cytokine. These parameters values are established based on initial predicted conditions of the model in the absence of S. aureus.
qTNF - The concentration of TNF-α in the absence of a pathogen
qIL6 - The concentration of IL-6 in the absence of a pathogen
qIL8 - The concentration of IL-8 in the absence of a pathogen
qIL10 - The concentration of IL-10 in the absence of a pathogen
Half-maximum value: The ηi parameters describe the effector cytokine concentration at which target cytokine activity would reach half-maximum with units of pg mL-1. They are included in sigmoidal Hill functions used to model up- or down-regulation of a specific target cytokine by a specific effector cytokine.
ηTNF6 - Half-maximum value associated with downregulation of TNF-α by IL-6
ηTNF10 - Half-maximum value associated with downregulation of TNF-α by IL-10
η610 - Half-maximum value associated with downregulation of IL-6 by IL-10
η66- Half-maximum value associated with the auto-negative feedback of IL-6
η6TNF - Half-maximum value associated with upregulation of IL-6 by TNF-α
η810 - Half-maximum value associated with downregulation of IL-8 by IL-10
η8TNF - Half-maximum value associated with upregulation of IL-8 by TNF-α
η106- Half-maximum value associated with upregulation of IL-10 by IL-6
Hill function exponent: ℎi represents the steepness of the Hill functions used to model the up- or down-regulation of each interaction.
h106- Hill function exponent associated with upregulation of IL-10 by IL-6
h6TNF - Hill function exponent associated with upregulation of IL-6 by TNF-α
h66 - Hill function exponent associated with auto-negative feedback of IL-6
h610 - Hill function exponent associated with downregulation of IL-6 by IL-10
h8TNF - Hill function exponent associated with upregulation of IL-8 by TNF-α
h810- Hill function exponent associated with downregulation of IL-8 by IL-10
hTNF10- Hill function exponent associated with downregulation of TNF-α by IL-10
hTNF6 - Hill function exponent associated with downregulation of TNF-α by IL-6
Parameters for half-maximum (ηi) and the respective exponents (hi) are fixed at their nominal values given in Table 2, while the rate constants k and source parameters q are estimated to fit the model to data.
2.4 Sensitivity Analysis
We have performed a sensitivity analysis to estimate how the model solution is affected by small perturbations to each model parameter. The sensitivity index was defined as the ratio
where j denotes a model output that depends on a parameter p, δ is a perturbation to the parameter p, and δJ is the resulting perturbation to the output j. The sensitivity index is a measure of the percentage of change in the output given a perturbation in each parameter (Quintela et al., 2014). The value of each parameter was varied by 10%, while other parameters were kept fixed at their baseline values.
2.5 Parameter Estimation
Parameters for the cell model are taken from literature, whereas we estimate rate constants and source parameters for the cytokine model minimizing the least squared error (j) between model predictions and data, given by
where rX = {Xm(t1) – Xd(t1), Xm(t2) – Xd(t2), …, Xm(tN) – Xd(tN)}/,and X = TNF, IL6, IL8, IL10 Superscript (m) refer to the model prediction and superscript (d) to the data. Parameter estimation was done using the MATLAB fminsearch function.
2.6 Statistical Analysis
A regression model was used to compare simulated results to experimental data. In each regression, simulated relative concentrations of each cytokine were compared to experimental concentrations of the same cytokine, induced by either LTA or PepG.
We assume linear regression of the form,
Where c1 represents the y-intercept and c2 represents the slope of each linear least squares (LLS) regression. This analysis was completed in MATLAB using the linear least squares package fitlm and shows the significance of the simulation relative to the experimental data (Table 3).
Table 3 LLS regressions of simulated results for each cytokine vs. experimental data of immune response to LTA and PepG.
3 Results
The results illustrate the model simulation in comparison to experimental data. Similar to LPS, S. aureus causes an initial spike in the concentrations of the cytokines. As the bacteria multiply, cytokine concentrations increase, reaching a plateau at which the bacterial growth rate equals the immune response rate. At this stage, the immune response decreases until the response rate overcomes the growth of the bacteria (Brennan and Zheng, 2007).
To test the validity of the model on a macroscale, the model was used to simulate a response to an infection by S. aureus. As expected, there were clear indications of the immune system response. From the start of infection, simulated macrophage and cytokine concentrations increase in proportion to the activated macrophages until a maximum concentration is reached at five hours (Figure 2). Once the bacteria were eliminated by the immune system response in the simulation, the macrophages and cytokines started reverting back to their basal values (Figure 3).
Figure 3 S. aureus average cell concentration in the tissue with and without immune response over a 24-hour period.
The microscale validity of the simulation was tested through close comparison of the simulation to the studies of ex vivo stimulation conducted by Wang et al. and Yao et al. (Yao et al., 1996; Wang et al., 2000; Wang and Deisboeck, 2014). Cytokine response by LTA and PepG stimulation in whole blood samples from those studies closely resemble the simulated cytokine response described herein (Yao et al., 1996; Wang et al., 2000). The peaks of TNF-α, IL-10, IL-6, and IL-8 coincide between the model results and ex vivo blood sample data (Yao et al., 1996; Wang et al., 2000) (Figure 4).
Figure 4 Comparison of the simulated TNF-⍺ and experimental TNF-⍺ activity in response to the introduction of 10 μg of PepG/mL or 100 μg of LTA/mL of human blood over the 24-hour period (A). Comparison of the simulated IL-6 and experimental IL-6 activity based on the introduction of 10 μg of PepG/mL or 100 μg of LTA/mL of human blood over a 24-hour period (B). Comparison of the simulated IL-10 and experimental IL-10 activity based on the introduction of 10 microgram of PepG/mL or 100 μg of LTA/mL of human blood over the 24-hour period (C). Comparison of the simulated IL-8 and experimental IL-8 activity based on the introduction of S. aureus-infected endothelial cells containing 10 μg of PepG/mL and 100 μg of LTA/mL over the 24-hour period (D). All cytokine concentrations are relative values as discussed in the methods.
Due to the different magnitudes of cytokine concentrations in the coupled simulation and blood sample data, a conversion factor to convert concentrations to relative values was introduced for meaningful comparison (Wang et al., 2000). Due to the limitations in the whole blood samples, TNF-⍺, IL-6, IL-8, and IL-10 were studied in response to low concentrations of PepG and LTA. On the contrary, the coupled simulation studied the immune response to antigen concentrations at significantly higher concentrations. Thus, the initial cytokine, antigen, and macrophage concentrations were modified to match the experimental conditions for a more accurate comparison (Figure 5).
Figure 5 Cytokine concentrations over a 24-hour period in response to stimulations with low dose S. aureus in our model.
Additionally, sensitivity analyses were conducted to validate the ex vivo experimental data with the simulated results. Each parameter in the cytokine equations were tested through its real-time concentration changing fold number, acquired by dividing the estimated parameters by their initial concentration and integrated over a 24-hour period (Figures 4, 5). The sensitivity index of 50 parameters in the model following a 24-hour simulation was conducted by implementing 10% parameter variations in each cytokine (Figure 6, 7). The parameters qTNF and kTNF in the TNF-⍺ equation, qIL8 and k8 in the IL-8 equation, qIL6 in the IL-6 equation, and k10M in the IL-10 equation were determined to be the most sensitive and therefore suitable for use in fitting the conditions of the experimental data.
Figure 6 Parameter adjustments of the individual cytokines. A ten-fold increase and decrease in the TNF-⍺ parameters (A), IL-6 parameters (B), IL-8 parameters (C), and IL-10 parameters (D).
Figure 7 Sensitivity indices denoting the most influential parameters to each cytokine after 24h simulation. Shown are the first 10 parameters that influence at least a 10% change in the resulting value of at least one of the four analyzed variables. Negative sensitivity indexes indicate reduced cytokine output while the omitted parameters trivially affected cytokine output.
To assess the performance of the model, a linear least squares (LLS) regression on each cytokine was utilized to compare the simulated results to the LTA- and Pep-G based immune response studies. The results of the regressions and corresponding regression validation parameters are given in Table 3. An F-test was performed on each regression to test the fit of the linear regression model. Each F-test resulted in a p-value significant to the ⍺ = 0.05 level for all comparisons of simulated and experimental data for cytokine-inducer interactions except for the TNF-⍺ and LTA interaction (significant to ⍺ = 0.10 level). This inconsistency was attributed to an outlier in the experimental TNF-⍺ LTA data, leading to a regression with a lower coefficient of determination. Based on the output of the statistical model, the cellular and cytokine mathematical model formulated under septic conditions accurately predicts whole human blood ex vivo conditions. Further development and clinical trial data will allow for additional cross-validated statistical algorithms to be run by splitting subjects into training, validation, and test sets for regressive modeling.
The linear regression models are given with 95% confidence intervals with a slope of 1 indicative of a perfect fit between the simulated concentrations and experimental results. Within one standard deviation, all of the cytokines fit the model except for IL-6 induced by LTA. The root mean squared error (RMSE), coefficient of determination, R2, and the adjusted coefficient of determination were also calculated to estimate the error distribution and variability of each regression (Table 3). The results of the statistical analysis indicate significant model similarities to the ex vivo experimental results and thus validate the accuracy of the simulation through clinical data.
To further ensure the authenticity of the simulation in a biological setting, the spatial domains and diffusion within the 24-hour period were measured from the initial injection of the antigen in a small portion of tissue. The overall change in the levels of bacteria, resting macrophages, and activated macrophages are shown in Figures 8–10. This discretized domain of the 3-D diffusion model is a hexahedron representing 10 mm3 of tissue. This initial injection of antigen was represented in the center of the hexahedral domain of simulation (between 3 mm and 7 mm over the axes). Initially, the presence of macrophages is equally distributed over the tissue. Following their initial interaction with S. aureus or its cell wall components, the resting macrophages are activated, resulting in the production of cytokines and chemokines. Acting as antigen-presenting cells, activated macrophages travel to the nearest lymph node and present the antigen to lymphocytes, which initiate the activation and differentiation of T-lymphocytes and B-lymphocytes into T-helper 2 lymphocytes and plasma cells, respectively. All images show a cut view of the volume along the x-axis in order to better visualize both the initial condition and the concentration and diffusion of bacteria and macrophages (Figures 8–10).
Figure 8 Diffusion of the bacteria at different time periods. Initial condition of bacteria injected only at the center of the domain (A), after 12h of simulation (B), and after 24h of simulation (C).
Figure 9 Diffusion of the resting macrophages at different time periods. Initial condition of the resting macrophages (A), after 3h of simulation where the values decrease as a result of change of state from resting to activated (B), after 12h of simulation (C), and after 24h of simulation (D).
Figure 10 Diffusion of the activated macrophages at different time periods. Initial condition of the resting macrophages (A), after 3h of simulation where the values increase as a result of change of state from resting to activated (B), after 12h of simulation (C), and after 24h of simulation (D).
4 Discussion
In this study, we explored the cellular-cytokine relationships of TNF-⍺, IL-6, IL-8, and IL-10 in response to S. aureus by utilizing mathematical modeling to predict cytokine levels in silico, and clinical literature and statistical analysis to validate the results of the simulation.
Cellular-cytokine interaction. The cellular model by Quintela et al. outlines the relationships between S. aureus (A), activated macrophages (MA), resting macrophages (MR), and antibodies (F) while the cytokine model by Brady et al. combines activated and resting macrophages with TNF-⍺, IL-6, IL-8, and IL-10 (Quintela et al., 2014; Brady et al., 2016). This model seeks to combine these models by the association of activated macrophages.
Sensitivity analysis and parameter estimation. This simulation validates the clinical experimentation from Wang et al. and Yao et al. (Yao et al., 1996; Wang et al., 2000). To scale the results of the clinical results, the parameters for each ODE was analyzed and modified using a min search nonlinear optimization to yield meaningful comparison.
Model assumptions and limitations. The limitations on the accuracy and precision of the model are discussed. This simulation does not account for factors such as neutrophil flux, complement response, humoral immune response effects, or variation between individuals, among other factors too complex to mathematically model and validate using clinical data.
Future studies. While this model seeks to expand the depth of current immune system models, the model remains incomplete due to the limitations in mathematical immune system research. The future expansions of the model open new pathways for new immune system research and may facilitate large-scale in silico pharmaceutical testing.
4.1 Cellular-Cytokine Interaction
This model combines the cellular model by Quintela et al. and the cytokine model by Brady et al. (Quintela et al., 2014; Brady et al., 2016). The cell model predicts the relationships between S. aureus (A), activated macrophages (MA), resting macrophages (MR), and antibodies (F), while the cytokine model studies how changes in activated and resting macrophages impact cytokine dynamics for TNF-⍺, IL-6, IL-8, and IL-10 (Quintela et al., 2014; Brady et al., 2016).
To reconcile the difference between the cell model, predicted as a function of time and space, and the cytokine model, which only varies with time, the activated macrophages were integrated into average activated macrophages and implemented into the cytokine model. In order to ensure both models operate within the same time frame, the concentrations of the initial bacteria and activated macrophages reported by Brady et al. (Brady et al., 2016) were scaled to match those reported in the ex vivo studies by Wang et al. and Yao et al. (Yao et al., 1996; Wang et al., 2000). Results of this parameter scaling demonstrated that we were able to match the model to all data metrics with statistical significance (pvalues <0.05) (Table 3) (Yao et al., 1996; Wang et al., 2000).
For future studies, a similar approach can be used to adapt the model to other experimental settings, e.g. to study the effect of varying macrophage destruction, flux, and phagocytosis within specific individuals (Emam et al., 2019). While the current simulation does not have the ability to accurately predict individual cytokine response, it has the capacity to model an average response to S. aureus, and it holds the potential to have a more profound impact with continued expansion in future studies.
4.2 Sensitivity Analysis and Parameter Estimation
Sensitivity analysis is an interpretable and adaptable tool used to provide insight into computational immunology studies investigating different components of the immune system to understand the extent of the spatial-temporal variables and parameters at play. This procedure has been shown to provide insight in computational immunology studies investigating different components of the immune system and their activation in response to a pathogen (Chen et al., 2019; Faro et al., 2019).
This coupled simulation is the product of sensitivity analysis using parameter estimation and model fitting to ex vivo data from Wang et al. and Yao et al. (Yao et al., 1996; Wang et al., 2000). Results of this analysis indicated that k8m, which describes the upregulation of IL-8 by the activated macrophages, as the most influential parameter and was subsequently modified to mirror the clinical data. Throughout the analysis, the initial concentrations of S. aureus and the cytokines were scaled to yield meaningful comparison.
4.3 Model Assumptions and Limitations
While a comprehensive model would be powerful on a global scale, the current simulation is limited by the lack of clinically backed mathematical models of the human immune system.
The reconciliation of utilizing S. aureus (a gram-positive bacteria), from the cell model, and LPS (a component in gram-negative bacteria), used in the cytokine model, are discussed previously. The body’s response to the bacteria may change based on the damage inflicted by the pathogen including sepsis and septic shock. From a mathematical and clinical property standpoint, there are minimal differences in the activation of the immune system cells and cytokine expression between S. aureus and LPS, allowing for a simplified mathematical model coupling (Feezor et al., 2003; Fournier and Philpott, 2005). However, future studies investigating each PAMP independently would aid in further confirming or validating our findings, and additionally may need to incorporate the complementary effects of LTA and PepG in S. aureus and LPS in gram-negative bacteria to model the array of pro-inflammatory (TNF-α and IL-6), anti-inflammatory (IL-10), and neutrophil chemoattractant (IL-8) responses under septic conditions (De Kimpe et al., 1995; Fournier and Philpott, 2005).
While the simulation results accurately predict ex vivo data, it should be noted that several assumptions were needed throughout the design of the model. The concentrations of S. aureus, LTA, and PepG are assumed to be proportional at low doses based on their activation of cytokines in the human immune response. Furthermore, toll-like receptors (TLRs) found in the membrane and cytosol of macrophages are the primary sentinels of PAMPs, and we assume a proportional relationship between macrophage concentration and PAMP recognition. Our model considers circulatory immune elements and does not consider intercellular or genetic regulatory aspects of immune response. (De Kimpe et al., 1995; Fournier and Philpott, 2005; Pigozzo et al., 2013).
Although the self-regulatory cytokine network model can respond to higher concentrations of S. aureus (6), this simulation utilized low concentrations of S. aureus similar to those found in ex vivo experimental conditions. Higher concentrations of S. aureus cause rapid increases in cytokine and cellular responses due to tissue damage corresponding to the endotoxicity of gram-positive bacteria. However, these situations were omitted due to the low concentrations of S. aureus simulated in this work and insufficient data to validate increased concentrations.
The neutrophil flux into the tissue was also omitted from this simulation due to the lack of clinical experiments and presence of high variability between subjects (Spaan et al., 2013). S. aureus-neutrophil interactions are human-specific and may influence the way this model predicts average cytokine levels. Future simulations may incorporate the effects of neutrophils in the IL-8 ordinary differential equation as additional studies and ex vivo experiments focusing on their interactions emerge (Fournier and Philpott, 2005).
To simplify the complexity of the bacterial complement response, this model does not incorporate the relevant complement proteins. Complement response of the human body to pathogens such as S. aureus plays a role in the ability of the human body to activate chemoattractants for phagocytosis of the bacteria (Fournier and Philpott, 2005; Laarman et al., 2010). The complement response has variable cytokine effects on the human immune response based on the extent and type of inflammatory condition, which further complicates its role in cytokine production and mediation (Fournier and Philpott, 2005). Moreover, complement, cytokine, and chemokine responses have overlapping biological effects on the body under septic conditions, and are therefore omitted from this model (Charchaflieh et al., 2012).
The effect of S. aureus on the humoral immune response is another limitation of this study given the insufficient findings of how that mechanism can be modeled (Smith et al., 2011). The pathogen is known to suppress the humoral immune response by means of the protein Sbi and therefore can allow extended survival of the bacteria (Smith et al., 2011). The effects of the protein Sbi on the concentration of bacteria inputted in the system were omitted from this model to solidify the cellular-cytokine interactions within the cellular model.
Macrophage destruction and flux are used as constants as a generalization in this model. Although these terms are human-specific and vary in response to infection by S. aureus, this model simulates the average response. Moreover, many components of the human immune response vary across populations. This complex problem of variability presents a challenge to mathematical modeling and is not addressed in this model. Therefore, certain variabilities such as macrophage destruction and flux were generalized to an average value to simplify the complexity of a human model. Understanding the different sources and mechanisms of action allows simulations to be used as predictive models to limit the scope of assumptions. One possible solution to account for patient-specific variation may be through estimating parameters individually for each dataset, as discussed in detail in Brady et al. (2016).
While this simulation can predict the average cytokine response to S. aureus, it does not account for all of the individualized variabilities within each immune response. This model aims to simplify the cytokine dynamic and become a foundation for future immune system expansion with increased immune system research. The model remains confined by the limited number of relevant mathematical immune system models and the availability of experimental data for validating more complex models. Continuous in vivo and ex vivo research into the role of the human immune system in response to a pathogen will provide additional data for validation and expansion of this model.
4.4 Future Studies
The human immune system, similarly, to many other body systems, is adversely affected by the aging process. Studies have shown age-associated alterations in immune system mechanisms such as a decreased T cell activation, reduced neutrophil efficiency, and altered toll-like receptor 1 (TLR1) expression (Oh et al., 2019) limit the ability of the immune system to eliminate foreign pathogens by reducing surface expression and effectivity (Khoa et al., 2003). In particular, the body’s cytokine and chemokine levels change with age, leading to chronic inflammation and progression of other degenerative diseases (Ponnappan and Ponnappan, 2011). This improved model along with other mathematical models simulating immune system functions may facilitate a better understanding of inflammatory responses and mechanisms that lead to the degradation of the human immune system and may lead to effective solutions towards preserving the functions of the immune system during the aging process.
Mathematical and computer models of the complex relationships within the immune system can open larger avenues for pharmaceutical and biochemical research that can be used to combat immune system diseases (Wang and Deisboeck, 2014). By incorporating both ordinary and partial differential equations, this simulation is able to provide a more realistic representation of the complex relationships within the immune system and, with further expansion, may serve as a vehicle for drug testing in silico.
While not incorporated into this model, pharmaceutical drugs such as FOM (Michalopoulos et al., 2011), CAM, (PubMedHealth, 2018) and DEX (Johnson et al., 2021) may significantly alter the results of this simulation. FOM inhibits PepG production in the cell wall of gram-positive bacteria, such as S. aureus, while also enhancing the production of IL-6 and IL-10. This may alter the cytokine production rate in this model as decreased S. aureus production eliminates the need for other cytokines in the immune response. With FOM, IL-6 and IL-10 production may increase.
The antibiotic, CAM, down-regulates the production of TNF-⍺, IL-6, and IL-8, which would lead to a measurable decrease in cytokine production. This cytokine inhibitor can significantly affect the simulation results, causing downward shifts of multiple cytokine values. Lastly, DEX, a steroid that hinders the immune system in the presence of inflammation, reduces the production of both IL-6 and IL-10, resulting in reduced TNF-⍺, IL-6, IL-8, and IL-10 concentrations.
Future expansions of this model may provide a clearer image of the various interactions within the immune system, facilitate a better understanding of the mechanisms that lead to the degradation of the immune system during the aging process, and may become a vehicle for in silico clinical trials (Swain and Nikolich-Zugich, 2009). This model is foundational and future clinical research is encouraged to characterize additional cytokines and cell behavior to create a more comprehensive and accurate understanding of the immune system. In particular, IL-12 is a proinflammatory cytokine that forms a vital link between the innate and adaptive immune system and future work would benefit greatly from accurately modeling IL-12 behavior (Trinchieri, 2003).
Data Availability Statement
The original contributions presented in the study are included in the article. Further inquiries can be directed to the corresponding author.
Ethics Statement
The studies involving human participants were reviewed and approved by the Regional Committee on Health Research Ethics (protocol ID H-3-2012-011) and the Regional Data Monitoring board (ID j-2007-58-0015, local 30-0766), and reported to clinicaltrials.gov (NCT01592526). The patients/participants provided their written informed consent to participate in this study.
Author Contributions
KT contributed to the design, analysis, drafts, format, simulation, and overall production of the paper. SG contributed to the design and conception of the study. BQ, RS, and ML contributed to the mathematical modeling and coding aspects of the study. MO contributed to the mathematical modeling, writing and proofreading of the paper. JC contributed to the design, analysis and drafts of the paper. JJ contributed to the research collection, review, and major editing of the paper. PB contributed to the research collection and references of the paper. ES contributed to the research collection of the paper. WP contributed to the statistical analysis and simulation of the paper. NM contributed to the equation review and overall model parameters. MM contributed to the sensitivity analysis and PDE diffusion model. All authors contributed to the article and approved the submitted version.
Funding
This study was funded by the Center for Research and Education in Aging, UC Berkeley Endowment.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We gratefully acknowledge the efforts of the following: Judith Segall, Paul Segall, Biotime, The Center for Research and Education in Aging, Renee Brady [Moffitt Cancer Center], Hal Sternberg, Ivan Labat, and George A. Brooks [Department of Integrative Biology, University of California, Berkeley]. We would also like to thank the following individuals for their efforts and the research that they conducted: Jasmine Cisneros, Cassandra Chung and Leanne Peres da Silva.
References
Alexander, C., Rietschel, E. T. (2001). Bacterial Lipopolysaccharides and Innate Immunity. J. Endotoxin Res. 7 (3), 167–202. doi: 10.1179/096805101101532675
Altan-Bonnet, G., Mukherjee, R. (2019). Cytokine-Mediated Communication: A Quantitative Appraisal of Immune Complexity. Nat. Rev. Immunol. 19 (4), 205–217. doi: 10.1038/s41577-019-0131-x
Alvarez, R. F., Barbuto, J. A. M., Venegeroles, R. (2019). A Nonlinear Mathematical Model of Cell-Mediated Immune Response for Tumor Phenotypic Heterogeneity. J. Theor. Biol. 471, 42–50. doi: 10.1016/j.jtbi.2019.03.025
Anderson, L. S., Reynolds, M. B., Rivara, K. R., Miller, L. S., Simon, S. I. (2019). A Mouse Model to Assess Innate Immune Response to Staphylococcus Aureus Infection. J. Vis. Exp. 144, 1–15. doi: 10.3791/59015
Badwey, J. A., Karnovsky, M. L. (1980). Active Oxygen Species and the Functions of Phagocytic Leukocytes. Annu. Rev. Biochem. 49 (1), 695–726. doi: 10.1146/annurev.bi.49.070180.003403
Brady, R., Do, F.-I., Tran, H. T., Janum, S., Møller, K., Brix, S., et al. (2016). Personalized Mathematical Model Predicting Endotoxin-Induced Inflammatory Responses in Young Men. Math. Model. Nat. Phenom. 13 (42), 1–20. doi: arXiv:1609.01570.
Brennan, K., Zheng, J. (2007). Interleukin 8. Elsevier Inc Xpharm: The Comprehensive Pharmacology Reference. xPharm 1, 1–4. doi: 10.1016/B978-008055232-3.61916-6
Caudill, L., Lynch, F. (2018). A Mathematical Model of the Inflammatory Response to Pathogen Challenge. Bull. Math. Biol. 80 (8), 2242–2271. doi: 10.1007/s11538-018-0459-6
Cavaillon, J. M. (2018). Exotoxins and Endotoxins: Inducers of Inflammatory Cytokines. Toxicon 149, 45–53. doi: 10.1016/j.toxicon.2017.10.016
Chandler, C. E., Ernst, R. K. (2017). Bacterial Lipids: Powerful Modifiers of the Innate Immune Response. F1000Res 6, 1–11. doi: 10.12688/f1000research.11388.1
Charchaflieh, J., Wei, J., Labaze, G., Hou, Y. J., Babarsh, B., Stutz, H., et al. (2012). The Role of Complement System in Septic Shock. Clin. Dev. Immunol. 2012:407324. doi: 10.1155/2012/407324
Chen, X., Wang, W., Xie, G., Hontecillas, R., Verma, M., Leber, A., et al. (2019). Multi-Resolution Sensitivity Analysis of Model of Immune Response to Helicobacter Pylori Infection via Spatio-Temporal Metamodeling. Front. Appl. Math. Stat. 5, 4 (4). doi: 10.3389/fams.2019.00004
Chow, C. C., Clermont, G., Kumar, R., Lagoa, C., Tawadrous, Z., Gallo, D., et al. (2005). The Acute Inflammatory Response in Diverse Shock States. Shock 24 (1), 74–84. doi: 10.1097/01.shk.0000168526.97716.f3
Cole, J., Aberdein, J., Jubrail, J., Dockrell, D. H. (2014). The Role of Macrophages in the Innate Immune Response to Streptococcus Pneumoniae and Staphylococcus Aureus: Mechanisms and Contrasts. Adv. Microb. Physiol. 65, 125–202. doi: 10.1016/bs.ampbs.2014.08.004
De Kimpe, S. J., Kengatharan, M., Thiemermann, C., Vane, J. R. (1995). The Cell Wall Components Peptidoglycan and Lipoteichoic Acid From Staphylococcus Aureus Act in Synergy to Cause Shock and Multiple Organ Failure. Proc. Natl. Acad. Sci. U.S.A. 92 (22), 10359–10363. doi: 10.1073/pnas.92.22.10359
de Waal Malefyt, R., Haanen, J., Spits, H., Roncarolo, M. G., te Velde, A., Figdor, C., et al. (1991). Interleukin 10 (IL-10) and Viral IL-10 Strongly Reduce Antigen-Specific Human T Cell Proliferation by Diminishing the Antigen-Presenting Capacity of Monocytes via Downregulation of Class II Major Histocompatibility Complex Expression. J. Exp. Med. 174 (4), 915–924. doi: 10.1084/jem.174.4.915
DiLeo, M. V., Kellum, J. A., Federspiel, W. J. (2009). A Simple Mathematical Model of Cytokine Capture Using a Hemoadsorption Device. Ann. BioMed. Eng. 37 (1), 222–229. doi: 10.1007/s10439-008-9587-8
Du, S. Q., Yuan, W. (2020). Mathematical Modeling of Interaction Between Innate and Adaptive Immune Responses in COVID-19 and Implications for Viral Pathogenesis. J. Med. Virol. 92 (9), 1615–1628. doi: 10.1002/jmv.25866
Eftimie, R., Gillard, J. J., Cantrell, D. A. (2016). Mathematical Models for Immunology: Current State of the Art and Future Research Directions. Bull. Math. Biol. 78 (10), 2091–2134. doi: 10.1007/s11538-016-0214-9
Emam, M., Tabatabaei, S., Sargolzaei, M., Sharif, S., Schenkel, F., Mallard, B. (2019). The Effect of Host Genetics on In Vitro Performance of Bovine Monocyte-Derived Macrophages. J. Dairy Sci. 102 (10), 9107–9116. doi: 10.3168/jds.2018-15960
Fan, H., Williams, D. L., Zingarelli, B., Breuel, K. F., Teti, G., Tempel, G. E., et al. (2007). Differential Regulation of Lipopolysaccharide and Gram-Positive Bacteria Induced Cytokine and Chemokine Production in Macrophages by Gαi Proteins. Immunology 122 (1), 116–123. doi: 10.1111/j.1365-2567.2007.02619.x
Faro, J., von Haeften, B., Gardner, R., Faro, E. (2019). A Sensitivity Analysis Comparison of Three Models for the Dynamics of Germinal Centers. Front. Immunol. 10:2038 (2038). doi: 10.3389/fimmu.2019.02038
Feezor, R. J., Oberholzer, C., Baker, H. V., Novick, D., Rubinstein, M., Moldawer, L. L., et al. (2003). Molecular Characterization of the Acute Inflammatory Response to Infections With Gram-Negative Versus Gram-Positive Bacteria. Infection Immun. 71 (10), 5803–5813. doi: 10.1128/IAI.71.10.5803-5813.2003
Feuerstein, G. Z., Liu, T., Barone, F. C. (1994). Cytokines, Inflammation, and Brain Injury: Role of Tumor Necrosis Factor-Alpha. Cerebrovasc Brain Metab. Rev. 6 (4), 341–360.
Foster, T. J. (2004). The Staphylococcus Aureus “Superbug”. J. Clin. Invest. 114 (12), 1693–1696. doi: 10.1172/JCI23825
Fournier, B., Philpott, D. J. (2005). Recognition of Staphylococcus Aureus by the Innate Immune System. Clin. Microbiol. Rev. 18 (3), 521–540. doi: 10.1128/cmr.18.3.521-540.2005
Hergott, C. B., Roche, A. M., Tamashiro, E., Clarke, T. B., Bailey, A. G., Laughlin, A., et al. (2016). Peptidoglycan From the Gut Microbiota Governs the Lifespan of Circulating Phagocytes at Homeostasis. Blood 127 (20), 2460–2471. doi: 10.1182/blood-2015-10-675173
Hoebe, K., Janssen, E., Beutler, B. (2004). The Interface Between Innate and Adaptive Immunity. Nat. Immunol. 5 (10), 971–974. doi: 10.1038/ni1004-971
Holmes, W. E., Lee, J., Kuang, W. J., Rice, G. C., Wood, W. I. (1991). Structure and Functional Expression of a Human Interleukin-8 Receptor. Science 253 (5025), 1278–1280. doi: 10.1126/science.1840701
Ip, W. K. E., Hoshi, N., Shouval, D. S., Snapper, S., Medzhitov, R. (2017). Anti-Inflammatory Effect of IL-10 Mediated by Metabolic Reprogramming of Macrophages. Science 356 (6337), 513–519. doi: 10.1126/science.aal3535
Iyer, S. S., Cheng, G. (2012). Role of Interleukin 10 Transcriptional Regulation in Inflammation and Autoimmune Disease. Crit. Rev. Immunol. 32 (1), 23–63. doi: 10.1615/critrevimmunol.v32.i1.30
Johnson, D. B., Lopez, M. J., Kelley, B. (2021). Dexamethasone (StatPearls. (Treasure Island (FL: StatPearls Publishing Copyright © 2021, StatPearls Publishing LLC).
Khoa, N. D., Montesinos, M. C., Williams, A. J., Kelly, M., Cronstein, B. N. (2003). Th1 Cytokines Regulate Adenosine Receptors and Their Downstream Signaling Elements in Human Microvascular Endothelial Cells. J. Immunol. 171 (8), 3991–3998. doi: 10.4049/jimmunol.171.8.3991
Kim, J., Bentley, P., Aiklelin, U., Greesmith, J., Tedesco, G., Twycross, J. (2007). Immune System Approaches to Intrusion Detection—A Review. Natural Computing. Natural Comput. 6, 413–466. doi: 10.1007/s11047-006-9026-4
Kwiecinski, J. M., Horswill, A. R. (2020). Staphylococcus Aureus Bloodstream Infections: Pathogenesis and Regulatory Mechanisms. Curr. Opin. Microbiol. 53, 51–60. doi: 10.1016/j.mib.2020.02.005
Laarman, A., Milder, F., van Strijp, J., Rooijakkers, S. (2010). Complement Inhibition by Gram-Positive Pathogens: Molecular Mechanisms and Therapeutic Implications. J. Mol. Med. (Berl) 88 (2), 115–120. doi: 10.1007/s00109-009-0572-y
Lan, F., Zhang, N., Holtappels, G., De Ruyck, N., Krysko, O., Van Crombruggen, K., et al. (2018). Staphylococcus Aureus Induces a Mucosal Type 2 Immune Response via Epithelial Cell-Derived Cytokines. Am. J. Respir. Crit. Care Med. 198 (4), 452–463. doi: 10.1164/rccm.201710-2112OC
Liu, X., Yin, S., Chen, Y., Wu, Y., Zheng, W., Dong, H., et al. (2018). LPS−induced Proinflammatory Cytokine Expression in Human Airway Epithelial Cells and Macrophages via NF−κb, STAT3 or AP−1 Activation. Mol. Med. Rep. 17 (4), 5484–5491. doi: 10.3892/mmr.2018.8542
Lowings, K. M., Christodoulides, M., Warner, J. A. (2009). Cytokine Response to Staphylococcus Aureus (SA) in Human Lung, Skin and Nasal Tissue. J. Allergy Clin. Immunol. 123 (2), S51. doi: 10.1016/j.jaci.2008.12.161
Marchuk, G. I. (1997). Mathematical Modelling of Immune Response in Infectious Diseases (Netherlands: Springer).
McNicholas, S., Talento, A. F., O’Gorman, J., Hannan, M. M., Lynch, M., Greene, C. M., et al. (2014). Cytokine Responses to Staphylococcus Aureus Bloodstream Infection Differ Between Patient Cohorts That Have Different Clinical Courses of Infection. BMC Infect. Dis. 14, 580. doi: 10.1186/s12879-014-0580-6
Meier-Schellersheim, M., Varma, R., Angermann, B. R. (2019). Mechanistic Models of Cellular Signaling, Cytokine Crosstalk, and Cell-Cell Communication in Immunology. Front. Immunol. 10, 2268 (2268). doi: 10.3389/fimmu.2019.02268
Michalopoulos, A. S., Livaditis, I. G., Gougoutas, V. (2011). The Revival of Fosfomycin. Int. J. Infect. Dis. 15 (11), e732–e739. doi: 10.1016/j.ijid.2011.07.007
Mirzaei, R., Ranjbar, R., Karampoor, S., Goodarzi, R., Hasanvand, H. (2020). The Human Immune System Toward Staphylococcus Aureus. Open Microbiol. J. 14, 164–170. doi: 10.2174/1874285802014010164
Möbius, W., Laan, L. (2015). Physical and Mathematical Modeling in Experimental Papers. Cell 163 (7), 1577–1583. doi: 10.1016/j.cell.2015.12.006
Monguió-Tortajada, M., Franquesa, M., Sarrias, M. R., Borràs, F. E. (2018). Low Doses of LPS Exacerbate the Inflammatory Response and Trigger Death on TLR3-Primed Human Monocytes. Cell Death Dis. 9 (5), 499. doi: 10.1038/s41419-018-0520-2
Moreillon, P., Majcherczyk, P. A. (2003). Proinflammatory Activity of Cell-Wall Constituents From Gram-Positive Bacteria. Scand. J. Infect. Dis. 35 (9), 632–641. doi: 10.1080/00365540310016259
Morel, P. A., Lee, R. E. C., Faeder, J. R. (2017). Demystifying the Cytokine Network: Mathematical Models Point the Way. Cytokine 98, 115–123. doi: 10.1016/j.cyto.2016.11.013
Morikawa, K., Watabe, H., Araake, M., Morikawa, S. (1996). Modulatory Effect of Antibiotics on Cytokine Production by Human Monocytes. vitro. Antimicrob. Agents Chemother. 40 (6), 1366–1370. doi: 10.1128/aac.40.6.1366
Nijhout, H. F., Best, J. A., Reed, M. C. (2015). Using Mathematical Models to Understand Metabolism, Genes, and Disease. BMC Biol. 13 (1), 79. doi: 10.1186/s12915-015-0189-2
Oh, S.-J., Lee, J. K., Shin, O. S. (2019). Aging and the Immune System: The Impact of Immunosenescence on Viral Infection, Immunity and Vaccine Immunogenicity. Immune Network 19 (6), e37–e37. doi: 10.4110/in.2019.19.e37
Osawa, Y., Nagaki, M., Banno, Y., Brenner, D. A., Asano, T., Nozawa, Y., et al. (2002). Tumor Necrosis Factor Alpha-Induced Interleukin-8 Production via NF-kappaB and Phosphatidylinositol 3-Kinase/Akt Pathways Inhibits Cell Apoptosis in Human Hepatocytes. Infection Immun. 70 (11), 6294–6301. doi: 10.1128/IAI.70.11.6294-6301.2002
Pigozzo, A. B., Macedo, G. C., dos Santos, R. W., Lobosco, M. (2013). On the Computational Modeling of the Innate Immune System. BMC Bioinf. 14 (6), S7. doi: 10.1186/1471-2105-14-S6-S7
Ponnappan, S., Ponnappan, U. (2011). Aging and Immune Function: Molecular Mechanisms to Interventions. Antioxidants Redox Signaling 14 (8), 1551–1585. doi: 10.1089/ars.2010.3228
Quintela, B. D. M., dos Santos, R. W., Lobosco, M. (2014). On the Coupling of Two Models of the Human Immune Response to an Antigen. BioMed. Res. Int. 2014, 410457. doi: 10.1155/2014/410457
Sarah, D. H., Richard, B. B. (2009). Pneumonia Caused by Staphylococcus Aureus. Curr. Respir. Med. Rev. 5 (1), 62–67. doi: 10.2174/157339809787354056
Smith, E. J., Visai, L., Kerrigan, S. W., Speziale, P., Foster, T. J. (2011). The Sbi Protein is a Multifunctional Immune Evasion Factor of Staphylococcus Aureus. Infection Immun. 79 (9), 3801–3809. doi: 10.1128/IAI.05075-11
Spaan, A. N., Surewaard, B. G., Nijland, R., van Strijp, J. A. (2013). Neutrophils Versus Staphylococcus Aureus: A Biological Tug of War. Annu. Rev. Microbiol. 67, 629–650. doi: 10.1146/annurev-micro-092412-155746
Swain, S. L., Nikolich-Zugich, J. (2009). Key Research Opportunities in Immune System Aging. Journals Gerontology. Ser. A Biol. Sci. Med. Sci. 64 (2), 183–186. doi: 10.1093/gerona/gln068
Tawfik, D. M., Lankelma, J. M., Vachot, L., Cerrato, E., Pachot, A., Wiersinga, W. J., et al. (2020). Comparison of Host Immune Responses to LPS in Human Using an Immune Profiling Panel, In Vivo Endotoxemia Versus Ex Vivo Stimulation. Sci. Rep. 10 (1), 9918. doi: 10.1038/s41598-020-66695-2
Torres, M., Wang, J., Yannie, P. J., Ghosh, S., Segal, R. A., Reynolds, A. M. (2019). Identifying Important Parameters in the Inflammatory Process With a Mathematical Model of Immune Cell Influx and Macrophage Polarization. PloS Comput. Biol. 15 (7), e1007172. doi: 10.1371/journal.pcbi.1007172
Trinchieri, G. (2003). Interleukin-12 and the Regulation of Innate Resistance and Adaptive Immunity. Nat. Rev. Immunol. 3 (2), 133–146. doi: 10.1038/nri1001
Wang, Z., Deisboeck, T. S. (2014). Mathematical Modeling in Cancer Drug Discovery. Drug Discovery Today 19 (2), 145–150. doi: 10.1016/j.drudis.2013.06.015
Wang, J. E., Jørgensen, P. F., Almlöf, M., Thiemermann, C., Foster, S. J., Aasen, A. O., et al. (2000). Peptidoglycan and Lipoteichoic Acid From Staphylococcus Aureus Induce Tumor Necrosis Factor Alpha, Interleukin 6 (IL-6), and IL-10 Production in Both T Cells and Monocytes in a Human Whole Blood Model. Infect. Immun. 68 (7), 3965–3970. doi: 10.1128/iai.68.7.3965-3970.2000
Winer, D. A., Winer, S., Shen, L., Wadia, P. P., Yantha, J., Paltser, G., et al. (2011). B Cells Promote Insulin Resistance Through Modulation of T Cells and Production of Pathogenic IgG Antibodies. Nat. Med. 17 (5), 610–617. doi: 10.1038/nm.2353
Winslow, R. L., Trayanova, N., Geman, D., Miller, M. I. (2012). Computational Medicine: Translating Models to Clinical Care. Sci. Transl. Med. 4 (158), 158rv111. doi: 10.1126/scitranslmed.3003528
Xue, X., Falcon, D. M. (2019). The Role of Immune Cells and Cytokines in Intestinal Wound Healing. Int. J. Mol. Sci. 20 (23), 6097. doi: 10.3390/ijms20236097
Yao, L., Lowy, F. D., Berman, J. W. (1996). Interleukin-8 Gene Expression in Staphylococcus Aureus-Infected Endothelial Cells. Infection Immun. 64 (8), 3407–3409. doi: 10.1128/iai.64.8.3407-3409.1996
Zhang, J. M., An, J. (2007). Cytokines, Inflammation, and Pain. Int. Anesthesiol Clin. 45 (2), 27–37. doi: 10.1097/AIA.0b013e318034194e
Zheng, Y., Q, C., Zhang, X., Zhu, Y., Li, A., Wang, M., et al. (2020). The Tst Gene Associated Staphylococcus Aureus Pathogenicity Island Facilitates its Pathogenesis by Promoting the Secretion of Inflammatory Cytokines and Inducing Immune Suppression. Microbial Pathogenesis 138, 1–27. doi: 10.1016/j.micpath.2019.103797
Keywords: cytokines, mathematical modeling, immune response, immune system, Staphycoccus aureus, cytokine response, cell activation
Citation: Talaei K, Garan SA, Quintela BdM, Olufsen MS, Cho J, Jahansooz JR, Bhullar PK, Suen EK, Piszker WJ, Martins NRB, Moreira de Paula MA, dos Santos RW and Lobosco M (2021) A Mathematical Model of the Dynamics of Cytokine Expression and Human Immune Cell Activation in Response to the Pathogen Staphylococcus aureus. Front. Cell. Infect. Microbiol. 11:711153. doi: 10.3389/fcimb.2021.711153
Received: 18 May 2021; Accepted: 14 October 2021;
Published: 10 November 2021.
Edited by:
Luciana Balboa, Academia Nacional de Medicina, ArgentinaReviewed by:
Maria Cernada, La Fe University Hospital, SpainSurya Prakash Pandey, University of Pittsburgh, United States
Copyright © 2021 Talaei, Garan, Quintela, Olufsen, Cho, Jahansooz, Bhullar, Suen, Piszker, Martins, Moreira de Paula, dos Santos and Lobosco. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Kian Talaei, S2lhbl90OTdAYmVya2VsZXkuZWR1