A concise network of the main components controlling T cell metabolism was constructed based on experimental information. A central actor in the lymphocyte metabolic activity is the AMPK complex, which is capable of sensing the intracellular AMP/ATP ratio, which represents the T cell energy pool availability, and of regulating the main metabolic pathways leading to the production of energy reserves (OXPHOS) or to the rapid generation of metabolites and structural proteins (glycolysis) (47–52). It has been demonstrated that two main effects related to metabolism take place upon TCR stimulation and CD28 co-stimulation (53). First, an increase in the basal activity of oxidative phosphorylation (OXPHOS) arises promoted by the action of the nutrient sensor AMPK, which is activated directly by the PI3k-AKT axis and calcium release (54). Afterwards, mTORC1 is activated with the consequent inhibition of AMPK and the activation of glycolysis (55).

In the network, AMPK is activated in several ways: signaling from the TCR and CD28 via the PI3k-AKT axis, calcium release, a high AMP/ATP ratio, the activation of the serine–threonine liver kinase B1 (LKB1) and the Foxp3 transcription factor (7, 47, 48, 56). Although there are several signaling intermediates proposed in these pathways, they were mathematically implied due to the lineal nature of the signaling trajectory in order to obtain a set of simplified logical propositions.

The balance between OXPHOS and glycolysis depends primordially on a negative feedback loop between AMPK and mTORC1, which plays the role of a metabolic polarization switch (50, 52, 57, 58). In the network presented here, TCR and CD28 activation, an elevated AMP/ATP ratio, LKB1 and Foxp3 activate AMPK, thus inhibiting the activity of mTORC1, finally leading to OXPHOS and inhibition of glycolysis.

Conversely, under a low AMP/ATP ratio, AMPK is inhibited, which in turns activates mTORC1 function and additional inhibition of AMPK, leading to glycolysis ( Figure 2 ). The metabolic module was then incorporated into the previously reported Boolean network of T CD4 cell activation through links associated to AMPK and mTOR ( Figure 3 ). The network previously constructed encompasses modules corresponding to an activation core derived from TCR and CD28 signaling, a feedback loop for IL-2 production though its high-affinity receptor (CD25), the role of the anergy factor NDRG1 (controlled by AKT), activation of the checkpoint by CTLA-4, and the induction of the effector phenotypes by external cytokines (59, 60). The construction of the network can be consulted in (40).

The interactions involved in the metabolism module were formalized as Boolean propositions as shown in Table 1 . After incorporation into the general T-cell activation network proposed in (40), an exhaustive analysis was performed to analyze its general behavior and congruence with our previously published results ( Supplementary Material 3 ). The robustness of the integrated model was also verified by introducing random noise in the initial states and measuring the distance between transition states and attractors. Robustness was also tested by either inducing perturbations in the network structure by random bit flipping of the Boolean functions or by random permutation of the output values (data not shown). The set of functions that defines the whole activation network are shown in Supplementary Material 1 and the model structure can be consulted in Wolfram Notebook format in https://github.com/DrDavidMM/TCD4cell-activation-model-supplementary-material-.git.

Fuzzy logic is characterized by a graded approach, so that the degree to which an object exhibits a given property is specified by a membership (or characteristic) function with truth values ranging from total falsity, µ[w _k ] = 0, to totally true, µ[w _k ] = 1. Here, w _k denotes a fuzzy logic proposition describing the interactions of node k with the rest of network nodes. By assuming that the state of the regulatory network at time t is described by the set {q1(t),…q _n (t)}, µ[w _k ], may be represented by a sigmoid function with continuous variation in the interval [0,1]:

where w _thr is an activation threshold, hereby considered as w _thr = 1/2, and β is a saturation rate (46, 61).

The Boolean interaction rules W_k [q ₁(t), q ₂(t), … q_n (t)] were translated into fuzzy logic expressions

This procedure can be straightforwardly implemented by replacing the Boolean logic connectors ‘and’, ‘or’ and ‘not’ by their fuzzy counterparts according to the following scheme:

An example of the translation from Boolean to a fuzzy framework is:

Within the continuous scheme, the dynamical behavior is determined by a set of ordinary differential equations describing the temporal change of the activity level of the network components. For the k - th node, this is written as

Where d _k is the decay rate of node k. In this work, we assume that and that the default value of d _k = 1, unless otherwise stated. It is important to notice that in absence of an input (w _k = 0), the activity of node k decays exponentially, that is, q k ∼ e − d k t ; therefore, the parameter τ _k = 1/d_k represents a characteristic expression time, so that d_k >1 (d_k <1) gives rise to a relatively rapid (slow) decay of the activity of the element k of the network. This kind of analysis allows to assess the influence of different time-scales of activity of the elements comprising the regulatory network.

In the continuous scheme, the equilibrium states (attractors) of the system are defined by the condition dq_k /dt = 0. This condition implies that, for a specific set of initial values {q ₁(0),…,q_n (0)}, the system dynamically evolves until reaching steady-state values given by:

This latter expression shows that the resulting set of asymptotic states, { q 1 s t , … q n s t } , is determined, besides the initial values q_k (0), by the actual values of the decay rates d_k . A consequence of the former results is that the emergent behavior of the system may conduce to alternative dynamic patterns depending on the specific values of i ) the set of initial concentrations, dosages, or expression levels, {q_k (0)}, and ii) differences in either stimulation times, τ _stim , or characteristic expression times of the network components given by their decay rates, τ _k ~ 1/d_k (46, 61, 62).

The extent and time span of stimulation of T CD4 cells due to MHC-antigen presentation to TCR and of CD80/86 binding to CD28, was broadly modeled by introducing two input nodes whose priming activity only lasts for a limited lapse of time τ _stim . This was described by means of functions with a step-like behavior associated to an initial and constant avidity strength, A _MHC/A and A _{CD 8086}, suffering an abrupt decay at time t = τ _stim by a factor D _MHC/A and D _{CD 8086}, respectively. The time variations of avidity can be related to changes in the number of TCR-MHC-peptide complexes, the presence of adhesion molecules, TCR internalization or degradation after initial engagement, etc. (10, 62–67). Now, by introducing the step-function H(t − τ _stim ) defined by

then the avidity variations are written as follows:

where D_MHC/A and D_{CD 8086} represent the magnitude of detachment reduction of TCR and CD28 from their ligands for times longer than τ _MHC/A and τ _{CD 8086}, respectively.

This approach does not only allow to analyze the behavior of the cell as a function of the extent of stimulation, but also the description of the competitive action between CD28 and CTLA-4 for binding to CD80/86. This is performed through the downstream interactions of CD28 with CTLA-4 (see Figure 3 ). Upon activation, CTLA-4 may down-regulate the engagement of CD28 with CD80/86. However, its inhibitory capacity depends on its decay rate, d_{CTLA 4}, which should be relatively small (<1), in order to have an expression time long enough to overwhelm the influence of factors that promote the transcription factors activity. In the model, longer interaction times of the antigen with the TCR and CD28 allows a sustained activation state before being arrested due to the activity of CTLA-4. This process is simulated by introducing diverse values for the decay rate of CTLA-4, d_{CTLA 4}, above and below the default value d_{CTLA 4}=1, combined with different temporal duration of antigen attachment, τ _stim ( Figure 5 ).

The set of differential equations that defines the dynamical system is shown in the Supplementary Material 2 . The model equations are presented in Wolfram Notebook format in https://github.com/DrDavidMM/TCD4cell-activation-model-supplementary-material-.git.

For the computation of the differential equations system, Wolfram Mathematica 12.2.0.0 and open-source R studio have been used with the packages BoolNet, deSolve and ggplot2. For visual display of the interaction network, we use yEd graph editor 3.20.1 from yWorks. A link to the wolfram cloud public code has been added in the readme file in the GitHub repository.

We have assumed that T-cells initially are in a naive state, that is, TCR and CD28 are expressed and not activated. These conditions were implemented by considering that at time t = 0, TCR = 0, and CD28 = 0. MHC-antigen and CD80/86 act as inputs activating the TCR and CD28, respectively (although with the course of time CD28 may be displaced by CTLA-4). Optimal stimulation by MHC-antigen and CD80/86 was represented as A _{MHC / A} = 1, and A _{CD 8086} =1. Similarly, sub-optimal stimulation was represented by introducing values smaller than unity for either A _{MHC / A} or A _{CD 8086} (or both). As shown below, the level of stimulation defined by these parameters can produce activation or anergy. Since the metabolic profile of naive T-cells is characterized by a basal level of OXPHOS, modeling was performed assuming an initial value of OXPHOS = 0.2 and a high AMP/ATP ratio = 1, consistent with a significant activity of the nutrient sensor AMPK = 1.

Under optimal stimulating conditions, the network dynamics conduces to alternative states of sustained activation or CTLA-4-mediated arrest, each differing in their metabolic profile ( Figure 4 ). The presence of MHC-antigen and CD80/86 (at time t = 0) induces the activity of TCR and CD28 ( Figure 4A ). These interactions were maintained at a maximal level at stimulation times τ _stim = τ _{MHC / A} = τ _{CD 8086} = 15 units, until disengagement was induced at longer times t > τ _stim . During the stimulation time, dimerized CTLA-4 is transiently expressed at a low level, which however, is not sufficient to compete with CD28 for binding to CD80/86. Under these conditions, activation leads to the expression of the AP1, NFAT, and NFκB transcription factors at later times ( Figure 4B ) and, accordingly, transcription of the IL-2 gene (IL2G) and MTORC1 are also fully expressed ( Figure 4C ).

The predicted dynamic behavior of the activation process is coherent with trends inferred from the network interactive relationships depicted in Figures 2 and 3 . As indicated above, we assumed an initial naive state characterized by a low basal level of OXPHOS, a large AMP/ATP ratio, and activity of AMPK ( Figure 4D ). Upon T cell stimulation by TCR and co-stimulatory molecules, AMPK is further activated by Ca²⁺, LKB1 and AKT, boosting an early increase of OXPHOS. After reaching a peak of activity, OXPHOS is undermined due to the decrease of the AMP concentration and the afterward contribution of mTORC1, which impairs the repressive activity of AMPK and leads metabolism polarization towards glycolysis. Glycolysis contributes to the synthesis of cell-growth metabolites and maintenance of the ATP pool. This predicted behavior is congruent with experimental data showing that AMPK is activated early after T cell stimulation, which indicates that the engagement of mitochondrial metabolism is important for exiting quiescence (68), whereas the expression of enzymes pertaining to the glycolytic pathway is dispensable at earlier times (69). Thus, the model is in agreement with the proposal that AMPK activation ensures sufficient ATP availability to progress through full activation (53, 70). It can be observed that, even if TCR and CD28 were stimulated during a limited time span, the production of the AP-1, NFAT and NFκB transcription factors is held longer.

After TCR and CD28 stimulation, in due course the expression level of dimerized CTLA-4 increases, displacing CD28 from co-stimulatory molecules and leading to a state of activation arrest, or checkpoint. Since arrested cells have decreased levels of protein synthesis and expansion, the action of CTLA-4 may also have implications on the regulation of metabolism (71). To model this process, we assumed initial optimal engagement of TCR with the MHC-antigen complex, and CD28 with CD80/86, during an antigen presentation time τ _MHC/A = τ _{CD 8086} = 15 units, combined with a high-level activity of CTLA-4, induced by a low decay rate d_{CTLA 4} = 0.5.

Similar to the previous case, TCR and CD28 are fully activated due antigenic stimulation ( Figure 4E ). However, with the course of time the inhibitory action of CTLA-4 and the anergy factor NDRG1 are expresses, coinciding with the diminution of the levels of activity of TCR and CD28. This leads to an only transitory expression of the AP-1, NFAT, NFκB transcription factors ( Figure 4F ), the IL-2 gene, and mTORC1 activity ( Figure 4G ). The initial metabolic activity shows a pattern similar to that displayed in sustained activation. Nevertheless, the regulatory action of CTLA-4 eventually induces the decay of glycolysis ( Figure 4H ).

The model predicts that cell activation depends on whether the extent of the stimulation time of TCR and CD28, τ _stim , is long enough to overcome the regulatory activity of CTLA-4, which persists during a characteristic time, defined by τ _{CTLA 4} = 1/d_{CTLA 4}. The functionality of CTLA-4 resides in its continuous turnover, cellular location, and membrane delivery (18). In combination with activation-promoting factors, in the present model the decay rate of functional CTLA-4 is able to encompass these processes. In the simulation shown in Figure 4 , sustained activation ensued by assuming that τ _{MHC / A} = τ _{CD 8086} = 15, and d_{CTLA 4} = 5, whereas regulated activation was associated to τ _{MHC / A} = τ_{CD 8086} = 15, and d_{CTLA 4} = 0.5.

An analysis of the values of d_{CTLA 4} leading to states of sustained or regulated activation as a function of the stimulation time, τ _stim , is presented in Figure 5 . We observe that no activation arises for τ _stim < 7 units, while sustained activation ensues for 7 ≤ τ _stim ≤ 10 units, and regulated activation occurs for τ _stim > 10 units. The stage of no activation is associated with insufficient priming time by TCR and CD28. Sustained activation, independent of d_{CTLA 4}, is due to the fact that signaling from TCR and CD28 propagates downstream during a certain time without activating CTLA-4. On the other hand, regulated activation is driven by two conditions: first, that τ _stim is long enough, and second, that CTLA-4 activity persists longer than those of activation-inducing elements ( τ _{CTLA 4} = 1/d_{CTLA 4} > 1). In other words, CTLA-4 should be functional long enough to overcome a threshold level and its activity should be maintained to perform inhibition. Therefore, regulated activation persists for even longer values of τ _stim during which d_{CTLA 4} displays a quasi-periodic behavior with slight variations centered at d_{CTLA 4} ~ 0.5. The oscillation of the regulation threshold can be explained by considering that the expression of CTLA-4 depends on activation of the TCR and CD28. Once expressed, CTLA-4 initiates inhibition of signaling, which reduces its own expression and therefore allows activation again; this induces a stimulation-inhibition cycle which is downstream-propagated throughout the network. As a consequence, the threshold value of d_{CTLA 4} leading to regulation will depend on the expression level attained by activation inducers at a given phase of the cycle.

Variable levels of stimulation may originate from the progressive engagement of TCR to MHC-antigen complexes, the lack of a minimal half-life of interaction, diverse levels of antigen concentration, mutations in both cell receptors, etc. (63, 72, 73). Incomplete stimulation leads T cells to anergic states (14, 15, 74) and has been associated with mechanisms of tolerance “(adaptive tolerance”) or antigen presentation failure. On the other hand, defective CD28 co-stimulation allows the expression of the N-Myc Downstream Regulated 1 (NDRG1) protein, leading to anergy (75–77). The effect of different levels of stimulation on the downstream expression of the network components was determined. To simulate these phenomena, we considered initial conditions in which either TCR or CD28 were subjected to a “strong” or a “weak” stimulation. A first case considers A_{MHC / A} = 1 and A_{CD 8086} =0.5, and a second one, A _{MHC / A} = 0.5 and A _{CD 8086} = 1.

The left-hand side of Figure 6 shows the activation dynamics arising from a strong signaling from TCR and a weak co-stimulation through CD28. TCR and CD28 show full and partial activation, respectively. Low levels of AP-1 and interleukin 2 are transiently expressed. NDRG1 is expressed at a later time, coinciding with the activity decay of TCR and CD28 ( Figure 6A ). As before, this conduces to a full but transient expression of NFAT and NFκB, and a very small level of AP-1 ( Figure 6B ). A similar behavior if shown by mTORC1, and IL-2 is strongly suppressed ( Figure 6C ). Notably, the metabolic pattern is identical to that obtained in the case of regulated activation although in this case activity of CTLA-4 is not induced ( Figure 4E and Figure 6D ). The former results are consistent with reports indicating that TCR binding in the absence of CD28 ligation results in either apoptosis or a state of anergy that does not involve CTLA-4. Such anergic T cells are unable to produce IL-2 or proliferate on subsequent stimulation, even in the presence of co-stimulation (74, 78).

In the right-hand side of Figure 6 we show the alternative dynamics arising from a weak signaling of TCR, but a strong co-stimulation by CD28. In this case, TCR and CD28 are transiently expressed at respective low and high levels, but they subsequently decay independently of the action of CTLA-4 or the anergy factor NDRG1 ( Figure 6E ). On the other hand, no transcription factors are induced ( Figure 6F ) neither IL-2, although mTORC1 is transiently expressed ( Figure 6G ). Here, the metabolism profile ( Figure 6H ) also coincides with that obtained for regulated activation ( Figure 4H ).

It is widely documented that the different effector phenotypes such as Th1, Th2, Th17, and Treg present different energy requirements (56, 79, 80), and studies aimed to elucidate the intricate links between lymphocyte activation and metabolic reprogramming are needed (27). We simulated the individual conditions required for the differentiation of naive T cells into Th1, Th2, Th17 and Treg phenotypes, determined by the presence of the appropriate external cytokines.

Differentiation toward the effector Th cell lineages Th1, Th2, and Th17 is known to be reliant on mTOR activity, while inhibition of mTOR with rapamycin has been shown to favor Treg cell differentiation (56, 80, 81). mTOR complex1 (mTORC1) is formed with the scaffolding protein regulatory associated protein of mTOR (RAPTOR), while mTOR complex 2 (mTORC2) uses Rapamycin-insensitive companion of mammalian target of rapamycin (RICTOR) as a scaffold. All effector lineages, including Th2 cells, require mTORC1 activation (82, 83). Treg cells are a particular case in which the cell uses glycolysis to grow size and replicate; however, at stable stages they predominantly express Foxp3, which is a direct activator of AMPK and OXPHOS activity. This establishes metabolism as a key factor for the correct function of Treg cells (56).

Modeling of differentiation shows that, after activation, production of the characteristic cytokines starts first for Th1 and Th17 ( Figures 7A, C ), while the Th2 cytokines are induced at later times ( Figure 7B ). As shown before, OXPHOS is transiently upregulated upon activation ( Figures 4D and 7D ). Next, the three cell phenotypes develops a predominant and stable glycolytic profile. Instead, Treg differentiation shows the production of IL-10 and TGFβ at later times compared to Th1 and Th17 cytokines ( Figure 7E ). Interestingly, the simulation shows that glycolytic activity is followed by a strong and sustained polarization to OXPHOS, corresponding with AMPK activity ( Figure 7F ). This metabolic behavior could allow Treg cells to increase the pool of amino acids required for the synthesis of structural proteins, cell growth for cell clonal expansion, and the production of diverse metabolites necessary to carry out function. However, as the expression of Foxp3 stimulates AMPK, it induces again a polarization to OXPHOS leading cells to a stable regulatory stage. On the other hand, the model is in agreement with the observation that the functional form of CTLA-4 (CTLA-4dim) is upregulated and is constitutively expressed on Treg cells (84) ( Figure 7E ). Thus, the model effectively integrates the pathways that lead to the adjustment of the metabolic profile of different effector phenotypes.