The model was developed following the four-step sequence proposed by the system dynamics methodology (25). First, experimental data and other evidence were used to create a mental modeling of the response of monocytes to endotoxins. Second, the model structure able to explain the polarization of the monocytes from the proinflammatory to the ET phenotype was represented as a Forrester diagram (Figure 1A). Third, the system was mathematically modeled as a continuous dynamic process represented by a set of differential and algebraic equations (Tables 1–3). Finally, the model was parameterized (Table 4), optimized, and validated to fit the experimental data from Figures 2–6 and findings published in previous articles (6, 33, 34). The parameters of the Hill function corresponding to the activation rate were adjusted by using a non-linear regression analysis (Figure 3) to fit the experimental data summarized in Table 5. Similarly, the parameters of the Hill function corresponding to the TNF synthesis rate were adjusted by using a non-linear regression analysis (Figure 4C) to fit the estimated TNF synthesis rates for each LPSs concentration and summarized in Table 6. The rest of the model parameters were adjusted using the Vensim’s optimizer and Matlab to get the best match between model behavior and the data (Figures 4B and 5C). These adjustments were performed by minimizing the total quadratic error, which measures the goodness of the model. Finally, the model was validated using datasets from Figures 5B and 6.

The modeling process, simulations, and optimization analyses were performed using Vensim DSS software, version 5.7a (Ventana Systems, Harvard, MA, USA) and Matlab R2015a (Mathworks, Inc., MA, USA).

Peripheral blood was obtained from three patients with sepsis (mean age ± SD: 68 ± 10.6 years) who had microbiologically confirmed Gram-negative bacteremia (positive blood cultures for Escherichia coli), secondary to a urinary tract infection. The patients who met the consensus conference definition of sepsis (35) were admitted to the Department of Internal Medicine Service at La Paz University Hospital. Blood samples were collected from the patients within 24 h of sepsis confirmation by blood culture. The following exclusion criteria were imposed: the presence of malignancy or chronic inflammatory diseases, treatment with steroids or immunosuppressive drugs during the last month, hepatic failure (serum aspartate aminotransferase and/or alanine aminotransferase > 100 IU/l; prothrombin time < 60%; total bilirubin < 60 μmol/l), renal insufficiency (plasma creatinine > 200 μmol/l), HIV/AIDS, hepatitis B or C, pregnancy, and age > 80 years. Peripheral blood from 18 healthy volunteers was also obtained. Of them, 12 were used in setting up the model in Figures 2, 4 and 5, and the other 6 were used in validating the model in Figure 6. All the procedures were in accordance with the Helsinki Declaration of 2000, and informed consent was obtained from all the participants. This study was approved by the La Paz University Hospital Ethics Committee.

Peripheral blood mononuclear cells were isolated from healthy volunteers or sepsis patients by centrifugation on Ficoll-Hypaque Plus (Amersham Biosciences), and monocytes were obtained by adherence as previously described (6). The purity of the monocyte cultures was tested by CD14 labeling and flow cytometry analysis (average 82% CD14 + cells). Other cell surface markers were also tested (average CD1a = 4.1%, and CD89 = 85%). All the reagents used for the cell cultures were endotoxin-free, as assayed with the limulus amebocyte lysate test (Cambrex).

To establish the model parameters, the monocytes were treated with LPS concentrations ranging from 5 ng/ml to 0.5 pg/ml during the times of stimulation (t_s) indicated in Figure 2. After LPS treatment, the cells were washed three times with PBS and kept in complete medium for 16 h to ensure the monocytes entered into ET. The cells were then restimulated with 5 ng/ml LPS for 3 h and cytokine levels were quantified. The control cells were not treated with LPS during the first stimulation.

To compare the model output with experimental data, monocytes from healthy volunteers and patients with sepsis were treated once or twice with LPSs at the concentrations indicated in Figures 4–6. Samples of culture supernatant were collected at times ranging from 1 to 24 h to quantify the levels of TNFα and CCL2.

The LPS from Salmonella abortus was a kind gift from Dr. Galanos (Max Planck Institute for Immunobiology, Freiburg, Germany). All work with LPS was performed in biosafety level (BSL) two facilities under appropriate work practices, and proper use of personal protective equipment. The medium used for the cell culture was DMEM from Invitrogen.

The cytokine levels (TNFα and CCL2) in the culture supernatants from the human samples were determined using the CBA Flex Set (BD Biosciences) following the manufacturer’s protocol. The data collected were analyzed by flow cytometry, using a BD FACSCalibur flow cytometer (BD Biosciences).

The data were collected and are expressed as mean ± SEM (n is indicated in each figure legend). The statistical significance was set at p < 0.05 and all statistical analyses were conducted using Prism 5.0 software (GraphPad) and Matlab.

A dynamic compartmental model can be used to properly parameterize and describe several components of the immune system (36, 37). In this sense, both physical sites and functions, such as the organs of the immune system and states of immune cell differentiation, respectively, can be schematized as a set of compartments.

To explain the transition of the peripheral blood monocytes through the various functional stages after an endotoxin challenge, we considered a dynamic model divided into two subsystems (Figures 1A,B). The main subsystem reproduces the activation/immunomodulation process by which monocytes transit after endotoxin challenge. It contains four compartments (states) and three transitions (physical flows) between them (Figure 1A), which are based on various assumptions. First, the monocytes are classified into three subgroups: resting monocytes, proinflammatory monocytes, and ET monocytes. Second, the transition from the resting monocytes to the proinflammatory monocytes (activation) only occurs when LPSs are present in the medium. Third, the flow between proinflammatory monocytes and ET monocytes (immunomodulation) is an LPS-independent process. Finally, we assumed that LPS concentration is not continually the same; thus, an outflow was included in the model (LPS removal).

The model also contains a second subsystem modeling the production of TNF (as a proinflammatory marker) and CCL2 (as an ET marker) (Figure 1B). Because these cytokines are differentially synthesized by the three subsets of monocytes, this subsystem is used to indirectly monitor the dynamic process followed by monocytes modeled in the main subsystem. It is composed of two states (TNF and CCL2) and four flows (TNF production, TNF degradation, CCL2 production, and CCL2 degradation), modifying the balance of these stocks.

The model structure, built with Vensim software (Figures 1A,B), contains the states and flows previously mentioned, two auxiliary variables (inside circle), and 14 parameters (in bold). These elements are linked by the physical flows (double line with arrow) and by the information transmissions (single line with arrow), according to the mathematical model represented by the set of differential and algebraic equations (Tables 1–3).

The six differential equations of Table 1 establish the mass balance (inflows minus outflows) in the compartments. The first three differential equations of Table 1 describe the changes in the number of monocytes in the three subgroups. Because LPS concentration in the medium is known to decrease over time, changes in LPSs are also represented by a differential equation (Table 1, Eq. 4). Equations 5 and 6 of Table 1 express the changes in TNF and CCL2 levels per time unit, respectively.

The activation of monocytes after LPS stimulation depends directly on the stock of resting monocytes and the auxiliary variable activation rate (Table 2, Eq. 1). The activation rate is modulated by a Hill activation function (Table 3, Eq. 1). This expression depends on: the maximum activation rate, which allows setting the top value that the auxiliary variable can take; the AC50, which sets the concentration of LPS needed to reach 50% of the maximum activation rate; the n Hill coefficient, which describes the steepness of the curve; and the LPS concentration expressed as a stock variable. The immunomodulation depends directly on the amount of proinflammatory monocytes and the parameter immunomodulation rate (Table 2, Eq. 2). Similarly, LPS removal at a given time is directly proportional to the LPS concentration at that time and the LPS removal rate parameter (Table 2, Eq. 3).

Equation 4 of Table 2 expresses the TNF production that depends on the amount of proinflammatory monocytes, the auxiliary variable TNF synthesis rate, and the TD time delay in TNF parameter. Because the signaling cascade initiated after the binding of LPSs to the receptors in the plasmatic membrane leads to a delay between the time in which the resting monocytes are being activated and the production of TNF, the parameter TD time delay in TNF was incorporated to model this fact. In addition, TNFα expression in response to LPSs is regulated by the translocation of NF-κB to the nucleus. NF-κB translocation is highly dependent on the LPS concentration in the medium and it has been shown that the LPS-dose–response curve of NF-κB translocation follows a Hill equation (34). Therefore, in the model the TNF synthesis rate is modulated by a Hill activation function (Table 3, Eq. 2). This expression depends on: the maximum TNF synthesis rate, which is the top value that the auxiliary variable can take; the SC50 which sets the concentration of LPS needed to reach half of the maximum TNF synthesis rate; the m Hill coefficient which defines the slope of the curve; and the LPS concentration.

Because CCL2 is produced by the three subtypes of monocytes at different rates, CCL2 production is represented by an equation that considers the three states, each of which is multiplied by its corresponding parameter (Table 2, Eq. 6). The TNF degradation and CCL2 degradation outflows are proportional to the mass in the source stock (TNF and CCL2, respectively) and the TNF degradation rate or CCL2 degradation rate parameters, respectively (Table 2, Eqs 5 and 7, respectively).

A model’s performance is highly dependent on its parameters. Thus, defining the values of the model’s parameters is one of the most important processes in system dynamics modeling. Although this model is relatively simple, the parameterization was challenging due to the uncertainty of the parameter values. As a consequence, there is a wide range of value combinations, many of which can result in unexpected model outputs far from the reality. To solve this problem, the parameters were sequentially estimated by isolating the different processes (e.g., activation, TNF production, CCL2 production, etc.) from the rest. Thus, the values were estimated by using the experimental data from Figures 2–5. They were optimized following various approaches, such as non-linear regression analysis, the Vensim’s optimization tool, or using discrete models in Matlab to set the final values for the model parameters that are summarized in Table 4. These approaches require to select initial values and ranging values for some parameters that are not being optimized in that moment. To this end, findings published in previous articles were used (6, 33, 34).

As described above, the Hill function to express the activation rate depends on the maximum activation rate, the AC50, and n parameters. To obtain the value for these three parameters, the activation rate was estimated for monocytes from healthy volunteers using various LPSs concentrations and following the experimental design described in Figure 2A. Since TNF is produced by the proinflammatory subset of monocytes, this cytokine was used to indirectly estimate the percentage of cells that transit to the proinflammatory state after an LPSs challenge. Thus, the relative TNFα expression levels after a second LPSs challenge with 5 ng/ml LPS were used to estimate the percentage of resting monocytes that were not activated during the first stimulus (Figures 2B–F). As expected, when the time of the first challenge with LPS increased, there was a decrease in the fraction of cells remaining in the resting state after the first challenge (Figures 2B–F). In addition, this decrease was faster when the LPS concentrations were higher (compare Figures 2B,F). The activation rates for each LPSs concentration were estimated from the adjusted one-step decay functions obtained by non-linear regression analysis of the experimental data (Figures 2B–F) and are summarized in Table 5. Since the activation rates at super-low LPS doses are longer, the experimental design was slightly changed. Thus, the monocytes were first challenged with increasing concentrations of LPS, ranging from 0.5 pg/ml to 5 ng/ml, following the scheme of Figure 2G. The first stimulus was maintained for 16 h to significantly reduce the relative expression of TNFα in the second challenge with 5 ng/ml LPS. The percentage of resting monocytes that were not activated during the first stimulus was estimated as described above. The activation rates for LPSs concentrations 0.5, 5, and 50 pg/ml were calculated from a one-step decay function adjusted to pass through the experimental data that are also summarized in Table 5. Overall, these results showed that the LPS concentration highly impacted the activation rate (Table 5). Interestingly, no priming effect was observed in our system even with LPSs concentrations as low as 0.5 pg/ml (Figure 2H) or 0.05 pg/ml (data not shown). Graphic representation of the estimated activation rates from Table 5 in relation to the LPS concentration (Figure 3) showed a sigmoidal tendency that match with an activation Hill equation. Thus, the maximum activation rate, the AC50, and n parameters were obtained from the activation rate equation adjusted to fit the activation rates for each LPSs concentration (Figure 3) estimated from the experimental data shown in Table 5, which were set to 4.499/h, 0.1889 ng/ml, and 1.4825, respectively. Thus, the half-maximal activation rate takes place when the LPS concentration equals AC50, and the almost maximal rate occurs when LPS quadruples the value of this parameter.

Similarly to the activation rate, the TNF synthesis rate expressed as a Hill function depends on three parameters. To optimize the maximum TNF synthesis rate, the SC50 and the m parameters, the TNF synthesis rate was estimated using the Vensim’s optimization tool to fit the experimental data obtained from scenarios using various amounts of LPS (Figures 4A,B). As expected, the synthesis of TNF was highly dependent to the LPSs doses (Figure 4B) which led to different estimated TNF synthesis rates depending on the LPS concentration (Table 6). Representation of the estimated TNF synthesis rates (Table 6) showed a sigmoidal trend that indicated that the process can be explained by an activation Hill function (Figure 4C). The estimated values were used to adjust the Hill equation (Figure 4C) from which the maximum TNF synthesis rate, SC50, and m parameters were determined which were 0.0071 pg/h/monocytes, 0.1407 ng/ml, and 1.767, respectively. However, when these values were used to simulate the various scenarios, significant deviations between the model output and the experimental data were observed for lower LPSs doses (data not shown). In this sense, the values for maximum TNF synthesis rate and m were considered as definitive. However, the SC50 was further adjusted together with the immunomodulation rate, LPS removal rate, and TNF degradation rate parameters (Table 4) using a discrete model built in Matlab to obtain the best fitting considering all scenarios from Figure 4B. Using the adjusted values (Table 4), the model highly reproduces the shape of the TNF curves obtained after simulation of isolated monocytes from healthy volunteers with a wide range of LPSs doses (Figure 4B).

Similarly, to optimize the model parameters related to CCL2 synthesis and degradation, a Matlab discrete model was used to fit the data from an experimental ET model obtained, following the design described in Figure 5A. In agreement with what has been previously published (19), ET monocytes showed a reduced capacity to produce TNF whereas CCL2 was upregulated (Figures 5B,C). As with the initial value of LPSs in the medium, factors related to the initial characteristics of the monocyte subsets have been set as user-defined. Therefore, setting different values for these variables allows the simulation of particular scenarios, such as the response of ET monocytes to LPSs. To simulate the ET scenario, the model was programmed with 0, 75,000, and 925,000 of the initial monocyte population in the resting, proinflammatory, and ET states, respectively, and then simulated with initial concentration of 5 ng/ml of LPS. The TNF curves obtained after simulation of the control and ET scenarios were similar to those observed in the cultures of monocytes (Figure 5B) which prove the fitness of the model and validate the optimization process of the TNF synthesis parameters explained above. The parameters CCL2 synthesis rate on resting, CCL2 synthesis rate on inflammation, CCL2 synthesis rate on ET, and CCL2 degradation rate were automatically optimized to fit the CLL2 levels measured in monocytes cultures (Figure 5C). As expected, both the TNF and CCL2 curves obtained after simulation of the ET scenario were markedly different from those obtained in a scenario programmed with all the monocytes in the resting state (Figures 5B,C, compare black and blue lines). Thus, the model showed a lower production of TNF and higher production of CLL2 when simulating the ET scenario compared with the control conditions. These results indicate that this model is able to reproduce the ET phenomenon.

Since initial monocyte population in the resting and ET levels can be considered dependent on each other, a tolerance index parameter was included in a modified version of the model. This parameter specifies the initial percentage of the initial monocytes that are in ET and is used to calculate the initial population in ET and resting states (Table 7). Since the transition of the monocytes through the proinflammatory state is a rapid process, the initial proinflammatory monocytes value can be considered equal to 0.

To further analyze whether our model is able to simulate intermediate ET scenarios, monocytes from healthy volunteers were challenged with 0.5 ng/ml LPSs for 30 min and then incubated in fresh medium for 16 h allowing the culture to reach the ET (Figure 6A). Similar to results in Figure 5, this led to a decrease in the production of TNF and an increase in the production of CCL2 in comparison with the control scenario (Figures 6B,C, compare gray triangles with black circles). By setting the Tolerance index value in the model to 68%, the TNF and CCL2 outputs reproduce the cytokine curves obtained after a second LPS challenge with 5 ng/ml LPSs (Figures 6B,C), indicating that the model is able to reproduce different ET intensities in monocyte cultures.

We and other authors have previously reported that monocytes derived from patients with sepsis are in an ET state in which they show a diminished inflammatory response after ex vivo LPSs challenge (7, 8). Taking this into account, the model was programmed with a Tolerance index of 68%, and then simulated with initial concentration of 5 ng/ml of LPSs. The TNF and CCL2 levels at 16 h of simulation were similar to those observed in the cultures of monocytes isolated from patients with sepsis and challenged with 5 ng/ml LPSs for 16 h (Figures 6D,E) reinforcing the idea that this model is able to reproduce the ET exhibited by cultured monocytes from patients with sepsis when they are challenged with LPSs, and that it could be used to determine the tolerance status of the monocytes of these patients.