### Statistical ensemble of gene regulatory networks

In statistical mechanics, a statistical ensemble refers to a large number (ideally infinite) of virtual copies of a system, each of which represents a possible and realistic state of the actual system. An ensemble therefore represents a probability distribution for the state of the system [39].

We generated copies of the above-constructed gene regulatory network of M1/M2 differentiation (i.e., the “system”) to build a statistical ensemble in order to calculate the probability distribution of the “state” of the differentiating macrophage, and to describe how this state changes with respect to biological stimuli such as bacterial infections or others.

Since a macrophage taken in isolation from the other immune cells and antigenic stimuli would not undergo any differentiation, we also incorporated in the model other cellular and molecular components so to have the main signals needed to represent both possible M1 and M2 differentiation pathways. This inclusion fully implements the mesoscopic (cell-cell interactions) of the model. The intra-cellular microscopic gene regulation level and the inter-cellular mesoscopic level, together constitute a multiscale system and enacts a full fledged adaptive immune response.

It is worth to note that for simplicity the model considers antigen-specific clonotypes of lymphocytes, i.e., it does not implements a full clonal selection process with antigen recognition.

In more details, we included in the ensemble virtual copies of B-cells (indicated *B*), macrophages (*M*), T helper lymphocytes (*T*), a generic antigen (*Ag*) bearing LPS as membrane molecules and cytokines IFNg (*IFNg*), IL-10 (*IL*10) and IL-4 (*IL*4).

Macrophages are subdivided in the undifferentiated phenotype (*M*
_{0}), in the terminally differentiated type-1 (indicated *M*
_{1}) and type-2 (*M*
_{2}).

Helper T lymphocytes are further divided in type-1 (indicated *H*
_{1}), type-2 (*H*
_{2}) and regulatory cells (*H*
_{
r
}) each of which secretes a different differentiation stimuli upon contact with the macrophage. Moreover, since cytokines are secreted by cells in particular activation states, we refined the representation of macrophages and lymphocytes including this information. Hence we defined macrophages in the resting (*M*
^{r}), active (*M*
^{a}) and presenting state (*M*
^{p}); helper lymphocytes (all three classes 1,2,*r*) in the resting (\(H_{1,2,r}^{r}\)) and active (\(H_{1,2,r}^{a}\)) state; B lymphocytes in the active (*B*
^{a}) and presenting state (*B*
^{p}). By “presenting” we mean that the antigen presenting cell B or macrophage has captured, internalised, processed and presented the antigen on the cell membrane together with the Major Histocompatibility Complex molecule. This process is not further detailed to keep the model simple. Table 2 contains a symbol legend.

The so constructed mesoscopic cell interaction level can be seen as a *bioreactor* in which each individual cell undergoes state transformations according to predefined immunologically-motivated rules [40]. Cells populate a 3D lattice where they freely diffuse and interact according to stochastic rules of the form

$$\alpha A + \beta B + \cdots \overset{p}{\longrightarrow} \gamma G + \delta D + \cdots $$

meaning that the rule applied to *α* instances of A, *β* instances of B etc. that are located in the same lattice site will produce *γ* instances of G, *δ* instances of D etc. with probability *p*.

The full list of rules is the following:

$$\begin{array}{*{20}l} M^{r} + IFNg & \overset{p_{1}}{\longrightarrow} M^{a} \end{array} $$

(1)

$$\begin{array}{*{20}l} M^{a} + Ag & \overset{p_{2}}{\longrightarrow} M^{p} \end{array} $$

(2)

$$\begin{array}{*{20}l} M^{p} & \overset{p_{3}}{\longrightarrow} M^{a} \end{array} $$

(3)

$$\begin{array}{*{20}l} M^{a} & \overset{p_{4}}{\longrightarrow} M^{r} \end{array} $$

(4)

$$\begin{array}{*{20}l} M^{p} + {H_{1}^{a}} & \overset{p_{5}}{\longrightarrow} M^{p} + {H_{1}^{a}} + IFNg \end{array} $$

(5)

$$\begin{array}{*{20}l} M^{p} + {H_{2}^{a}} & \overset{p_{6}}{\longrightarrow} M^{p} + {H_{2}^{a}} + IL4 \end{array} $$

(6)

$$\begin{array}{*{20}l} {M_{2}^{p}} + {H_{r}^{a}} & \overset{p_{7}}{\longrightarrow} M^{p} + {H_{r}^{a}} + IL10 \end{array} $$

(7)

$$\begin{array}{*{20}l} B^{a} + Ag & \overset{p_{8}}{\longrightarrow} B^{p} \end{array} $$

(8)

$$\begin{array}{*{20}l} B^{p} & \overset{p_{9}}{\longrightarrow} B^{a} \end{array} $$

(9)

$$\begin{array}{*{20}l} B^{p} + {H_{1}^{a}} & \overset{p_{1}0}{\longrightarrow} 2 B^{p} + 2 {H_{1}^{a}} + IFNg + Ab \end{array} $$

(10)

$$\begin{array}{*{20}l} B^{p} + {H_{2}^{a}} & \overset{p_{11}}{\longrightarrow} 2 B^{p} + 2 {H_{1}^{a}} + IL4 + Ab \end{array} $$

(11)

$$\begin{array}{*{20}l} H_{1,2,r}^{r} + M^{p} & \overset{p_{12}}{\longrightarrow} H_{1,2,r}^{a} + M^{p} \end{array} $$

(12)

$$\begin{array}{*{20}l} H_{1,2,r}^{a} & \overset{p_{13}}{\longrightarrow} H_{1,2,r}^{r} \end{array} $$

(13)

$$\begin{array}{*{20}l} {H_{1}^{a}} + B^{p} & \overset{p_{14}}{\longrightarrow} {H_{1}^{a}}+ B^{p} + IFNg \end{array} $$

(14)

$$\begin{array}{*{20}l} Ab + Ag & \overset{p_{15}}{\longrightarrow} \emptyset \end{array} $$

(15)

The first “reaction”, or rule R1 of Eq. (1), stands for macrophage activation by IFNg.

R2 represents antigen phagocytation, digestion and presentation by macrophages.

R3 stands for macrophages stopping presenting the antigen on cell surface.

R4 stands for macrophages returning to the resting state.

R5 indicates macrophage presenting the antigen to helper T cells of type-1 and inducing release of IFNg.

R6 is the same as before but for type-2 the released cytokine is IL-4.

R7 likewise, upon contact with *M*
_{2} cells, Treg lymphocytes release IL-10.

R8 is antigen phagocytation, digestion and presentation by B-cells.

R9 stands for B-cells stopping presenting the antigen on cell surface thus going back to the active state.

R10 is B-cells presenting the antigen to helper T cells of type-1 and inducing release of IFNg as well as antibodies (with the implicit assumption of B-cell differentiation to plasma B cells secreting antibodies).

R11 is the same as before but for type-2 the released cytokine is IL-4.

R12 antigen presenting macrophage activation of helper T-cells.

R13 stands for helper T-cells going back to resting.

R14 is releasing of IFNg by class-1 T-helper upon recognition of antigen peptides presented by B-cells.

R15 is the neutralisation of antigen by the antibodies.

These rules are executed in randomised order at each simulation step (to avoid biases). In principle *p* should be different for each rule, however, to keep the model simple we set ∀*i*=1,…,15,*p*
_{
i
}=1, that is, we allow all reactions to take place whenever all the necessary molecules are met on the same lattice point.

We can say that the ensemble consists of a Reactive Lattice Gas Automata (RLGA) on a three-dimensional cartesian lattice *L*×*L*×*L* with periodic boundary conditions. At the start, a large number of copies of cells populate the simulated volume. The relative proportion of lymphocytes and monocytes is set as in the standard human white blood cell counts [41]. Macrophages, that are the only agents whose internal dynamics is fully specified by the network show in Fig. 1, are initialised with the gene-state equal to zero meaning that the nodes of the network are in the inactive form (see below)). The rules R1–R15 (respectively Eq. (11)–(15)) are then applied in each lattice point. After that, all entities diffuse randomly on the lattice with equal speed. In this respect the lattice represents a homogenous media and a lack of differences in the diffusion speed of the various cells and molecules does not affect the “logic” of the model, which is what we are more interested in studying.

### The rule for macrophage differentiation

The last rules of the mesoscopic model implements the differentiation of the macrophages from the undifferentiated *M*
_{0} phenotype to either *M*
_{1} or *M*
_{2}. They realise the connection between the two levels of the multiscale description of the model.

As already mentioned, the networks in the ensemble are initialised as \(\bar {x_{0}} = \bar {0}\), that is, all genes are not expressed in the initially populating phenotype *M*
_{0}. The differentiation is then performed according to the following rule

$$\begin{array}{*{20}l} M_{0}^{a,p} & \overset{f_{1}}{\longrightarrow} M_{1}^{a,p} \end{array} $$

(16)

$$\begin{array}{*{20}l} M_{0}^{a,p} & \overset{f_{2}}{\longrightarrow} M_{2}^{a,p} \end{array} $$

(17)

where the macrophage has to be either active (superscript *a*) or presenting (superscript *p*) and the *f*
_{1} and *f*
_{2} are two probability functions defined as follows

$${} f_{1} \!\equiv\! \left\{ \begin{array}{lr} 1, & S(\bar{x}_{t} \wedge I,k) \,=\, X_{1}\\ 0, & \text{otherwise}\\ \end{array} \right. \qquad f_{2} \equiv \left\{ \begin{array}{lr} 1, & S(\bar{x}_{t} \wedge I,k) = X_{2}\\ 0, & \text{otherwise}\\ \end{array} \right. $$

with

$$S(\bar{x}_{t} \wedge I,k) = \underbrace{F \circ F \circ \cdots \circ F}_{\mathrm{k}} (\bar{x}_{t} \wedge I). $$

Thus, if the network state attained after *k* synchronous application of the Boolean dynamics *F*(·) reaches an attractor *X*
_{1} or *X*
_{2}, then the macrophage is considered committed, that is, terminally differentiated to the type-1 or type-2. Otherwise the state remains \(\bar {x}_{t+k}\) and the macrophages is still undifferentiated (i.e., *M*
_{0}). Specifically, \(S(\bar {x}_{t} \wedge I,k)\) is the state of the network after *k* application of the function *F*(·) starting from \(\bar {x}_{t}\) and after the input nodes have been updated with *I*=(*IFNg,IL*10,*IL*4,*LPS*) where LPS here means Ag with LPS on its surface membrane. In other words \(\bar {x}_{t} \wedge I\) is the vector \(\bar {x}_{t}\) with the values corresponding to IFNg, IL10, LPS and IL4 updated with the values in *I* (see extracellular stimuli layer in Fig. 1). Now if

$$\forall k>0, k\in \mathbb{N}, \quad S(\bar{x}_{t} \wedge I,k) = \bar{x}_{t} \wedge I $$

this means that the dynamics has converged to a fixed point. We call *X*
_{1} the fixed point associated to the phenotype *M*
_{1} and *X*
_{2} the fixed point associated to the phenotype *M*
_{2}. By construction *X*
_{1}≠*X*
_{2} and therefore it follows that in rules R16 and R17 (respectively in Eqs. (16) and (17)) it must be that at any time *f*
_{1}+*f*
_{2}≤1 (i.e., none or only one of the two steady state is reached). In other words, each network of the ensemble representing a macrophage initially in the undifferentiated state has the chance to differentiate at each discrete step on the basis of its current state \(\bar {x}_{t}\) and the input *I*, which in turn is given by the local (i.e., in the same lattice site) “concentration” of IFNg, IL10, IL4 and antigen with LPS in the current lattice location (the lattice location has not been indicated in the formalism above for simplicity).