The mathematical methodology employed in this work is summarized in Figure 1. An overall assumption is that activation, differentiation and effector processes in cells comply with logical rules comprising a regulatory network of node interactions, underlying a (metaphorical) epigenetic landscape (20). Figure 2 displays the regulatory network for macrophage differentiation put forth in this work. It includes 128 components that incorporate current experimental evidence of the molecular interactions involved in the differentiation of the macrophage phenotypes M1, M2a, M2b, and M2c. The network consists of inputs (environmental stimulating agents) and nodes (representing receptors, enzymes, second messengers, transcription factors, genes, etc.) linked through activating or inhibitory interactions.

In a first level, the network interactions are described by Boolean rules that allow to set up the network architecture, including feedback and switching modules. In the Boolean scheme, logic rules are defined by dichotomous variables with truth values, 0 (unexpressed), or 1 (expressed), and the network dynamics is described in terms of discrete-time mappings whose steady states (attractors) define cell phenotypes. A Boolean approach is fundamental for building the set of regulatory interactions; however, this procedure may be insufficient to account for the complex phenomenology observed in specific experimental studies. A more suitable modelling of biological systems may be achieved by introducing a continuous logical analysis (21, 22). In this approach, the expression levels and concentrations of the signaling elements can acquire any value within a continuous range limited only by functionality constraints and concomitant continuous dynamics. Then, we extended our initial analysis to consider a regulatory network characterized by Fuzzy Logic, aimed to provide formal foundation to logical inferences involving uncertain or vague reasoning and has found applications in physical, control, biomedical, and linguistic sciences (23). In this approach, node interactions are described by continuous-valued logical propositions where the expression values of the variables, q_i, can be any real number ranging from completely false, 0, to totally true, 1 (13, 14, 24). The corresponding fuzzy interaction rules may be straightforwardly derived from their Boolean counterparts by replacing ordinary logic connectors by algebraic expressions, as shown in Figure 1. The kind of approach has been previously employed by some of the present authors to describe differentiation and plasticity processes in floral organs (13), pancreatic beta cells (25), and CD4 T-cells (26–28).

To proceed, we denote the expression levels of the network agent i at a time t by a continuous variable q_i(t), with 0 ≤ q_i ≤ 1. The network dynamics is then described by a set of ordinary differential equations (ODEs):

Here, the input μ[wi]=1/(1+exp[−β(wi(qi,…qn)−θth)]) is a sigmoid function representing the truth value of the fuzzy interaction rule w_i with respect to a threshold θ_th, and β is a saturation rate. The second term yields an interaction decay at a rate αi=1/τi, where τ_i is a characteristic expression time. In this work, we assume β = 10, θ_th = 1/2, and α_i = 1. The differentiation dynamics induced by any considered cell microenvironment can be simulated by solving the ODE system (1) for a fixed set of initial conditions {q_i(0)}. Solutions exhibit a long-time steady behavior, qiT=qi(t=T≫1), giving rise to a set of steady-state solutions {qiT}, as it will be shown in Results. In the simulations reported here, T = 30 is time long-enough long time long-enough. Alternatively, the set of steady-state solutions {qiss}, may be derived from the condition dq_i/dt = 0, yielding

and we may identify qiss=qiT. Accordingly, the steady-state expression level of any network’s agent is determined by the truth value of its steady logical connections, modulated by its decay rate α_i. In the case αi≫1, then qiss≈0.

Simulations of the regulatory network dynamics were performed using a Python code involving numerical solvers for the ODEs described above; specifically, the ‘odeint‘ function from the ‘scipy‘ library. This allowed the implementation of an interactive program Final monocyte network.ipynb that may be used as a tool in the Google Colab application to test the dynamical expression of the network components under different microenvironmental and initial conditions. In particular, the initial conditions defining an unstimulated monocyte, {q_k(0)}, were set to zero for most input variables (except for AMPK and GSK3β activities, which were set to 1). Then, different stimulation conditions were applied, including LPS to activate TLR4, ssRNA to activate TLR7, IgG immune complexes to engage Fcγ receptors, and combinations of these with cytokines (e.g. IL-4, IL-10). Each simulation was carried out over a predefined time interval, capturing the evolution of cell-produced cytokine levels, metabolic state, and transcription factor activities. The final steady states were analyzed to evaluate macrophage polarization markers under different stimulation conditions. Key variables, including cytokine expression, metabolic activity, and transcriptional factors dynamics, were plotted as time-series data to visualize the trajectories of monocyte activation and differentiation. Visualization was performed using the ‘matplotlib‘ library. Comparative analyses across conditions highlighted the influence of receptor activation and microenvironmental factors on the macrophage phenotype.

The results provided by the model have been validated against experimental data reported in the literature, ensuring consistency in cytokine profiles, metabolic shifts, and receptor-mediated differentiation. Sensitivity analyses were conducted to assess the robustness of the results to variations in parameter values. The Python code used for simulations is available upon request, ensuring reproducibility and transparency for future studies.

Robustness implies the capacity of a complex network to maintain its functionality subject to the action of topological or dynamical perturbations. Here, we assume that the network topology is already defined by selective forces, and we thus constrain our analysis to two different perturbation sources: i) the influence of randomly distributed initial values of the microenvironmental and cell variables on the differentiation stability. ii) the effect of noisy perturbations on the intrinsic differentiation dynamics, for fixed values of the initial microenvironment and cell variables.

The network, displayed in Figure 2, includes nodes and green or red arrows that represent activating or inhibitory node interactions. The specific exogenous inputs are listed in Table 1, while the resulting Boolean interaction rules are presented in Table 2. For a visual assessment of the built-in elements, the network is presented as interconnected modules, whose different colors denote specific functions: extracellular stimuli, TLR’s, activating IgG immune complexes, cytokine receptors, downstream signals, metabolism, transcription factors, and secreted cytokines (output cytokines). Details of particular interactions can be visualized using the Yed software as MONOCYTE NET.graphml in the https://grci.mx web site.

Inflammatory pathways are represented in Figure 3A, which displays the convergence of activating and inhibitory signals on the NFκb node, as well as the consequent production of inflammatory cytokines (TNFα, IL-1β, IL-12, IL-6, IL-18, and IL-33). This module includes the activity of metabolic elements such as mTOR and mTORC1, as well as transcription factors and cofactors such as CREB1, AP-1, IRF3 and IRF7 that are either activated by the TNFα and FCγ receptors, or by TLR engagement. Since the effects of Akt on macrophage polarization appear to be Akt isoform-specific (30), its isoforms Akt1, Akt2 and Akt3, have been incorporated in the network regulatory pathways (Figure 3B). In particular, experiments with macrophages deficient in Akt1 or Akt2 revealed that this default gives rise to M1 and M2 macrophages, respectively (31).

The pathways leading to the anti-inflammatory role of CREB1 are shown in Figure 3B, highlighting GSK3β as inhibitor of the CREB1 function. It can be observed that this module displays the downstream pathway:

which involves four consecutive inhibitory interactions, leading to the activation of the last downstream element; therefore, the presence of IFN-γ in the cell microenvironment promotes the activation of NF-κB. Conversely, in the absence of IFN-γ, an anti-inflammatory macrophage profile may arise, as may be inferred from the middle factors in the route: upon activation (by PI3K) the Akt1 and Akt3 mediators hinder the action of the constitutively expressed factor GSK3β (32), leading to anti-inflammatory effects of CREB1 through NF-κB inhibition.

Figure 3C shows a module for AMPK-controlled glycolysis and oxidative phosphorylation (OXPHOS), reflecting the metabolic profiles of M1 and M2 phenotypes. In this module, AMPK is a central energy sensor of the AMP/ATP ratio and involves a negative feedback loop with mTORC1, defining a switch driving either OXPHOS or glycolytic activity.

The signaling pathways induced by the exogenous immunoregulatory cytokines IL-4 and IL-10 are shown in Figure 4. Receptor signaling by these cytokines leads to inhibition of NFκb, promotes the production of IL-10 and the activation of AMPK, directing the macrophage metabolism toward OXPHOS.

We applied the mathematical model to the analysis of a set of relevant initial conditions introduced as inputs of the system (see Table 1). Figures 5-8 depict the time evolution of network components induced by diverse cell microenvironments, and Figure 9 summarizes the stable states obtained.

The network faithfully reproduces the spectrum of macrophage functional states, known as M1, M2a, M2b, and M2c. Figure 5 displays the dynamics of metabolic changes, transcription factor activity, cytokine production and markers associated with the M1, M2a, M2b and M2c phenotypes. It can be observed that co-stimulation with LPS and IFN-γ drives a full M1 polarization (Figure 5A). The system converges to a sustained activation of NF-κB, a marked increase in glycolytic flux, and robust production of all pro-inflammatory cytokines (IFN-α, IL-1β, IL-6, IL-12 and TNF-α). This outcome reproduces the metabolic and secretory profile typical of classically activated macrophages (3, 33). On the other hand, IL-4 alone (Figure 5B) promotes a M2a state: AMPK-driven oxidative phosphorylation prevalence over glycolysis, activity of PPAR-γ, JMJD3, and STAT6, and production of IL-10. These features mirror the anti-inflammatory and tissue-repair phenotype described for IL-4–stimulated macrophages (34, 35). When LPS is paired with IL-1β (Figure 5C), the model yields the M2b response: a transient burst of TNFα is induced (similarly to the case of LPS alone (Figure 6A) Gycolysis remains dominant, and a mixed set of pro- and anti-inflammatory cytokines is produced, with activation of Erk and STAT3 and production of IL-10. This phenotype recapitulates observations that LPS+ IL-1β stimuli can generate hybrid profiles (36, 37). Finally, IL-10 alone (Figure 5D) drives a full M2c program: sustained OXPHOS metabolism, strong GSK3β activation (not shown), and production of IL-10 through the activity of STAT3; pro-inflammatory cytokines are fully suppressed. This result aligns with the immunoregulatory profile induced by IL-10 in vitro and in vivo, and are consistent with tissue repair functions (3, 36, 37). Thus, the monocyte model reproduces the hallmark metabolic and cytokine signatures of canonical macrophage phenotypes.

Signaling events related to immunomodulation by IL-4 and IL-10 were incorporated in the network, as seen in Figure 4 and Table 2. IL-4 receptor signals for stimulation of the histone demethylase JMJD3, STAT6, and PPAR-γ, a known inhibitor of NF-κB. STAT6 activation leads to the expression of the IL-10 gene (34, 38–40). Likewise, STAT6 activates the Krüppel-like factor 4 (KLF4), which contributes to the inhibition of NF-κB and therefore of the TNF-α synthesis. Thus, IL-4 can inhibit NF-κB via PPAR-γ and KLF4 even in the presence of IFN-γ.

The addition of IL-4 to the LPS+IFN-γ stimulus promotes partial differentiation to the anti-inflammatory phenotypes M2a and M2b, due to the positive action of JMJD3 on the IRF-4 transcription factor and inhibition of STAT5 (Figure 6C). On the other hand, activation of STAT6 by IL-4 activates SOCS-1, an inhibitor of STAT1. In the model, inhibition of STAT5 and STAT1 reduces IL-12 production (Figures 6C), as reported (41, 42). M2 macrophages tend to depend more on mitochondrial oxidative phosphorylation and fatty acid oxidation than glycolysis (8, 35, 36, 43, 44). The preferential use of OXPHOS over glycolysis by IL-4 or IL-10-stimulated M2 macrophages is reproduced by the model (Figures 6C, D). Thus, the model reproduced de role of IL-4 and IL-10 in steering differentiation away from a strictly pro-inflammatory state.

The CREB1 inhibitor GSK3β is constitutively expressed in monocytes (19). Upon TLR4 stimulation, GSK3β is inhibited via the PI3k-Akt-dependent pathway, specifically through the Akt1 and Akt3 kinases, so that CREB1 remains active, inhibiting the activity of NF-κB and promoting the production of IL-10 (Figure 3B) (reviewed by (45). The constitutive expression of GSK3β was incorporated in the model by asignyinig it an initial 1 value. Simulations showed that TLR4 activation by LPS alone drives the activity of Akt1 and Akt3, inhibition of GSK3β and induction of CREB1 activity (Figure 7). As a result, an incomplete pro-inflammatory profile is obtained, marked by a transient activation of NF-κB, with a decay coincident with the induction of CREB1 activity. Accordingly, IL-1β and TNF-α are transiently produced. The production of IFN-α and IL-6, together with the activity of the transcription factors AP-1 and IRF7 and a glycolytic metabolism are observed.

The effect of variable levels of INF-γ on the M2b/M2c phenotype induced by LPS is shown in Figures 7A-D. INF-γ inhibits the protein kinases Akt1 and Akt3 activities, allowing the function of GSK3β, leading to CREB1 inhibition. In this condition, NF-κB is stably active and a complete M1 cytokine profile is induced, as shown in (Figure 7A). Thus, modeling results are compatible with a relevant role of CREB1 in the inhibition of the pro-inflammatory response via inhibition of NF-κB and promotion of IL-10 synthesis. The model predicts that a full activity of NF-κB is maintained even when the levels of INF-gamma decreases to half of the optimum value (Figures 7B-C). However, lower amounts of this cytokine (0.25) allow the activity of AKT1 and AKT3, with the concomitant depletion of GSK3β, which in turn allows the activity of CREB1 and the subsequent inhibition of NF-κB. As a result, the M1 polarization is disrupted, and mixed M2 phenotypes are generated. The absence of IFN-γ yields the same result (Figures 7D, E). Thus, the model shows that the stabilization of NF-κB activity by IFN-γ can be reached through inhibition of the pathways leading to activation of the CREB1 inhibitor GSK3β, and highlights that GSK3β is necessary for the TLR4-mediated M1 inflammatory response, as suggested before Xia et al. (32)Ko and Lee (46).

The binding of IgG-antigen immune complexes (IC) to Fc receptors induces strong phagocytosis, antigen presentation capabilities and production of pro-inflammatory cytokines in monocytes (47–49). Figure 8A shows the result of modelling the addition of IgG-antigen immune complexes to LPS+IFN-γ-activated monocytes. Besides the induction of all the pro-inflammatory cytokines, signaling through the Fcγ receptor induces the promotion of phagocytosis, antigen processing, and MHC class II presentation through the Syk-VAV-Rac pathway (47–49). As shown above, the presence of IFN-γ keeps CREB1 inactive and so NF-κB is active.

Figure 8B shows that the input of IL-4 to the LPS+IFN-γ+IgG (IC) condition, blocks the NF-κB function and TNF-α synthesis, while the activation of JMJD3 promotes the function of the IRF4 transcription factor. IRF4 inhibits STAT5, thus inhibiting IL-12 secretion; however, a transient production of this cytokine is still supported by STAT1, NF-κB, and AP-1 activities (Figure 8B). IL-1 is not produced, as described (50). The production of IL-10 is induced by IL-4 respect to the previous condition and the metabolism is shifted to oxidative phosphorylation. As a whole, IL-4 leads to an M2 profile. The results of the modelling agree with the notion that IL-4 can counteract the stimulation by LPS, IgG-immune complexes, and IFN-γ in monocytes (51–53), whereas the phagocytosis and antigen processing capabilities are maintained.

Intracellular TLR7 and TLR8 recognize ssRNA and signal though the TRAF3 and MyD88 elements. Whereas signals downstream MyD88 induce the production of pro-inflammatory cytokines, TRAF3 induces the activity of interferon regulatory factors 3 and 7 (IRF3 and IRF7), which promote the production of type I interferons (IFNα and IFNβ), essential components of the anti-viral response (54). These rules were included in the model (Table 2) and the effect of ssRNA stimulation was analyzed. Results showed that ssRNA-only stimulation induces a glycolytic metabolism and only a transient expression of NF-κB and IRF-3 (see the interactive program Final monocyte network.ipynb in the Google Colab application). IFN-γ greatly potentiate the M1 pro-inflammatory profile (55–57), with a stable induction of NF-κB and IRF3. IgG (IC)) added phagocytosis and antigen presentation activities to the ssRNA+IFN-γ situation (Figure 8C). The expression of M2 markers was not obtained. Therefore, the model integrates the view that IFN-γ supports the production of pro-inflammatory cytokines, whereas IgG immune complexes induce phagocytosis and antigen presentation capabilities during the ssRNA+IFN-γ response. As in the case shown in Figure 8B, the addition of IL-4 produces the inhibition of NF-κB, a transient production of IL-12, and the expression of M2 markers. In addition, the expression of IRF3 is inhibited (Figure 8D).

Stochastic models have been instrumental in exploring system stability in the face of fluctuations arising from random variations of the exogenous cell microenvironment or inner signaling pathways. A robustness index, R, was calculated for perturbations of the initial levels of phenotype-inducing microenvironments: M1 (LPS and IFNG-γ), M2a (IL-4), M2b (LPS and IL-1β) and M2c (IL-10). It may be observed that up to perturbation amplitudes ϵ_R∼ 20%, the phenotype differentiation processes show a strong robustness with respect to initial-condition variations. In the case with ϵ_R∼ 30%, only the robustness of the M2b phenotype is significantly reduced (Table 3). A similar analysis performed for the β parameter indicated that up to ϵ_R∼ 30%, the robustness index R ≥ 0.98 in all the above mentioned cases. As an alternative stability proof, we translated the deterministic model into a stochastic scheme (Figure 1 and 2.3.2) to explore the effect of intrinsic noise on the system dynamics (13, 28). Different levels of noise were introduced into the equations as described in Section 2.3.2, and 1,000 iterations were performed to obtain an average percentage of differentiated cells under each tested condition. Figure 10 shows differentiation efficiencies (ϵ_D) under M1, M2a, M2b and M2c input conditions, considering 0% to 50% noise. It can be seen that the differentiation process of all phenotypes is very robust, with ϵ_D ≥ 0.9 (90%), except for M2c cells for noise levels higher than 30 %, where ϵ_D exhibits a drastic reduction.

The model captures the dynamics of macrophage functional polarization in agreement with experimental observations, providing a solid framework for the comprehensive description of monocyte behavior. Further work should include interactions improving the prediction of functional features, like those driven by metabolic changes mediated by mTORC2, cytoskeletal dynamics and chemokine production. Other relevant functions of GSK3β in addition to inhibition of CREB1 (46) may be considered. The stimulatory effect of NF-κB on the expression of phagocytosis receptors different from FcγR, like Dectin and the mannose receptor are incorporated in the model; mathematical simulation of signaling from them is a work in progress.

The modular structure of the monocyte/macrophage network allows its integration with regulatory networks pertaining to other cells of the immune system, or to tissues interacting with them through soluble factors. Work in progress involves connecting the present macrophage network with a CD4 T cell lymphocyte differentiation network (26) with the aim of simulating the immune response to respiratory infections (64), and that of diseases with chronic inflammatory origin, such as type 2 diabetes (25).