A multi-approach and multi-scale platform to model CD4+ T cells responding to infections

Immune responses rely on a complex adaptive system in which the body and infections interact at multiple scales and in different compartments. We developed a modular model of CD4+ T cells, which uses four modeling approaches to integrate processes at three spatial scales in different tissues. In each cell, signal transduction and gene regulation are described by a logical model, metabolism by constraint-based models. Cell population dynamics are described by an agent-based model and systemic cytokine concentrations by ordinary differential equations. A Monte Carlo simulation algorithm allows information to flow efficiently between the four modules by separating the time scales. Such modularity improves computational performance and versatility and facilitates data integration. We validated our technology by reproducing known experimental results, including differentiation patterns of CD4+ T cells triggered by different combinations of cytokines, metabolic regulation by IL2 in these cells, and their response to influenza infection. In doing so, we added multi-scale insights to single-scale studies and demonstrated its predictive power by discovering switch-like and oscillatory behaviors of CD4+ T cells that arise from nonlinear dynamics interwoven across three scales. We identified the inflamed lymph node’s ability to retain naive CD4+ T cells as a key mechanism in generating these emergent behaviors. We envision our model and the generic framework encompassing it to serve as a tool for understanding cellular and molecular immunological problems through the lens of systems immunology.


Introduction
Immune responses mediated by CD4+ T cells involve complex interactions among immune cells and molecules. Resting CD4+ T cells are activated by antigen-presenting cells and cytokines, further differentiate, and secrete cytokines to act against pathogens and abnormal cells. They also recruit other immune cells to the sites of infection. Depending on the cytokine milieu, activated CD4+ T cells may differentiate into various phenotypes, including T helper type 1 (Th1), T helper type 2 (Th2), T helper type 17 (Th17), and induced T regulatory cells (Tregs) [1][2][3]. Activated CD4+ T cells utilize certain signaling and metabolic pathways, such as aerobic glycolysis, to produce the energy and molecular precursors required to achieve a specific mixture of phenotypes [1,4]. To fully understand these complex interactions underlying the dynamics of CD4+ T cell immune response, we must integrate events taking place at various spatial, temporal, and organizational scales, such as immune cell proliferation, development, and migration; cell-cell and cell-molecule interactions; intracellular signaling; and intracellular metabolism. In other words, it is a systems immunology problem.
Multi-scale modeling aims to integrate spatial, temporal, and organizational scales of biological systems. Such integration could be achieved by combining different modeling approaches, such as ordinary differential equations (ODEs) and partial differential equations (PDEs) for the chemical kinetics and transport of molecular species (in terms of concentrations) in and across different cells, organs, or tissues; agent-based modeling (ABM) for cells and molecules interacting in heterogeneous spatial environments; and ODE/PDE-based, rulebased, logic-based, or constraint-based models for intracellular dynamics. Genome-Scale Metabolic Models (GSMMs) have been developed for various human tissues and cell types [5,6] and can be used to model whole-cell metabolism in a context-specific manner. Logic-based models can be used to model cell signaling and gene regulation where kinetic parameters are unavailable [7]. This integrative strategy has been used extensively in immunology [8]. For instance, multi-scale models have been developed to study infections and inflammatory processes [9] and the immune response to the Helicobacter pylori infection [10]. In a study of the role of TNF-α in the formation of granulomas during Mycobacterium tuberculosis infections [11], a macrophage and T cell ABM was combined with a PDE model of the extracellular dynamics of soluble molecules and an ODE model of the intracellular signaling events. A coarse-grained version was used to model multiple granulomas as agents, each comprising a system of ODEs to represent the concentration dynamics of cells, cytokines, and bacteria in one granuloma [12]. Another extension was built by placing the ABM in a multi-compartment setup to shed light on the trafficking of antigen-presenting cells from the lung to the lymph node as a control mechanism for immunity against tuberculosis [13]. Beyond the study of tuberculosis, compartmental modeling was used to integrate ODE models of different parts of the immune system [14]. Multi-scale models have also been used to study the heterogeneous tumor immune microenvironment [15]. Our study was motivated by the desire to build on these studies by adding more mechanistic details about CD4+ T cells in terms of signal transduction (via logical models) and genome-wide metabolism (via constraint-based models).
Herein, with a view to comprehending the immune system as an integrated whole, we present a new multi-approach and multi-scale modeling framework that can be used to model diverse immune responses at molecular, cellular, and systemic scales. We demonstrated its capabilities by modeling the dynamics of CD4+ T cells in response to influenza infections (real and hypothetical) in different cytokine milieus, considering heterogeneous populations of Th0, Th1, Th2, Th17, and Treg subtypes and 11 cytokines in three spatial compartments (an infection site or target organ, lymphoid tissues, and a circulatory system). For each subtype, we modeled the intracellular signaling and gene regulation network underpinning the cell's proliferation, differentiation into that subtype, as well as the whole-cell metabolism of that subtype. Using the framework, we reproduced experimental results available in literature, identified and explained new behaviors of CD4+ T cells, proving it is a tool wherein new regulations within the immune system may be implemented, and whose predictions may be experimentally validated.

Four different modeling approaches represent processes spanning three spatial scales in three compartments
The multi-scale model is divided into three compartments (Fig 1). The first target organ compartment represents the site of initial insult (e.g., infection). Within the context of the validation and robustness studies presented below, pertaining to influenza infections, this compartment represents the lungs. The second compartment is an abstract representation of the lymphoid tissues connected to the target organ: a micro-environment where CD4+ T cells can be stimulated by antigens. In the case of influenza, this compartment represents the draining lymph nodes nearest the lungs. The final compartment is the circulatory system, which is an abstraction of the rest of the body. Each compartment is well-stirred with respect to both cells and cytokines. The CD4+ T cells in each compartment can sense and/or secrete 11 cytokines, and they are stimulated by a study-dependent and time-varying signal which represents the antigen load and the activity of the rest of the immune system. A CD4+ T cell is created as a naive cell with the Th0 phenotype and typically goes through three major stages: activation, expansion, and contraction (including memory formation and cell death). The behaviors and attributes of each cell are dynamic and depend on the outputs of the intracellular models of signal transduction, gene regulation, and metabolism. The processes related to signal transduction and gene regulation include the canonical pathways regulating the transdifferentiation of effector CD4+ T cells (interconversion among mixtures of Th0, Th1, Th2, Th17, and Treg phenotypes) [2], and events such as aerobic glycolysis, memory formation, and cell death. These pathways are detailed in S1 Text. Although the metabolic models consider whole-cell metabolism at the genome-scale, different subsets of metabolic fluxes are used to model different parts of the cell cycle in the multi-scale model. For instance, the fluxes related to glutaminolysis are active during the S phase when DNA is synthesized. Overall structure of the model. Schema representing the arrangement of models at three spatial scales. Cells are represented by discrete, autonomous agents that can differentiate into different subtypes, divide, and die. Each cell contains a constraint-based metabolic model for each phenotype and a logical model of signaling and gene regulatory pathways. The The aforementioned processes are represented using four different modeling frameworks. At the systemic scale, 11 sets of three coupled compartmental ordinary differential equations model the temporal dynamics of the 11 cytokine concentrations in the three compartments. Each cytokine is secreted and/or sensed by the T cells and is carried between compartments by blood and lymph (Fig 2A). At the cellular scale, CD4+ T cells are described by an agent-based model. In the activation stage, a naive agent migrates from the target organ to the lymphoid tissues (Fig 2A), where it is stimulated until activation into an effector agent [16][17][18]. In the expansion stage, it changes its phenotype adaptively in response to changes in cytokine concentrations, undergoes cell cycling, produces cytokines, and migrates to the target organ [2,19,20]. In the contraction stage, it undergoes activation-induced cell death (AICD) [21], activated cell-autonomous death (ACAD) [21], or memory formation, giving rise to a memory agent in the final case [22]. After the contraction stage, an agent representing a memory CD4+ T lymphocyte behaves as a naive agent and goes through the same three stages, albeit with different activation conditions and locations. At the molecular scale, signal transduction and gene regulation in each agent are described by a logical model consisting of 73 Boolean variables representing proteins and genes and 156 pairwise interactions between these variables. This model is an extended version of a previously published model used to study the plasticity of CD4+ T cell differentiation [2]. Based on the generic human metabolic model Recon 2.2.05 [23], together with 159 microarray datasets and 20 proteomic datasets, five phenotype-specific constraint-based models were built. The metabolic networks of effector Th0, Th1, Th2, Th17, and Treg CD4+ T cells comprise 4234, 4160, 4674, 5223, and 3854 reaction fluxes, respectively, and 2452, 2377, 2623, 2833, and 2178 metabolites, respectively. Cytokine concentrations are used as inputs for the cellular and molecular models; they are conversely influenced by the cellular models. The models at the molecular scale determine most of the attributes of the cellular agents.
Different parameter sets can be used to model and simulate diverse immune phenomena. At the highest level, the simulation algorithm is a series of fixed time steps, the number and duration of which are determined by the duration and resolution of the modeled event. Fig 2B illustrates the information flow between the four components of the multi-scale model. Every time step, each agent is evaluated for all behaviors according to the relevant logical model outputs, such as activation (naive agents only), cell death, memory formation, migration, and cytokine secretion. Its progression through the mammalian cell cycle is dependent on various metabolic/synthesis rates. For each agent, the logical model is implemented by using information from the agent and the compartment it is in to parameterize the input components, then calculating the activity level of each non-input component at convergence (estimated by the Kolmogorov-Smirnov test [24]), and finally using the outputs to update the agent's attributes and inform their behaviors. Some of the outputs are used to parameterize the constraint-based models in order to update the agent's metabolic/synthesis rates. Cytokine concentrations are then calculated by considering the attributes of every agent. A more detailed account can be found in S1 Text.

The model reproduced observed responses of naive CD4+ T cells to different combinations of cytokines
Eizenberg-Magar et al. [25] used different combinations and dosages of cytokines to differentiate naive CD4+ T cells into effector cells of various phenotypes. IL12 produced the Th1 agents themselves reside in and travel between three isotropic compartments. Each compartment is associated with a set of cytokine concentrations, modeled by ordinary differential equations, that affect and are affected by the agents.
https://doi.org/10.1371/journal.pcbi.1009209.g001 PLOS COMPUTATIONAL BIOLOGY phenotype, IL2 and IL4 together produced the Th2 phenotype, IL6 and TGFβ together produced the Th17 phenotype, and finally, IL2 and TGFβ together resulted in a Treg phenotype. A summary of this study can be found in S1 Text. We validated the multi-scale model by reproducing these results. We conducted four Monte Carlo simulations, each comprising 50 realizations of the modeled system (standard specifications, see Methods). In each case, the initial cytokine concentrations were configured to match a set of concentrations used in the experimental study. All the other parameters were set according to the standard specifications. We also conducted a fifth simulation with a combination of all five cytokines plus IFNγ. For each simulation, the 50 realizations were averaged to generate the results presented in Fig 3A. Every reported cell or event count is the average value from an ensemble of 50 realizations. Therefore, despite the discrete nature of the agents, we report continuous values.
As illustrated in Fig 3A, the observed differentiation patterns, as summarized in the last paragraph, were reproduced in silico. Such differentiation patterns were considered for non-physiological cytokine concentrations at the single-cell level in our previous study [2]. The multiscale model replicated our previous results for physiologically relevant cytokine concentrations at the cell-population level. The sensitivity of the multi-scale model varies for different cytokine combinations. For example, IL2 and IL4 triggered a strong response, converting all naive CD4+ T cells into a uniform population of Th2 effector cells. In contrast, the Treg response to IL2 and TGFβ was much weaker, with most effector cells adopting the Th0 phenotype.  The flow of information between the models during simulations. The intracellular events described by the logical model and the metabolic models are assumed to occur instantaneously. The cytokine concentrations described by the ordinary differential equations are updated every hour, which is the antigen load's time-resolution. The agent attributes are also updated every hour, but some agent-associated events take multiple hours to complete, such as cell cycling. Red numbers indicate the sequential steps of a complete update at the systemic/cellular time scale. Labels on the arrows list the variables exchanged between models.
https://doi.org/10.1371/journal.pcbi.1009209.g002 where 50 realizations of the modeled system were generated and averaged, thus explaining the continuous values (see Methods and S1 Text for details). The peak is defined as the time point with the maximum number of effectors and reactivated memory cells. In the column labels, 0 represents Th0; 1, Th1; 2, Th2; S, Th17; and R, Treg. A combination of symbols represents a mixed phenotype; for example, 12 represents the Th1-Th2 hybrid phenotype. (a) Phenotypic distributions of CD4+ T cells in different cytokine media. The initial cytokine concentrations were set to match an experimental study [25]. The first row represents the case with IL2, IL4, IL6, IL12, IFNγ, and TGFβ. (b) Phenotypic distributions of CD4+ T cells given different initial concentrations of IL6 (multiples of 20 ng/mL). In each simulation, the initial concentration of TGFβ was set to 10 ng/mL, while the other cytokine concentrations were set to zero. S:R denotes the ratio of Th17 count to Treg count (c) Same as (b), except the initial concentration of TGFβ was set to 100 ng/mL. (d) Total counts of naive cell activation and maximum cell counts given different initial concentrations of IL2 (multiples of 5 ng/mL). The initial concentrations of the remaining cytokines were set to zero.

PLOS COMPUTATIONAL BIOLOGY
https://doi.org/10.1371/journal.pcbi.1009209.g003 PLOS COMPUTATIONAL BIOLOGY cytokines. It presents the activity levels of selected nodes in estimated attractors, which are consistent with the results in Fig 3A. For example, GATA3 (activity level of GATA3) is 1 (maximal) for the IL2-IL4 combination, thus explaining the strong Th2 response presented in Fig  3A. Details regarding how these estimated attractors were computed can be found in S1 Text.
During our simulation corresponding to the IL6-TGFβ combination, both Th17 and Treg phenotypes were adopted by effector cells. Although TGFβ alone induced T cell differentiation into the Treg lineage, when IL6 was also present, T cells differentiated into Th17 or Treg in an IL6 concentration-dependent manner. The simulation with IL6 and TGFβ was repeated with multiples of the original IL6 concentration (20 ng/mL) used by Eizenberg-Magar et al. [25]. The results show a positive correlation between the Th17:Treg ratio and IL6 concentration (Fig 3B), in agreement with our previous dose-response study [2]. Repeating the simulation with a tenfold higher concentration of TGFβ resulted in higher cell counts but the same general trend (Fig 3C).

The model reproduced IL2-mediated metabolic regulation of CD4+ T cells
IL2 is a key cytokine in the metabolic regulation of CD4+ T cells because it induces both their proliferation and activation-induced cell death [26]. Therefore, this validation study centered on the regulatory effects of IL2 on the expansion stage of the cell. Monte Carlo simulations were carried out with different initial IL2 concentrations (0, 5, 50, 500, and 5000 ng/mL, multiples of the concentration used in experiments [25]). All other initial cytokine concentrations were set to zero. The other parameters were set according to the standard specifications (Methods). We calculated two metrics for each of the five Monte Carlo simulations. First, we calculated the average number of naive cells activated at the peak of each simulation, with the peak meaning the time point with the highest average number of effector and reactivated memory cells. Second, we calculated the average number of effector and reactivated memory cells at the peak of each simulation. Henceforth, we will call the two metrics 'number of activation events' and 'peak cell count'. As explained in the Methods, our agents representing effector cells and reactivated memory cells behave in the same way.
The results ( Fig 3D) show that with no IL2 in the compartments, 116 activation events (out of 433 potential activation events, one for each naive T cell) had occurred when the simulated immune response peaked, but 229 effector cells and reactivated memory cells were present at the peak, suggesting a basal level of proliferation even in the absence of IL2. At the highest IL2 concentration, the maximum number of cells (674) is almost three times that in the absence of IL2 (229). These results can be compared to the experimental results in a study [27] where CD4+ T lymphocytes were stimulated with anti-CD3 and IL2 at various concentrations (0.625 to 25 ng/mL). At the highest IL2 concentration, the extent of DNA replication in the cells undergoing the first round of division (a proxy for cell count) was measured at almost three times that at the lowest IL2 concentration, in agreement with our simulation results.
The population-level results shown in Fig 3D can be explained by considering the models embedded in each agent at the molecular scale. Fig 4C shows the activity levels of the metabolic nodes in the estimated attractors corresponding to three initial IL2 concentrations (0, 5, and 50 ng/mL). The results indicate that the activity levels of metabolic nodes should either increase or stay constant in response to higher IL2 concentrations. Details regarding how these estimated attractors were computed are presented in S1 Text. The activity levels of metabolic nodes in each estimated attractor were used to alter the constraints on the metabolic fluxes in the Th0 metabolic model before the fluxes were optimized for both the G1 and S phases of the cell cycle. Mathematically and computationally, it means that an objective function representing the goal of each phase was defined by combining the fluxes, and it was maximized by finding an optimal set of fluxes under the aforementioned constraints (linear programming problem solved by the simplex algorithm [23]). In G1, the objective function to be maximized is the rate of biomass synthesis, excluding DNA. In the S phase, it is the rate of DNA synthesis. The results are shown in Fig 5A-5D. Each scatter plot shows multiple fluxes off the identity diagonal, indicating changes in flux values. They suggest that IL2 is a metabolic regulator. For example, it boosts the biosynthetic fluxes related to aerobic glycolysis and lipid synthesis. IL2 not only enhances or reduces metabolic reactions, but it also turns them on or off as attested by the many fluxes lying on the axes. Some fluxes even lie on the opposite diagonal, indicating reversals of the corresponding reactions. Biologically, these reversals point to significant effects of IL2 on metabolism. The reversal of a metabolic flux means a catabolic reaction is turned into an anabolic reaction.

The model reproduced the population dynamics of CD4+ T cells during influenza infection
We used the standard specifications of the model to simulate the CD4+ T cell response to influenza infection. Due to the high computational costs associated with linear programming

PLOS COMPUTATIONAL BIOLOGY
and exponential population growth, the initial number of naive agents in the draining lymph node was set to 10 to represent a portion of the node only. The initial cell counts in the other two compartments were set to respect the relative population sizes in the three compartments. There were no more sophisticated reasons behind the choice of 10. As for the influenza infection, the viral trajectory from an experimental study [20] was converted into the input signal ( Fig 6A, the curve labeled 0.8) as described in the Methods and, with more detail, in S1 Text. The signal ranges from 0 to 1 in arbitrary units. Biologically, it represents the viral load and the PLOS COMPUTATIONAL BIOLOGY activity of the immune system triggered by the virus and relevant to CD4+ T cell activation. Computationally, it is the probability of its stimulating a CD4+ T cell agent.
The simulated population dynamics are shown in Fig 6B. In the first four days post-infection, the trajectory representing the draining lymph node rises nearly five-fold (10 to almost 50). In comparison, over the same period, the experimental study reported a 16-fold expansion of the corresponding population in the draining mediastinal lymph node. Although the rise in cell count in Fig 6B is smaller than the observed expansion, the same qualitative trend is present in both sets of results. The simulations also indicate that the CD4+ T cell population in the lungs should rise more slowly than that in the lymphoid tissues: a third of the peak after four days of infection. The predicted time delay is in qualitative agreement with the experimental study, where no CD4+ T cells were detected in any lung tissues at four days post-infection. In the experimental study, all cell populations were observed to peak at six days post-infection and decline between six and eight days post-infection, again agreeing with Fig 6B.

The model predicted switch-like and oscillatory emergent behaviors of CD4+ T cells in response to infection
After validating the multi-scale model, we performed further simulations to test its ability to make predictions. To this end, the default input signal was modified in terms of its peak strength and long-term behavior, as explained in S1 Text. In other words, we changed the signal quantitatively and qualitatively in order to test the model's robustness.
The first robustness test centered on the input signal strength. We generated five input signals to model acute infections with different peak strengths (Fig 6A). The Monte Carlo simulation presented in the previous subsection was repeated with each signal. Fig 6C and 6D summarize the results of averaging the six sets of 50 realizations. During the simulations, the input signals peaking at 0.05, 0.2, and 0.4 were too weak to activate any naive cells. In the model, a naive cell must migrate to the draining lymph node before activation can occur, but it can only move and stay there in response to a strong input signal, effectively forming a filter for weak input signals. On the other hand, when the input signal was set to peak at or above 0.6, the effector cells amplified its effects on proliferation nonlinearly. In Fig 6A, as the input signal peak rises from 0.6 to 0.8 to 0.95 (less than two-fold), the number of activation events increases from 4.60 to 118.46 to 318.58 (almost 70-fold), and the maximum cell count increases from 7.70 to 831.84 to 14705 (over 1900-fold). As before, each number of activation events or maximum cell count is an average over the ensemble of 50 realizations, thus explaining the continuous values. In the model, each effector cell's activation level is an attribute reflecting the input signal, and it is passed to the nonlinear logical model, which ultimately regulates proliferation. Cooperatively, the effector cells secrete IL2 to boost further proliferation in the population, so as they proliferate, the rate of proliferation increases. This effect is enhanced by the simple fact that each cell produces multiple daughter cells. Taken together, our results suggest that CD4+ T cells demonstrate a switch-like behavior in response to infection.
The second robustness test centered on the input signal's pattern. To simulate chronic infection, we produced a steady input signal, which does not decrease after peaking (Fig 6E). Unlike the acute input signal (Fig 6A), which led to a decrease in the simulated immune response (Fig 6B) after the signal was removed, the chronic input signal produced oscillations in the population size of CD4+ T cells (Fig 6F) despite the fact that the cells were not synchronized. This behavior emerges from the interplay between persistent stimulation and episodic replenishment of naive CD4+ T cells. In the presence of persistent stimulation, naive cells are continuously activated into effector cells. As more effector cells enter the contraction stage, cell death occurs more and more frequently, while the supply of naive cells continues to drop, slowing down the growth of the effector cell population. Eventually, their death rate exceeds their replacement rate, leading to a decline in their population. When the naive cell population

PLOS COMPUTATIONAL BIOLOGY
drops below a threshold, a fixed number of naive cells are added to the compartments. Although naive cells are replenished continuously in a lifelong process [28], an episodic mechanism is representative of the model's context presented in this work: when the thymus transiently produces more naive T cells during viral infections [29].
Oscillations can also be seen in the concentrations of secreted cytokines (Fig 7). Comparison with the first row of Fig 4A explains the variations between the six sets of concentration dynamics. IL6 is extremely low in the relevant estimated attractor. This, together with the fact that the cytokine-secreting ability of an effector CD4+ T cell increases with its division count [30,31], explains why the target organ concentration is higher than the lymph node concentration in Fig 7C (more divided cells in the target organ and low secretion rate in the lymph node). In the other panels of Fig 7, the lymph node concentrations are higher because of the lymph node's small volume, while the activity levels of the relevant cytokine nodes are not low enough to offset this effect. Fig 7A seems to contradict IL2 being zero in the estimated attractor, but the relatively plastic effector cells in the lymph node secrete IL2 by default [30,31]. Finally, the concentrations of IL17 (Fig 7D) and IL21 (Fig 7E) are an order of magnitude or two lower than the rest but fluctuate more. It is because their initial concentrations were set to zero. However, the oscillations in cytokine concentrations did not cause the cell count oscillations shown in Fig 6F. As stated above, the chronic input signal interplaying with episodic replenishment of CD4+ T cells caused the cell count oscillations. In turn, the oscillating population of cytokine-secreting cells caused the oscillations in cytokine concentrations.

Sensitivity analysis and the generality of results
To ensure that the results reported in the previous sections reflect the behavior of our computational platform, where information is passed between the constituting models, rather than just one parameter set, a sensitivity analysis of the model's parameters was performed. Given the large number of parameters, a global sensitivity analysis was infeasible. Instead, we performed a local analysis focused on a few key parameters. We deemed the initial number of agents in each compartment highly important because a small cell count may artificially inflate random events' impact. We also chose to assess the time interval between updates (step_size) due to the multiple time scales used in the framework. Due to the importance of the cells' migration rates in producing the switch-like behavior, we tested the probability of an effector agent migrating from the circulatory system to the target organ (k_effector_CtoT), as well as the fold increase in the probability of a naive agent migrating from the target organ to the lymph node when there is an infection (k_toln_boost).
Initial number of agents in each compartment. In a simulation using the default parameters, there are only 10 naive agents in the lymph node compartment initially, so the small cell count may inflate the impact of stochastic events spuriously, leading to drastic results such as the switch in Fig 6C and 6D. Using 10 times the initial numbers of agents in the three compartments, we repeated the simulations illustrated by those figures. The results are summarized in Table K in S1 Text. As expected, more active CD4+ T cells were triggered when there were more naive cells at the start of the simulations. However, the switch-like behavior appeared at the same threshold.
Time interval of each update (step_size). The default value of step_size is one hour. As explained in the caption of Fig 2, the intracellular events are assumed to occur instantaneously, while cellular events such as cell death may take several hours. To justify this separation of time scales, this parameter should not be set below the default. In principle, the user can experiment with increasing step_size without violating this principle because the algorithm will adjust the other parameters accordingly. For example, the probability that an effector cell migrates from the circulatory system to the target organ (k_effector_CtoT) is automatically doubled when step_size is changed by the user to two hours, thereby preserving the migration pattern.

PLOS COMPUTATIONAL BIOLOGY
We carried out a sensitivity analysis by repeating twice the simulation presented in Fig 6B after setting step_size to two and 0.5 hours, respectively. We present the results in Figs D and E in S1 Text. In agreement with Fig 6B, the population dynamics peaked at around six days postinfection in all simulations. The steep decline between six and eight days post-infection also occurred in all simulations. Crucially, the time delay of the target-organ dynamics relative to the lymph-node dynamics appeared in all simulations. However, the simulations activated different numbers of cells. For example, if step_size is one hour, meaning it takes 400 steps to complete the simulation, a naive cell will need 24 stimuli to activate. Assuming the probability of receiving one stimulus at any step is 0.1, we can compute from the binomial distribution a probability of 1.40e -3 that a naive cell will activate in such a simulation. If step_size is doubled, it will take 200 steps to complete the simulation, and because a step is longer, a naive cell will only need 12 stimuli to activate. Assuming the probability of being stimulated at any step is still 0.1, the probability that a naive cell will activate is higher, at 1.53e -2 . If step_size is halved, the probability that a naive cell will activate drops to 1.66e -5 .
Effector agents' probability of migrating from the circulatory system to the target organ (k_effector_CtoT). We increased and decreased this parameter by an order of magnitude and, in each case, we repeated the simulations illustrated by Fig 6C and 6D. The results are shown in Table L in S1 Text. The switch-like behavior appeared at the same threshold as the default case in both sets of simulations. However, as the migration probability was increased, the simulated immune response became weaker. This can be explained in terms of AICD. In the simulations with an increased k_effector_CtoT, more effector cells gained entry into the target organ per time step, thus increasing the frequency of restimulation and AICD. Fig F in S1 Text shows the population dynamics of effector cells in the two simulations with the default input signal. Similarly to the default case, the population dynamics peaked at around six days post-infection in both simulations and declined steeply between six and eight days post-infection. Somewhat surprisingly, the time delay of the target-organ dynamics relative to the lymph-node dynamics appeared in the new simulations despite the changes in this migration probability. From this result, the mechanism causing the time delay can be identified: it lies in the transition from the lymph node to the circulatory system, not that from the circulatory system to the target organ.
Naive agents' probability of migrating from the target organ to the lymph node when there is an infection (k_toln_boost). We halved and doubled this parameter and in each case, we repeated the simulations illustrated by Fig 6C and 6D. The results are shown in Table M in S1 Text. Contrary to our expectation-a decrease and increase in cell migration into the lymph node should make the switch harder and easier to flip, respectively-the switch-like behavior appeared at the same threshold as the default case in both sets of simulations. Based on this result, the underlying mechanism can be identified: the key is the inflamed lymph node's ability to retain naive cells, not its ability to attract them. Fig G in S1 Text shows the population dynamics of effector cells in the two simulations with the default input signal. These results are in agreement with the default case, indicating that the migration of naive cells from the target organ into the lymph node is not the reason for the time delay.
In summary, the results reported throughout this section are robust with respect to parametric sensitivity. They do reflect the multi-scale model's general behavior in passing information between the constituting models. In a model like ours, it is possible for a small cell count to inflate random events' impact. However, when we increased the initial number of naive cells in each compartment, the switch-like behavior appeared at the same threshold. The presence of multiple time scales means step_size is critical. The major trends in the population dynamics of CD4+ T cells during influenza infection, including the time delay of the targetorgan dynamics relative to the lymph-node dynamics, were conserved as we varied this parameter. The switch-like behavior and time delay were also successfully reproduced after changing effector agents' probability of migrating from the circulatory system to the target organ and their probability of migrating from the target organ to the lymph node upon infection. This result confirms the robustness of our results. In addition, the lack of sensitivity also identifies the mechanisms responsible for the results by elimination. Since changing naive agents' migration probabilities did not affect the switch-like behavior in the simulations, we can now rule out the inflamed lymph node's ability to attract naive cells as the cause. A likely mechanism is its ability to retain them. Neither of the two changes affected the time delay of the target-organ dynamics relative to the lymph-node dynamics. We can therefore rule out the movement of effector cells into the target organ as the cause, as well as the entry of naive cells into the lymph node. The most likely reason is the exit of effector cells from the lymph node into the target organ.

Discussion
As reviewed in the introduction, numerous computational approaches have been used to model the immune system and its functions, particularly in relation to T cells. Although impressive and informative, those studies did not employ together methods such as logical modeling and genome-scale metabolic models. Furthermore, many mechanistic details about the metabolic regulation and plasticity of CD4+ T cells were not considered. As such, they could not capture essential nonlinear aspects of immune responses, such as emergent properties. Therefore, a more complex model is an important addition to the systems immunology toolkit. Many models are phenotypic in that they lack precise molecular and cellular details, so they cannot readily be informed by-and their outputs compared with-experimental and patient data [32]. Although phenomenological models have been used to explain trends in population dynamics, such as an influenza infection kinetic model [33], mechanistic models can capture stochastic effects, spatial heterogeneity, and cause-effect relations at the cellular and molecular scales, such as a hybrid ABM-ODE model identifying TGF-β1 as a major inhibitor of effector T cells [34]. Some modeling studies used detailed agents but each focused on a given compartment such as germinal centers [35]. A few agent-based models were built by including different immune cell types in multi-compartment setups not dissimilar to the one presented in this paper [36,37]. However, the agent behaviors were not linked through the use of adaptive, intracellular models of biochemical networks. While they made the individual agents adapt to the changing cell populations and cytokine concentrations, they did not capture the resulting feedback control with the overall models. Logical models informing the phenotype of CD4 T cells have been developed [38] and even been used to build multi-cellular models [39]. However, if used for our study, the lack of differentiated anatomical compartments would not allow time delays caused by cell migration at the systemic scale to be captured, such as the time delay in Fig 6B. Therefore, their simulations would not display the switch-like behavior which depends on the fact that CD4+ T cells must migrate to the draining lymph node before being activated. Finally, while metabolism has been shown to be essential for T cell activation [40] and genome-scale models of CD4+ T cells exist [6], they were not used in any dynamical multi-scale modeling studies.
To the best of our knowledge, our model is the first to integrate ABM, ODEs, GSMMs, and a logical model within a compartmental setup. In addition to tackling the drawbacks outlined above, the multi-approach strategy confers the advantage of modularity on the model. New agent types can be developed to represent other cell types (such as B cells), and integrated with the existing framework. This was how we expanded our existing logical model of CD4+ T cell differentiation [2], and integrated it with the framework to reproduce simulation results (at the single-cell level and for abstract cytokine levels) at the population level using physiological cytokine concentrations. Modularity also permits reusability of existing models; new versions of intracellular models, using similar or different modeling approaches, can be integrated with the framework with minimal effort. In addition to the gains in computational performance brought by using different levels of granularity for different parts of the model, one could decouple the parts and cache simulation results. For example, we built a library of steady states for our metabolic models to match different estimated attractors of the logical model. Finally, modularity also helps with data integration. For instance, if a dataset related to memory formation is available at the cellular level, while another dataset about metabolism is available at the molecular level, one can easily work with both datasets within our framework. To facilitate validation at the modular-and systems-level-an important step after integrating a new agent, replacing an intracellular model with another, or integrating data-we have compiled the experimentally observed trends used for validation in Table N in S1 Text.
The multi-scale nature of our model is very important because immune responses are multi-scale. The concentrations of cytokines at the systemic scale, the collective behavior of the cell population, and the heterogeneity within the population are interwoven due to nonlinear structures such as feedback loops. For example, the difference between our results in Fig 3D and an experimental study [27] seems to indicate that the cells in the experiment were more sensitive to IL2 than the agents in our model. However, the effects of activation-induced cell death reported in [41] were not measured in the experiment which only considered the first round of division. In contrast, our model captures the complex relationship between IL2 (systemic scale), cell expansion (cellular scale), and AICD (molecular scale): IL2 stimulates both cell expansion and death. The dampening effects of AICD on cell expansion can be seen in the saturating behavior illustrated by Fig 3D. At the systemic level, compartmental ordinary differential equations were chosen to model cytokine dynamics since the compartments are considered isotropic and molecules are indistinguishable. Because of the very large number of molecules, they can be solved using deterministic algorithms. At the cellular level, an agent-based approach was selected to serve the needs to model heterogeneity in the population and track each cell's unique combination of attributes. While ordinary differential equations can describe how the sizes of several populations evolve, even cells within the same population (naive cells, for instance) have different attributes (cell cycle stage, for instance). We decided that cell expansion should be treated as a process (mammalian cell cycle) rather than a single event (division) to provide realistic links between the agent-based model and the models at the molecular level. In our simulations, for most agent behaviors (such as differentiation), the activity levels of the logical model were normalized before being passed to the agent-based models as probabilities. On the contrary, the activities of the nodes interfacing with the genome-scale metabolic models were used to alter the constraints on metabolic fluxes before the fluxes were optimized. A constraint-based approach was chosen over ordinary differential equations for the metabolic models to circumvent the lack of kinetic parameters and to allow a comprehensive (genome-scale) representation of metabolism. A disadvantage is the lack of transient dynamics. However, our algorithm mitigates the problem by updating the steady-state for each cell at every time step. Because this solution comes with a high computational cost, we built a library of optimized fluxes corresponding to the estimates of common logical model attractors.
Our model was able to reproduce diverse experimental observations made on CD4+ T cells, so we are confident in its validity. The model reproduced the ability of CD4+ T cells to differentiate into different phenotypes in a context-dependent manner, as well as the regulation of CD4+ T cell metabolism by IL2. We also reproduced the population dynamics of CD4+ T cells during influenza infection. These experimental observations were made in completely independent studies, and we emphasize that, except for the cytokine concentrations, they were not used to parameterize the model. The model was also able to make new predictions, as demonstrated by the simulation of emergent behaviors in virtual CD4+ T cell populations. The switch-like and oscillatory behaviors arise from interactions among biological phenomena taking place at different spatial scales. The switch-like behavior might help the body prevent chronic inflammation in the continuous presence of a small and non-threatening amount of antigen. In fact, switch-like behaviors have been observed in populations comprising other immune cell types. For instance, the chronic inflammation triggered by mast cells is a bistable switch [42]. Similarly, oscillations are not uncommon in the immune system [43], possibly as side effects of homeostasis. However, no one has observed or predicted such nonlinear behaviors in relation to CD4+ T cells. Although the predictions have not yet been validated experimentally, they highlight our model's ability to tease out new knowledge from the complexities of the immune system. Through a series of sensitivity studies, we confirmed the robustness of the switch-like behavior and identified the reason for its occurrence to be the ability of the inflamed lymph node to retain naive cells. At the same time, we also confirmed the robustness of the target-organ dynamics' time delay relative to the lymph-node dynamics, and identified the key migration step behind it as the exit of effector cells from the lymph node into the circulatory system. These results stress the ability of our model to identify the mechanisms underlying emergent behaviors.
Several opportunities to expand the presented modeling framework exist. The choice of a generic signal to represent pathogens/antigens and the rest of the immune system, with no feedback from the triggered response, can be improved upon by adding more mechanistic details. These unobserved components may also be a cause for non-identifiability concerns. When different parameter sets result in the same phenotype, a user might not distinguish redundancies intrinsic to CD4+ T cells from the unobserved components' contributions. Another possible consideration is the use of well-stirred compartments due to the computational costs involved in solving a cell-based model considering cytokine diffusion and physical cell-cell interactions within each compartment. In our simulations, the agents demonstrated different sensitivity levels to different cytokine combinations, in contrast to the experimental evidence indicating a clear relation between the cytokine combinations and phenotypes [25]. With the sole exception where only IL2 and IL4 were in the system, many effector cell agents exhibited the Th0 phenotype, as summarized in Fig 3A. This inconsistency is, however, a matter of interpretation. We used the low activity levels of transcription factor nodes in Fig 4A as the probabilities that an effector cell would adopt the corresponding phenotypes. Considering the qualitative nature of logical model outputs, we can simply rescale them to enhance sensitivity. Each cytokine concentration is converted to a logical model input using the same parameters regardless of the cytokine type. Using different parameters may also resolve the inconsistency. Another exciting area for further study is the negative feedback loop regulating the replenishment of naive cells. In addition to building a more elaborate model of the control system-perhaps by adding the effects of IL7 [29]-one could contrast the outcomes of using different replenishment thresholds and prolonging each replenishment, as these are the parameters crucial for the oscillations. As explained in the Methods, our agents representing effector cells and reactivated memory cells behave in the same way, so another extension is to confer different properties (such as proliferation rate) on the two groups, and then contrast the emergent behaviors of effector cells activated from replenished naive cells and reactivated memory cells.
We believe our framework provides the basis of a virtual immune system, well-suited to exploring the fundamental properties of the immune system and to modeling immune responses to specific diseases. Its multi-scale nature allows researchers to address cellular and molecular immunological problems with a systems approach by unraveling counterintuitive, nonlinear properties in the immune system.

Structure of the multi-scale model
Systemic scale. The multi-scale model is divided into three isotropic compartments. In each compartment, IL2, IL4, IL6, IL12, IL17, IL18, IL21, IL23, IL27, IFNγ, and TGFβ are represented as time-dependent homogeneous concentrations. An agent representing a CD4+ T cell can sense all cytokines except IL17 and IL21, and it can secrete IFNγ, IL2, IL4, IL6, IL17, and IL21. The concentrations of the former are passed to the models at the cellular and molecular scales (see below), while the concentrations of the latter are influenced by those models. This setup is consistent with our previous study at the molecular scale [2]. Furthermore, in each compartment, a study-dependent and time-varying input signal represents the antigen load and the activity of the rest of the immune system. Cytokine concentrations are governed by 33 ordinary differential equations: C i, j is the concentration of cytokine j in compartment i (M), P i, j is the baseline production rate (M�h -1 ), k deg i, j is the degradation rate constant (h -1 ), s t, i is the input signal (antigen load and the rest of the immune system, dimensionless) in the compartment at time t (hours), P in i, j is the production rate in response to the input signal (M�h -1 ), N eff i, j is the number of cells/agents producing the cytokine in the compartment, P eff i, j is the rate of production by one cell/agent in the compartment (M�h -1 ), and R i is a transport term (M�h -1 ). Compartment 1 is the target organ, 2 the lymphoid tissues, and 3 the circulatory system. The transport terms are described by the following equations: and where V 1 , V 2 , and V 3 are the compartments' volumes in dm 3 , Q a is the rate of blood flow into the target organ (dm 3 �h -1 ), and Q b is the rate of lymph flow into the circulatory system (dm 3 �h -1 ). Cellular scale. Cells are represented as discrete and autonomous agents defined by a set of attributes and behaviors. These agents respond to changes in both the cytokine concentrations and the input signal. A newly created agent represents a naive CD4+ T cell. They are created in all three compartments. Then, they undergo three major life cycle stages: activation, expansion, and contraction. An important point is that the agents are not synchronized, so any oscillatory behaviors must be explained in terms of driving forces external to the agents. For example, when the naive cell population drops below a threshold, a fixed number of cells are added to the compartments.
During the activation stage, a naive agent exhibits three behaviors. First, in response to the input signal in the lymphoid tissues (s t,2 ), the naive agent migrates from the target organ to the lymphoid tissues at a higher rate and stays there. Second, s t,2 is the probability that the naive agent is stimulated in a time step, and if it gets stimulated, its TCR_stimuli increases by one, and it records the signal magnitude (between zero and one) and cytokine concentrations. Third, if the naive agent has been stimulated a number of times (act_thres), it is activated and becomes an effector agent. Computationally, activation involves using the logical model and appropriate metabolic models to determine the effector agent's attributes, such as its phenotype (Th0, Th1, Th2, Th17, Treg, and hybrid phenotypes comprising these five states).
During the expansion stage, an effector agent exhibits four behaviors. First, it continuously monitors the environment for changes in the input signal and cytokine concentrations and updates internal attributes, such as the amount of biomass, which are used to track its progress through the cell cycle. Second, at the end of the cell cycle, it divides into two identical copies. Third, it secretes cytokines in accordance with its changing cellular phenotype. Fourth, after several rounds of division, it leaves the lymphoid tissues, enters the circulatory system, migrates to the target organ, and stays there.
During the contraction stage, an effector agent exhibits four behaviors. First, if it is re-stimulated by the input signal, it may undergo activation-induced cell death (AICD). Second, in the absence of restimulation, it may become a memory agent. Third, in the absence of restimulation and memory formation, it may undergo activated cell-autonomous death (ACAD). Fourth, if none of the first three cases applies, its activation level (in terms of TCR_strength and CD28_strength) will decrease by an amount determined by the input signal's strength in the time step.
After the contraction stage, a memory agent undergoes the same three stages described above for naive agents. However, fewer stimuli are required for its activation, and the resulting activation level is higher [44,45]. Unlike a naive agent, a memory agent can get stimulated and activated in all three compartments. Upon activation, it becomes a reactivated memory agent and behaves exactly like an effector agent [46]. Unless otherwise stated, all properties of effector cells and effector agents apply to their reactivated counterparts too.
Molecular scale. Most of the attributes of each agent are determined by a logical model of signal transduction and gene regulation and constraint-based models of metabolism.
Signal transduction pathways and gene regulation processes captured in the model include the canonical pathways regulating the differentiation of CD4+ T cells into various cell phenotypes [2] and events such as aerobic glycolysis, memory formation, and cell death. These pathways are detailed in S1 Text. The network diagram of the logical model is presented in Fig  A in S1 Text. Its structure is detailed in S1 Text. The model itself is accessible in Cell Collective [47], where it can be downloaded as an SBML file or directly modified and implemented in simulations.
The logical model has three classes of inputs corresponding to 1) cytokine concentrations, provided by the ordinary differential equations, 2) receptor signaling (TCR and CD28), determined by the agent's activation level, and 3) agent attributes transformed into logical model inputs, such as the activity of mTORC1, the concentration ratio of AMP and ATP, abundance of ribosome, and the agent's resistance to apoptosis.
The outputs of the logical model affect agent attributes and behaviors, such as phenotypes (Th0, Th1, Th2, Th17, and Treg), apoptosis, and cytokine secretion. They also parameterize the metabolic models. There are four classes of outputs representing 1) the transcription factors Tbet, GATA3, RORgt, and Foxp3, 2) the cytokines produced by the agent, 3) seven metabolic components providing the interface with the genome-scale, constraint-based models of metabolism, and 4) six cellular events (cell cycle progression, autophagy, memory formation, ACAD, AICD via the Fas pathway, and AICD via the B-cell lymphoma 2 (BCL2) pathway).
The metabolic networks of effector Th0, Th1, Th2, Th17, and Treg CD4+ T cells are formalized in five phenotype-specific constraint-based models. 159 microarray datasets were integrated with 20 proteomic datasets to define the gene activity profiles (matrices) of the five phenotypes. The profiles and the generic human metabolic model Recon 2.2.05 [23] were used to build the constraint-based models by the GIMME method in the COBRA toolbox [48]. In summary, the Th0 model contains 4234 metabolic fluxes and 2452 metabolites; Th1, 4160 and 2377; Th2, 4674 and 2623; Th17, 5223, and 2833; and Treg, 3854 and 2178.
During a simulation, at each time step at the cellular and systemic scales, multiple simulations at the molecular scale (logical and constraint-based models) are run. At each discrete update of the logical model, the activation probability of each input component serves as an input to the logical model. Some of the outputs (seven metabolic components) are used to set the constraints in the metabolic models, which are then used to calculate the metabolic/synthesis rates in the agents. A detailed account of the simulation algorithm and the corresponding pseudocode can be found in S1 Text.

Variables and parameters of the multi-scale model
The multi-scale model has four classes of variables and parameters. They are agent attributes, variables/parameters specific to the compartments, settings for the numerical methods, and four parameters that can be varied for model calibration. All the variables and parameters are presented in Tables A-J in S1 Text. A detailed description of all four classes and the parametric estimation process is in S1 Text, where it is also explained how the model can be used to study diverse immune phenomena by changing the parameters. Unlike most multi-scale models, where the lack of datasets typically makes it hard to constrain the many parameters uniquely, most of the parameters are-or can be-derived from literature values, while the remaining uncertainty is lumped into the four parameters for model calibration.

Convergence study
Because the model contains stochastic elements, such as the logical model and most agent behaviors, multiple simulations (realizations of the model) must be performed to obtain robust averages of model outputs. A convergence study helps the modeler decide how many simulations are required to guarantee robustness. This step is necessary, following parameterization, to describe a specific immune phenomenon. The convergence study we ran before the validation and robustness studies presented in this paper, as carried out using the default, acuteinfection-like input signal (Fig 6A, curve 0.8), is described in S1 Text. Based on the results, we decided that 50 simulations were required for the average model outputs to converge. The specifications detailed in S1 Text, namely the input signal, initial cytokine concentrations, cytokine production rates associated with the input signal, and the need for 50 model realizations, form the standard specifications of the studies presented in this paper. In the Results section, it is understood that a simulation refers to a set of 50 realizations of a parameterized model.