The role of T-helper/T-suppressor ratio in the adaptive immune response: a dynamical model

TL;DR

This paper presents a dynamical model of the adaptive immune response demonstrating that the T-helper/T-suppressor ratio critically influences immune responsiveness, supporting its use as an immunosuppression index.

Contribution

The study introduces a robust dynamical systems model showing the T-helper/T-suppressor ratio's pivotal role in immune response transitions, aligning with experimental observations.

Findings

The model is robust against noise and parameter changes.

A critical T-helper/T-suppressor ratio triggers a transition to immunosuppression.

The ratio's impact aligns with experimental data on immune response.

Abstract

Recent experimental studies have suggested the ratio between T-helper and T-suppressor lymphocytes as an index of immunosuppression in HIV, cancer, immunosenescence and inflammatory and auto-immune diseases. However, a quantitative understanding of the impact of this ratio on the immune response has lagged behind data and its validity as a tool for prognostic monitoring or therapeutic target remains an open question. In this work, we use statistical physics and dynamical systems approaches to analyze the time-dependent response to an antigen, of a simplified model of the adaptive immune system, which comprises B, T-helper and T-suppressor lymphocytes. The model is remarkably robust against changes in the noise level and kinetic parameters, but it is very sensitive to changes in the ratio between T-helper and T-suppressor lymphocytes, exhibiting, in particular, a transition from a responsive to an immuno-suppressed phase, as the ratio is lowered below a critical value, which is in line with experiments. This result supports the validity of the T-helper/T-suppressor ratio as an index of immunosuppression and may provide a useful theoretical benchmark to interpret and compare experiments.

Equations

$$Ag{\to }^{r}2Ag$$

$$Ag+B{\to }^{{\pi }^{+}}APB$$

$$T+APB{\to }^W{T}^{\star }+APB$$

$$APB{\to }^{{\pi }^{-}}B$$

$$T^{\star }+B{\to }^{{\lambda }^{+}}{T}^{\star }+2B$$

$$T^{\star }+B{\to }^{{\lambda }^{-}}{T}^{\star }$$

$$B+B{\to }^{\delta }B.$$

$$p({\xi }_i^{\mu })=\frac{c_{\mu }}{N^{\gamma }}{\delta }_{{\xi }_i^{\mu },1}+\left(1-\frac{c_{\mu }}{N^{\gamma }}\right){\delta }_{{\xi }_i^{\mu },0}$$

$$\frac{d}{dt}{p}_{\mu }$$

$$\frac{d}{dt}{\psi }_{\mu }$$

$$P(\eta )=\frac{1+ϵ}{2}{\delta }_{\eta ,1}+\frac{1-ϵ}{2}{\delta }_{\eta ,-1}$$

$$R=\frac{1+ϵ}{1-ϵ}.$$

$${\sigma }_i(t+\Delta )=\theta (\sum _{\mu =1}^P{\xi }_i^{\mu }{p}_{\mu }(t)-z(t))$$

$$\frac{d{b}_{\mu }}{dt}$$

$$\frac{d}{dt}{b}_{\mu }$$

$$m_{\mu}({\mbox{\boldmath$\sigma$}})=\frac{1}{c_{\mu}N^{1-\gamma}}\sum_{i=1}^{N}\sigma_{i}\eta_{i}\xi_{i}^{\mu}$$

$$p_{t+\Delta }({\sigma }_i)=\mathcal{P}\left(z\le (2{\sigma }_i-1)\sum _{\mu }{\xi }_i^{\mu }{p}_{\mu }(t)\right)$$

$$\mathcal{P}(z\le x)=\frac{1}{2}\left(1+\tanh \frac{\beta x}{2}\right),$$

$$W_t({\sigma }_i)=\frac{1}{2}\left[1+(1-2{\sigma }_i)\tanh \frac{\beta }{2}\sum _{\mu }{\xi }_i^{\mu }{p}_{\mu }(t)\right]$$

$$\frac{d}{dt}p_{t}({\mbox{\boldmath$\sigma$}})=\sum_{i}[p_{t}(F_{i}{\mbox{\boldmath$\sigma$}})W_{t}(1-\sigma_{i})-p_{t}({\mbox{\boldmath$\sigma$}})W_{t}(\sigma_{i})]$$

$$\frac{d}{dt}⟨{\sigma }_j⟩=⟨(1-2{\sigma }_j)W_t({\sigma }_j)⟩$$

$$\frac{dm_{\mu}}{dt}=\frac{N^{\gamma}}{2c_{\mu}}{\langle}\xi^{\mu}\eta[1+\tanh\frac{\beta}{2}\sum_{\nu}\xi^{\nu}p_{\nu}]{\rangle}_{\eta,{{\mbox{\boldmath$\xi$}}}}-m_{\mu}$$

$$P(\eta,{\mbox{\boldmath$\xi$}})=\frac{1}{N}\sum_{j=1}^{N}\delta_{\eta,\eta_{j}}\delta_{{\mbox{\boldmath$\xi$}},{\mbox{\boldmath$\xi$}}_{j}}$$

$$\frac{dm_{\mu}}{dt}=\frac{\epsilon N^{\gamma}}{2c_{\mu}}{\langle}\xi^{\mu}[1+\tanh\frac{\beta}{2}\sum_{\nu}\xi^{\nu}p_{\nu}]{\rangle}_{{\mbox{\boldmath$\xi$}}}-m_{\mu}.$$

$$\frac{dm_{\mu}}{dt}\!=\!\frac{\epsilon}{2}\!\left[1\!+\!{\langle}\tanh\frac{\beta}{2}\Big{(}p_{\mu}\!+\!\sum_{\nu\neq\mu}^{P}\xi^{\nu}p_{\nu}\!\Big{)}{\rangle}_{\{\xi^{\nu}\}}\!\right]-m_{\mu},$$

$$\frac{d{b}_{\mu }}{dt}$$

$$p_{\mu }=a_{\mu }{b}_{\mu }{\psi }_{\mu }$$

$$\frac{d{m}_{\mu }}{dt}$$

$$({\psi }_{\mu },p_{\mu },m_{\mu },b_{\mu })=\left(0,0,\frac{ϵ}{2},\max \left\{0,\frac{{\kappa }_{\mu }ϵ}{2}\right\}\right).$$

$$\frac{d{m}_{\nu }}{dt}$$

$$\frac{d{p}_{\nu }}{dt}$$

Full text

1 Introduction

Recent years have seen a surge of mathematical and computational approaches aimed at modelling the immune system from a systemic perspective. Many popular models have been formulated in terms of ordinary differential equations [1, 2, 3, 4, 5, 6, 7, 8] and partial differential equations [9, 10, 11], as these provide intuitive frameworks to understand the dynamics of average quantities (e.g. cell concentrations). However, they normally ignore microscopic details and stochasticity, due to noise in the biological environment or fluctuations in cellular densities, and generally require the estimation of a large number of unknown parameters. Agent-based simulations [12, 13, 14, 15] and machine learning approaches [16, 17, 18, 19] have been successful in incorporating statistical noise and microscopic information (e.g. cellular interactions, antibodies sequences etc.), however, they usually require more significant computational efforts. In recent years, statistical mechanical approaches [20, 21, 22, 23, 24, 25] have supplied useful models to handle the immune system’s intricate patterns of interactions while remaining analytically tractable. However, most of these studies are performed at equilibrium: statistical mechanical models looking at the dynamics of the immune system after antigenic stimulus are still lacking or are at their infancy [26, 27]. One of the aims of this work is to advance non-equilibrium statistical mechanical models of the adaptive immune system.

The adaptive immune system is a complex network of cells that work together to defend the body against pathogens such as bacteria, virus or tumor cells. The major agents of adaptive immunity are lymphocytes, which are broadly divided into T and B cells. Each B cell is equipped with a large number $n={\mathcal{O}}{(10^{5})}$ of identical B cell receptors (BCR) on its surface, which are able to bind to specific pathogens. When a B cell comes across a pathogen that can bind its receptors, it ingests it, partially degrades it, and exports fragments of it, i.e. antigens, to the cell surface, where they are presented in association with proteins known as MHC molecules. T cells are also equipped with receptors (TCR), but in contrast to B cells, they do not bind to the antigen directly, but to the MHC-antigen complex on the surface of antigen presenting cells, such as B cells, dendritic cells and other immune cells. T cells can be divided into CD4+ (helper) and CD8+ (cytotoxic/suppressor) cells, which bind different classes of MHC-molecules. CD4+ cells bind to class II MHC-molecules and have the role of signaling B cells to initiate an immune response. CD8+ cells bind to class I MHC-molecules and can be divided into cytotoxic T cells, which destroy infected or cancer cells, and suppressor T cells, which switch off the immune response. In healthy hosts, the ratio between CD4+ and CD8+ cells is normally above one, while low or inverted CD4+/CD8+ ratios have been associated with impaired immune function in inflammation and autoimmune diseases [28, 29, 30, 31], cancer [32, 33] and immunosenescence [34, 35, 36]. Furthermore, recent studies have suggested that the CD4+/CD8+ ratio may affect the progression of HIV infections and response to antiretroviral therapies [37, 38, 39, 40, 41], and may be a marker for viral reservoirs [42, 43, 44, 45] in HIV-positive patients. However, a sturdy relation between CD4+/CD8+ ratios and viral reservoirs has not been proven and the general validity of the CD4+/CD8+ ratio as a tool for prognostic monitoring or therapeutic target remains an open question [46].

With the advance of molecular immunology, the existence of several subtypes of CD4+ and CD8+ lymphocytes has been documented and a subtype of CD4+ cells functioning as suppressor cells has been uncovered [47, 48]. These were named T-regulatory cells, and their biomarkers have been identified recently after a long-standing debate [49] and localized in the cell [50]. Similarly to the CD4+/CD8+ ratio, the T-helper/T-regulatory cell imbalance has been suggested as an index of immunosuppression in cancer patients [33, 51, 52], HIV patients [53, 54, 55], immunosenescence [56] and in inflammation diseases [57]. Results generally suggest that lower values of T-helper relative to T-regulatory cells are associated with unfavourable prognosis. However, different studies have used different biomarkers to identify T-regulatory cells (resulting in populations with different degrees of enrichment in T-regulatory cells) and have focused on different populations of T-helper cells (sometimes all CD4+ cells, other times only subtypes of them). In addition, cell imbalance has been measured in different ways (sometimes in terms of differences, others in terms of ratios between cell concentrations), and is subject to fluctuations in T cell counts, meaning that a global consensus on its reliability as a biomarker for immunosuppression has not been reached.

To add complexity, the CD8+/T-regulatory ratio has also been associated with clinical outcome in cancer recently [58], and with the increasing number of cell subtypes being uncovered, many different combinations of cell subtypes can be considered. However, assessing the significance of any observed correlation between the imbalance of different combinations of cell subtypes and prognosis in a particular disease, based on data from a limited number of patients, cannot be done reliably without a suitable theoretical model.

The aim of this work is to provide a theoretical framework to understand the impact of the T-helper/T-suppressor ratio on the response of the immune system to an antigen. We will move away from a detailed characterization of T-helper and T-suppressor cells in terms of biomarkers, and we simply assume that there are two broad categories of T cells: T-helper, that activate B cells, and T-suppressor, that inhibit B cells, in more or less direct ways, including cytotoxic and regulatory activity. Although our definition of T-suppressor includes different cell lineages i.e. CD8+ cytotoxic, CD8+ suppressor and CD4+ regulatory, we note that CD4+ regulatory cells are present in normal peripheral blood in low numbers, typically 5-10% of CD4+ T cells [59], so the T-helper/T-suppressor ratio considered in this work, does not deviate significantly from the CD4+/CD8+ ratio considered in the literature. This paper will articulate as follows: In Sec. 2 we provide a brief description of the adaptive immune system and summarize the main cellular reactions involved in its response to an antigen. In Sec. 3 we define a network model of the adaptive immune system, comprising T-helper, T-suppressor and B cells, and study its response to a single antigen, for a large range of its control parameters, in the thermodynamic limit. In Sec. 4 we simulate the reactions between immune cells using a Gillespie algorithm and check the accuracy of our model. In Sec. 5 we look at possible generalizations and extensions of our model. Finally, in Sec. 6 we summarize our results and propose pathways for future work.

2 The adaptive immune system: a brief description

T-helper cells get activated upon binding antigen presenting cells, like B cells, dendritic cells or other immune cells presenting antigens on their surface. When active, T-helper cells proliferate and release proteins called cytokines, which activate B cells. Active B cells undergo clonal expansion, i.e. form many copies of identical cells sharing the same antigen receptors, and secrete antibodies, i.e. a free form of those receptors, that can recognize and neutralize the antigen [60]. Clonal expansion usually involves migration of B cells to the germinal centre111Germinal centres are transient structures that form within peripheral lymphoid organs in response to T cell-dependent antigen, where B cells proliferate at a rate that is unparalleled in mammalian tissues [61]. Random mutations during clonal expansion cause the production of antibodies with a broad range of binding affinities for their antigen. B cells with unfavourable mutations will not get sufficiently activated by the antigen and T-helper cells and will die, while those with improved affinity will be stimulated to clone themselves. This leads to an effective selection processes within the germinal centre, known as affinity maturation, which results in the production of antigen-specific B cells within one week, in mammalians [62].

T-suppressor cells switch off the immune response. Although the mechanics of regulatory T cells is less clear than that of cytotoxic T cells [63, 64], their ability to directly suppress B cells has been well-documented in recent studies [65, 66, 67, 68, 69, 70]. From a modelling perspective, we will not distinguish between cytotoxic and regulatory cells, and we will simply regard them as T-suppressor cells, with the ability to reduce the number of effector 222Effectors B cells are B cells that secrete large volumes of antibodies. B cells, by either killing them or inhibiting their production of antibodies.

Below, we list the reactions that we shall include in our model. These occur in different parts of the organism, e.g. T and (antigen presenting) B cells mostly interact in the lymph nodes (organs of the lymphatic system), while lymphocytes and antigens also interact in infected tissues and other parts of the lymphatic system. In this study, we will not consider spatial effects and we will regard all the reactions below as taking place in the broad lymphatic system, which communicates lymphatic organs (such as lymph nodes) with infected tissues and with the vascular system. For simplicity, we suppose that each pathogen has a single epitope, that we will loosely call antigen. We have:

• •

Pathogens replicate in the host at rate $r$, leading to proliferation of antigens

[TABLE]

• •

B cells bind to antigens, thus becoming antigen presenting B cells (APB), at a rate $\pi^{+}$, which depends on the affinity between their receptors and the antigen

[TABLE]

• •

T cells get activated when they bind an APB, at a rate $W$

[TABLE]

We note that this is not the only mechanism for T cell activation, as the latter can also be driven by dendritic cells. We will show in Sec. 3 that the inclusion of dendritic cells has only a small quantitative effect on the model, hence, for simplicity, we will mostly negect it.

• •

Antigen presentation is a reversible process, that can be switched off [71], at a rate $\pi^{-}$

[TABLE]

Given that B cell activation requires antigen degradation, we will assume that antigen fragments freed up are not able to replicate, and do not contribute to the population of replicating antigens.

• •

Activated T cells induce B clonal expansion

[TABLE]

or contraction (via e.g. cytolysis [65, 66]),

[TABLE]

depending on the nature, helper or suppressor, of the T cell (we will introduce a binary variable to discern helpers and suppressors in the next section). We denoted $\lambda^{+}$ and $\lambda^{-}$, the rates of clonal expansion and contraction, respectively.

• •

B cells are kept in a state of activated apoptosis while undergoing clonal expansion in the germinal centre and compete for survival signal from T-helper cells [72]. We assume that when two B cells compete for the same resources, one will be selected at a rate $\delta$

[TABLE]

To keep the model simple, we will neglect several other processes such as antigen mutation, clonal expansion of T cells, activation of macrophages and other types of immune cells, production of antibodies by B cells, and affinity maturation of B cells in the germinal centre. We will simply assume that, due to the latter process, the affinity $\pi^{+}$ between BCR and antigen is an increasing function of time.

Despite the many simplifying assumptions, we believe that the above reactions capture the basic principles of an immune response to an antigen. In the next section, we will use them to build an analytically tractable model that links the microscopic dynamics of cellular activation with the macroscopic dynamics of clonal concentrations. We will assume that T and antigen presenting B cells interact via a network, while all the other species are well-mixed within the lymphatic system. In contrast to more traditional phenomenological approaches formulated in terms of differential equations, the resulting model will incorporate information on microscopic interactions and stochasticity, and will be formulated in terms of a minimal number of equations (i.e. four), which do not require the estimation of a large number of unknown parameters and are extremely robust to parameters variation. Results will highlight the role played by the ratio between T-helper and T-suppressor lymphocytes on immunosuppression, and will expose the role played by stochasticity on homeostasis. Both the T-helper/T-suppressor ratio and the mechanics of homeostasis have been the subject of intensive experimental investigation in recent years.

3 A network model of the adaptive immune system

We consider a population of T cells, each labelled by $i=1,\ldots,N$, and a population of B clones333A B clone is the ensemble of all B cells that have the same receptors and thus respond to the same antigen., each labelled by $\mu=1,\ldots,P$. Experimental lymphocytes counts in humans estimate the total number of T cells to be ${\mathcal{O}}{(10^{11})}$ [73], the total number of B cells to be roughly of the same order [74], and the number of B clones to be ${\mathcal{O}}{(10^{7})}$ or higher [75]. Hence, we set the number of B cells to be $\phi N$ with $\phi={\mathcal{O}}{(1)}$ and the number of B clones to be sub-extensive in $N$, i.e. $P=\alpha N^{\gamma}$ with $\gamma\in(0,1)$ and $\alpha={\mathcal{O}}{(1)}$, so that each B clone contains, on average, $N^{1-\gamma}\phi/\alpha$ cells. 444Clone sizes are known to be heterogeneous [76], here the average over all clones is given. For simplicity, we ignore cross-reactivity effects between B clones and antigens and we assume that each B clone $\mu$ is able to bind only to one type of antigen (assumed to have a single epitope) also labelled by $\mu$.

Next, we define the statistics of the interactions between T cells and B clones (see Fig. 1 for a schematic representation). We introduce a binary variable $\xi_{i}^{\mu}\in\{1,0\}$ to indicate whether ($1$) or not ([math]) T cell $i$ can bind cells in B clone $\mu$. We regard the variables $\{\xi_{i}^{\mu}\}$ as random i.i.d., with distribution

[TABLE]

where $c_{\mu}={\mathcal{O}}{(N^{0})}~{}\forall~{}\mu$, so that each B clone $\mu$ has, on average, ${\langle}\sum_{i=1}^{N}\xi_{i}^{\mu}{\rangle}=c_{\mu}N^{1-\gamma}$ conjugate T cells. If all the T cells that react with a B clone belong to the same T clonotype, then $c_{\mu}N^{1-\gamma}$ is the average size of the T clonotype conjugate to B clone $\mu$, at rest. However, this assumption is not necessary and we admit the possibility that a B clone interacts with different T clonotypes. In this case the $c_{\mu}N^{1-\gamma}$ T cells interacting with B clone $\mu$ may belong to different T clonotypes. In either case, choice (8) corresponds to the scenario where the number of T cells that can signal a B clone is of the same order as the number of B cells in that clone. Biologically, this would seem the most plausible scaling, as it avoids, on the one hand, a redundancy of B cells that would not get sufficiently signaled by a comparatively small number of conjugate T cells, and on the other hand, a waste of T cells in signaling a comparatively small number of B cells.

We denote with $\psi_{\mu}$, $b_{\mu}$ and $p_{\mu}$, the population densities of antigens of type $\mu$, B clone $\mu$ and APB of type $\mu$, respectively. It is convenient to define clonal densities as the number of cells in a clone divided by the average resting number of conjugate T cells, i.e. $b_{\mu}\!=\!B_{\mu}/c_{\mu}N^{1-\gamma}$, where $B_{\mu}$ is the number of cells in B clone $\mu$ and similarly for $p_{\mu}\!=\!P_{\mu}/c_{\mu}N^{1-\gamma}$ and $\psi_{\mu}\!=\!\Psi_{\mu}/c_{\mu}N^{1-\gamma}$. We can write the following equations for reactions (1), (2) and (4)

[TABLE]

with $p_{\mu},b_{\mu}$ and $\psi_{\mu}$ dimensionless and the kinetic coefficients $\pi_{\mu}^{+},\pi_{\mu}^{-},r_{\mu}$ having the unit of inverse time. In the absence of B cells, (10) gives $r_{\mu}=\ln 2/t_{\mu}^{\star}$, where $t_{\mu}^{\star}$ is the doubling time of antigen $\mu$. Doubling times vary across different diseases, however, typical values in human hosts, at the early stage of an infection, when an immune response has not started yet, are estimated to range between a few hours (bacteria) and a day (virus), see e.g. [77, 78]. It is then convenient to measure time in days, so to work with ${\mathcal{O}}{(1)}$ replication rates.

T cells can be helpers or suppressors. For each T cell $i$, we introduce a variable $\eta_{i}$ which takes values $1$ if $i$ is helper and $-1$ if it is suppressor. We assume each $\eta_{i}$ to be identically and independently sampled from

[TABLE]

The parameter $-1<\epsilon<1$ quantifies the imbalance between T-helper ($\eta=1$) and T-suppressor ($\eta=-1$) cells and is directly related to the T-helper/T-suppressor ratio $R$, measured in experiments

[TABLE]

For $\epsilon=0$, T-helper and T-suppressors cells are present in equal proportions, both equal to $1/2$ (i.e. $R=1$), while for $\epsilon>0$ ($\epsilon<0$) T-helper cells are more (less) than T-suppressor cells.

We represent the state of each T cell $i$ (active or inactive) with a variable $\sigma_{i}$ which takes value $1$ if $i$ is active and [math] otherwise. Helper T cells get activated via reaction (3), i.e. when their receptors bind to an APB. For simplicity, we assume that T cells update their state at regular time intervals of duration $\Delta$, according to the stochastic rule

[TABLE]

where $\theta(x)=0$ for $x\leq 0$ and $\theta(x)=1$ for $x>0$, and $z(t)$ is a zero-averaged random variable with suitably normalised variance, drawn, at each time $t$, from a symmetric distribution $p(z)$, which mimicks ”fast” noise in the biological environment or stochasticity in T cell activation. In the absence of noise ($z=0$), equation (13) tells that T cells activate in the presence of antigen presenting cells and become inactive in the absence of antigens. However, biological noise may occasionally lead the system to deviate from the expected behaviour. For example, some T cells may fail to activate even in the presence of their conjugate antigen presenting cell, or, conversely, some cells may activate randomly, even in the absence of antigens (for example due to activation by self-antigens in an auto-immune response). Such stochastic effects are modelled by the random variable $z$: positive values of $z$ damp the immune response while negative values enhance it. The time step $\Delta$ will eventually be sent to zero to retrieve the continuous time dynamics.

We note that we could easily account for activation of T cells mediated by dendritic cells by adding a term $\sum_{\mu}\xi_{i}^{\mu}d_{\mu}$ in the argument of the step function, with $d_{\mu}$ denoting the density of active dendritic cells, however, we will show later that this extra term has only a small effect, hence, for simplicity, we will neglect it here.

We note that equation (13) models reaction (3) in the presence of noise, at the microscopic level of individual T cells. Alternatively, one could model (3) at the population level, via reaction kinetics equations, similar to (9) and (10), for the densities of active and inactive $T$ cells, valid under the assumptions of well-mixed system and negligible fluctuations due to discreteness of cells. Noise could be included at population level, by introducing a reaction for spontaneous activation of $T$ cells of the type $T\to T^{\star}$ and one for spontaneous deactivation $T^{\star}\to T$, the rates of which would represent free parameters of the model. Our approach starts instead from stochastic equations for the microscopic cell states, which do not require the above assumptions and keep the number of free parameters to a minimum. Macroscopic cell densities as those involved in reaction kinetics approaches, can be obtained within our approach, as sums of microscopic variables, e.g. the density of active T cells binding APB $\mu$ can be obtained from $\sum_{i}\xi_{i}^{\mu}\sigma_{i}$, similarly the density of active helper T cells binding APB $\mu$ is obtained from $\sum_{i}\xi_{i}^{\mu}\sigma_{i}(1+\eta_{i})/2$.

Finally, B clones expand (contract) when they receive excitatory (inhibitory) signals from active T cells and compete for survival, so that each B clone follows a logistic dynamics

[TABLE]

Here, the first term represents clonal expansion, via reaction (5), triggered by active ($\sigma_{i}=1$) helper ($\eta_{i}=1$) T cells, specific for antigen $\mu$ ($\xi_{i}^{\mu}=1$). The second term accounts for clonal contraction, via reaction (6), triggered by active, specific suppressor ($\eta_{i}=-1$) T cells. The third and fourth term accounts for the binding and unbinding of B cells with antigens, respectively, via reaction (2) and (4). Since antigen binding typically engages only one (or at most a few) of the $n$ BCRs, the resulting APB is effectively still a B cell, with one spare receptor less. This leads to a decrease (increase) of the effective number of B cells by a fraction $1/n$ only, upon binding (unbinding) an antigen. The last term, a loss term proportional to the square of B cells population density, accounts for competition betweeen B cells. Different molecular mechanisms of clonal suppression have been described with no consensus on one universal mechanism [79]. For mathematical simplicity, we will assume that clonal suppression takes place at the same rate as clonal expansion, and set $\lambda_{\mu}^{+}=\lambda_{\mu}^{-}=\lambda_{\mu}$. Generalizations to the case $\lambda_{\mu}^{+}\neq\lambda_{\mu}^{-}$ are straightforward and we will comment on them later. This simplifies (14) to

[TABLE]

where we have neglected ${\mathcal{O}}{(n^{-1})}$ terms, bearing in mind that $n={\mathcal{O}}{(10^{5})}$, and defined the density of the net excitatory signal received by B clone $\mu$ from T cells as

[TABLE]

with {\mbox{\boldmath\sigma}}=(\sigma_{1},\ldots,\sigma_{N}) representing the microstate (active or inactive) of all T cells.

3.1 Macroscopic dynamics

The dynamics of the immune system model defined above entails the stochastic update (13) of the microstate $\sigma$ of all T cells. However, cell concentrations are seen to depend on $\sigma$ only through the variables {\bf m}({\mbox{\boldmath\sigma}})=(m_{1}({\mbox{\boldmath\sigma}}),\ldots,m_{P}({\mbox{\boldmath\sigma}})), via (15). In this section, we will use non-equilibrium statistical mechanical techniques [80, 81, 82] to derive, from the law of the microscopic system $\sigma$, equations for the macroscopic variables {\bf m}({\mbox{\boldmath\sigma}}). To this purpose, we write below the master equation for the probability density p_{t}({\mbox{\boldmath\sigma}}) to observe a microstate $\sigma$ at time $t$, and we will derive from this a dynamical equation for the macroscopic variables {\bf m}({\mbox{\boldmath\sigma}}).

Denoting by ${\mathcal{P}}(z\leq x)=\int_{-\infty}^{x}dz\,p(z)$ the cumulative distribution function of the noise distribution $p(z)$, the likelihood to observe configuration $\sigma_{i}$ at time $t+\Delta$, for any symmetric distribution $p(z)=p(-z)$, is [83]

[TABLE]

A natural choice for $p(z)$ would be a Gaussian distribution, with variance $\beta^{-1}$, which leads to ${\mathcal{P}}(z\leq x)={\frac{1}{2}}(1+{\rm erf}(\beta x/\sqrt{2}))$. An alternative choice is the so-called Glauber distribution leading to

[TABLE]

which is qualitatively very similar to the cumulative distribution function for the Gaussian distribution and is easier to work with analytically. The parameter $\beta$ has to be interpreted as an inverse noise level: for $\beta\to 0$, the dynamics (13) is fully stochastic, whereas for $\beta\to\infty$, the dynamics is deterministic. For the choice (17), the probability that helper T cell $i$ changes, in a single time step, its state $\sigma_{i}$ at time $t$ (to $1-\sigma_{i}$ at time $t+\Delta$) is

[TABLE]

where we used $1-2\sigma_{i}=\pm 1$ and $\tanh(\pm x)=\pm\tanh x$. Assuming that the update of T cells is sequential, i.e. at each time step one helper cell $i$, drawn at random, is updated with likelihood $W_{t}(\sigma_{i})$, one obtains, for $\Delta=1/N$ and $N$ large, the following master equation for the probability density p_{t}({\mbox{\boldmath\sigma}}) to observe microstate $\sigma$ at time $t$,

[TABLE]

where $F_{i}$ is a “cell-flip” operator that changes the configuration of T cell $i$ from $\sigma_{i}$ to $1-\sigma_{i}$ and has no effect on any other cell $j\neq i$. From (19) one can derive equations of motion for expectations {\langle}\cdot{\rangle}=\sum_{\mbox{\boldmath\sigma}}\cdot\,p_{t}({\mbox{\boldmath\sigma}}). Multiplying (19) by $\sigma_{j}$ and summing over $\sigma$, we obtain the rate of change

[TABLE]

of the average activity of T cell $j$, intuitively given by its variation $1-2\sigma_{j}$ upon a single cell flip $F_{j}$, times the rate $W_{t}(\sigma_{j})$ at which the cell is flipped.

Then, multiplying (20) times $\eta_{j}\xi_{j}^{\mu}$, summing over $j$, dividing by $c_{\mu}N^{1-\gamma}$, and using a mean-field approximation (see A), we obtain the following equation of motion for the average signal strength m_{\mu}(t)={\langle}m_{\mu}({\mbox{\boldmath\sigma}}(t)){\rangle} on clone $\mu$:

[TABLE]

In the above {\langle}\dots{\rangle}_{\eta,{\mbox{\boldmath\xi}}} denotes the average \sum_{\eta,{\mbox{\boldmath\xi}}}\dots P(\eta,{\mbox{\boldmath\xi}}) over the distribution of regulatory patterns in the system

[TABLE]

where $\delta_{x,y}=1$ for $x=y$ and $\delta_{x,y}=0$ otherwise, and {\mbox{\boldmath\xi}}_{j}=(\xi_{j}^{1},\dots,\xi_{j}^{P}) encodes the regulatory interactions between B cells and T cell $j$. Assuming that P(\eta,{\mbox{\boldmath\xi}})=P(\eta)P({\mbox{\boldmath\xi}}) i.e. the ability of a helper cell $i\,$ to bind to clone $\mu$ does not depend on whether $\,i\,$ is a helper or regulator, we have

[TABLE]

Finally, averaging over $\xi^{\mu}$, using the independence of the $\xi^{\mu}$’s and ${\langle}\xi^{\mu}{\rangle}=c_{\mu}/N^{\gamma}$, we get

[TABLE]

where one has still to carry out the average over the $\{\xi^{\nu}\}$ other than $\xi^{\mu}$. The second term in the round brackets represents clonal interference, i.e. contributions to clone $\mu$ arising from different clones $\nu$’s, due to the ability of T cells to bind different B clones, as illustrated in Fig. 1. Due to the specificity of the interactions between B and T cells, however, the $\xi$’s are extremely diluted, i.e. the probability for each $\xi^{\mu}$ to be non-zero is ${\mathcal{O}}{(N^{-\gamma})}$. Since only non-zero $p_{\nu}$’s contribute to clonal interference, as long as the number of different antigenic threats in the host is finite, the sum over $\nu$ consists of a finite number of terms, each ${\mathcal{O}}{(N^{-\gamma})}$, and vanishes in the thermodynamic limit. Clonal interference becomes instead ${\mathcal{O}}{(1)}$ when the number of different antigens in the host is ${\mathcal{O}}{(N^{\gamma})}$, i.e. of the same order as the number of B clones. In the remainder of the paper, we will focus on the immune response when a single antigen is present in the host, so cross-reactivity effects between B and T cells can be mostly neglected, however we will see in Sec. 5.2, that they may cumulate with cross-reactivity effects betweeen B cells and antigens, when present.

In A we show that for $\gamma<1$ fluctuations of m_{\mu}({\mbox{\boldmath\sigma}}) about its thermodynamic average $m_{\mu}$ vanish for large $N$. In this regime, the mean-field approximation becomes exact and we are allowed to replace m_{\mu}({\mbox{\boldmath\sigma}}) in (15) with $m_{\mu}$

[TABLE]

which enables us to express the dynamical evolution of our model in terms of a closed set of first order differential equations, namely (9), (10), (24) and (25).

We conclude this section by noting that B clonal dynamics (25) has one fixed point at $b_{\mu}=\lambda_{\mu}m_{\mu}/\delta_{\mu}$, which is stable for $m_{\mu}>0$, and one fixed point at $b_{\mu}=0$ which is stable for $m_{\mu}\leq 0$, suggesting that the immune system must keep a basal activity ($m_{\mu}>0$), even in the absence of antigens, to sustain cell numbers. In the absence of antigens, non-zero values of (16) are achieved via stochastic fluctuations in T cell activation. Such activity, whereby T cells keep sending survival signals to B clones in the absence of antigens, is experimentally observed and is believed to be one of the mechanisms to accomplish a homeostatic control of cell numbers [84]. In addition, lymphocyte homeostasis requires that clonal expansion during the immune response is matched by a comparable decrease in lymphocyte numbers once the antigen is removed. Its mechanics is not fully understood, and is believed to require a sophisticated regulation of the rates of cellular proliferation and programmed cell death. In this model, however, clonal contraction upon antigen removal, will emerge as a simple consequence of stochasticity in T cell activation, as we will see in the next section. As the antigen concentration is decreased, stochasticity effects become more important and eventually dominate the dynamics (13), thus restoring the activity (16) to its basal value, which determines the equilibrium sizes of B clones.

3.2 Response to a single antigen

In this section, we study the activation of the immune system when a single antigen of type $\nu$ is present, i.e. $\psi_{\nu}\neq 0$ and $\psi_{\mu}=0~{}\forall~{}\mu\neq\nu$. We assume that the system is initially at equilibrium in the absence of antigens, when antigen $\nu$ is introduced.

The equilibrium state in the absence of antigen is easily found by noting that $\psi_{\mu}=0$ is a fixed point of the $\mu$-antigen dynamics (10), and equation (9) implies that at stationarity the density of B cells presenting antigen $\mu$ is

[TABLE]

where $a_{\mu}=\pi^{+}_{\mu}/\pi_{\mu}^{-}$ can be thought of as a measure of the affinity between B clone $\mu$ and antigen $\mu$.555The affinity can be experimentally measured as the ratio between the equilibrium concentration of APB and the product of the equilibrium concentrations of antigen and conjugate B cells. Hence, at equilibrium, in the absence of antigen, one has $\psi_{\mu}=0$, $p_{\mu}=0~{}\forall~{}\mu$ and, from (24), $m_{\mu}=\epsilon/2~{}\forall~{}\mu$. Then (25) yields $b_{\mu}=\kappa_{\mu}\epsilon/2$, with $\kappa_{\mu}=\lambda_{\mu}/\delta_{\mu}$, for $\epsilon>0$, and $b_{\mu}=0$ for $\epsilon<0$.

Upon introducing antigen $\nu$, equation (24) gives, for any $\mu\neq\nu$

[TABLE]

showing that, for large $N$, $m_{\mu}=\epsilon/2$ is a stable fixed point for all clones $\mu\neq\nu$ not responding to antigen $\nu$, which thus permane in the state

[TABLE]

In contrast, for the $\nu$-clone, conjugate to the antigen introduced in the host, (24) gives

[TABLE]

which has to be solved together with

[TABLE]

leading to a fourth order dynamical system, with Jacobian

[TABLE]

The $\nu$-clone will return to its resting state, meaning that the antigen will be cleared, if the state $(\psi_{\nu},p_{\nu},m_{\nu},b_{\nu})=\left(0,0,\frac{\epsilon}{2},{\rm max}\left\{0,\frac{\kappa_{\nu}\epsilon}{2}\right\}\right)$ is a stable fixed point of the $\nu$-clone dynamics (28, 29, 30, 31). This requires the eigenvalues of the Jacobian, evaluated at the resting state, to be negative. For $\epsilon>0$, the resting state is $\left(0,0,\frac{\epsilon}{2},\frac{\kappa_{\nu}\epsilon}{2}\right)$ and the eigenvalues

[TABLE]

are all negative for $\epsilon>2r_{\nu}/(\pi_{\nu}^{+}\kappa_{\nu})$, meaning that the antigen will be cleared when the affinity $\pi_{\nu}^{+}$ of the conjugate B cells is greater than a value that increases as the pathogen replication rate increases and $\epsilon$ decreases. For $\epsilon<0$, the Jacobian must be evaluated at the state $\left(0,0,\frac{\epsilon}{2},0\right)$ and the corresponding eigenvalues

[TABLE]

show an instability in the antigen dynamics. It is easy to see from (28) and (31), that in this regime $m_{\nu}$ will converge to negative values, as $\tanh(x)\geq 0~{}\forall~{}x\geq 0$, and $b_{\nu}$ will vanish in the long-time limit, leading to an immuno-suppressed host, where the antigen grows indefinitely according to

[TABLE]

Finally, for $0<\epsilon<2r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+})$ the fixed point $(\psi_{\nu},p_{\nu},m_{\nu},b_{\nu})=\left(0,0,\frac{\epsilon}{2},\frac{\kappa_{\nu}\epsilon}{2}\right)$ is unstable and different types of dynamics may arise, depending on the range of the kinetic parameters. Equations (28), (29) and (31) show that, in this regime, $m_{\nu}$ will approach positive values in the long-time limit and $b_{\nu}$ and $p_{\nu}$ will thus evolve towards $\kappa_{\nu}m_{\nu}$ and $a_{\nu}\kappa_{\nu}m_{\nu}\psi_{\nu}$, respectively. Hence, in the long-time limit, the dynamics of clone $\nu$ can be described in terms of the second order dynamical system

[TABLE]

This has a fixed point at

[TABLE]

whose stability is determined from the eigenvalues

[TABLE]

where

[TABLE]

We have

[TABLE]

as $x\,{\rm atanh}x-{\rm atanh}x+1>0~{}\forall~{}x$, hence both eigenvalues $\lambda_{1,2}$ have negative real part for $A_{\nu}>0$ i.e. for $\epsilon>r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+})$. Therefore, in the regime $r_{\nu}/(\pi_{\nu}^{+}\kappa_{\nu})<\epsilon<2r_{\nu}/(\pi_{\nu}^{+}\kappa_{\nu})$, the fixed point $(m_{\nu}^{\star},\psi_{\nu}^{\star})$ is stable and the antigen will approach the non-zero value (43), which is higher, the higher the noise level $\beta^{-1}$, while $m_{\nu}$, $b_{\nu}$ and $p_{\nu}$ will approach the $\beta$-independent values $r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+})$, $r_{\nu}/\pi_{\nu}^{+}$ and $a_{\nu}r_{\nu}\psi_{\nu}^{\star}/\pi_{\nu}^{+}$, respectively. Conversely, for $0<\epsilon<r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+})$ the fixed point $(m_{\nu}^{\star},\psi_{\nu}^{\star})$ is unstable. It is easy to see from (41), that for long times $0\leq m_{\nu}\leq\epsilon$, as $-1\leq\tanh x\leq 1~{}\forall~{}x$, then (42) implies that when $m_{\nu}$ reaches its maximum value $\epsilon$, the viral concentration will keep growing in the regime $\epsilon<r_{\nu}/(\pi_{\nu}^{+}\kappa_{\nu})$.

Numerical solutions of the full dynamical system (28, 29, 30, 31) are plotted in Figure 2 for $r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+})<\epsilon<2r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+})$ and in Figure 3 for $0<\epsilon<r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+})$, and confirm these predictions, which can also be verified by evaluating numerically the eigenvalues of the full dynamical system, in these regimes.

In conclusion, four different regimes are identified: (i) $\epsilon>2r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+})$ where the antigen is cleared and the responding B clone returns to its resting size; (ii) $r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+})<\epsilon<2r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+})$ where the viral concentration evolves to a non-zero value, which, interestingly, depends on the noise level. In this regime the antigen permanes indefinitely in the host, although its growth is limited by the action of the immune system. The responding B clone remains in an expanded state $m_{\nu}^{\star}>\epsilon/2$ and fails to contract and return to its resting size, even at large times; (iii) $0<\epsilon<r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+})$ where the immune system responds to the antigen by expanding the conjugate B clone, but the response is too weak and the antigen proliferates indefinitely in the host; (iv) $\epsilon<0$ where the immune system is irresponsive and B cells decrease over time, while the antigen proliferates indefinitely. Finally we note that including the contribution from dendritic cells to T cells activation, would simply add a term in the argument of the hyperbolic tangent in (28). Since $-1\leq\tanh x\leq 1~{}\forall~{}x$, the inclusion of dendritic cells would not alter the different phases and would only have a small quantitative effect on the transient response.

3.2.1 The role of affinity maturation

We have so far regarded the binding rate between antigen and B cells $\pi_{\nu}^{+}$ as a constant, however, as B cells undergo clonal expansion, they increase their affinity with the antigen by several orders of magnitude [85], via affinity maturation, so $\pi_{\nu}^{+}$ is an increasing function of time [62]. Several studies have modelled in detail B cell maturation in the germinal centre, over the last few decades [86, 87, 88]. Here, we simply assume that the affinity is an increasing function of time that saturates at large time, and choose $\pi_{\nu}^{+}(t)$ as the sigmoid function

[TABLE]

where $v,t^{\star},\pi_{M}$ and $\pi_{0}$ are parameters that control, respectively, how fast, early and large the affinity grows and its initial value. Experimental findings suggest that high affinity B cells bind their target antigen within a few minutes [89], hence we estimate $\pi_{\nu}^{+}(\infty)=\pi_{M}(\pi_{0}+1)$ to be ${\mathcal{O}}{(10^{2})}/{\rm day}$. We choose the remaining parameters in such a way that the ratio $\pi_{\nu}^{+}(\infty)/\pi_{\nu}^{+}(0)=(\pi_{0}+1)/(\pi_{0}-\tanh vt^{\star})$ between high and low affinity is within the estimated physiological range ${\mathcal{O}}{(10^{5})}$ [90]. It is also suggested that birth and death rates of B cells in the germinal centre are of the same order of magnitute, yielding $\kappa_{\nu}={\mathcal{O}}{(1)}$ [91]. This implies that regimes (ii) and (iii) take place on a narrow range of values of $\epsilon\in(0,2r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+}(\infty)))$: above this range, the antigen is cleared, while below it, the immune system is irresponsive and the antigen replicates indefinitely in the host, no matter how fast or large $\pi_{\nu}^{+}$ grows. This results in a rather abrupt transition of the immune system from a responsive to a suppressed state, as $\epsilon$ approaches zero from above, in line with experimental findings.

Numerical solution of (28), (29), (30) and (31), in the regime $\epsilon\!>\!2r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+}(\infty))$ are shown in Figs. 4, 5, 6 and 7. Plots of $b_{\nu}$ and $\psi_{\nu}$ versus time, are shown for different choices of the inverse noise level $\beta$ (Fig. 4), time-dependence of $\pi_{\nu}^{+}$ (Fig. 5), kinetic parameters $\lambda_{\nu}$, $\delta_{\nu}$ (Fig. 6), and $\pi_{\nu}^{-}$, $r_{\nu}$ (Fig. 7).

Plots show an immune response that follows closely the behaviour of the antigen, i.e. it expands B clones while the antigen concentration is increasing and contracts them when the antigen concentration is decreasing, so that both B clone and antigen concentrations are unimodal functions of time, meaning that the system is able to accomplish homeostasis. Interestingly, the time-dependent $b_{\nu}$ concentration shows an initial plateau, which is consistent with the experimentally observed lag-time between an infection and a detectable immune response [79]. The stochastic noise $\beta^{-1}$ has the effect of mildly reducing the height of the peak in B clones concentration, thus increasing the peak in viral concentration, however the system will be able to remove the antigen even at high noise levels (see Fig. 4). The time-dependence of $\pi^{+}_{\nu}$, affects both the location and the height of the peaks, consistently with the intuition that the faster the affinity grows, the earlier and the smaller the peak in viral concentration. Fig. 5 shows that affinities increasing faster with time ($t^{\star}\!=\!7$, $v\!=\!1$) outperform those that increase slower ($t^{\star}\!=\!20$, $v\!=\!0.1$). The kinetic parameters $\lambda_{\nu}$ and $\delta_{\nu}$ have a mild effect on the shape of the B clonal concentration (see Fig. 6), while Fig. 7 (top panel) shows that $\pi_{\nu}^{-}$ affects the decay of B cells after the infection peak, potentially sustaining long-term memory. Fig. 7 (bottom panel) shows that although faster replicating antigens attain much higher concentrations, they are still removed, for $\epsilon\!>\!2r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+}(\infty))$, on similar timescales, by triggering a stronger immune response. Finally, Fig. 8 shows that at the critical value $\epsilon\!=\!2r_{\nu}/(\kappa_{\nu}\pi_{\nu}^{+}(\infty))$ the model predicts a very slight increase in the B cell population and successful clearance of the antigen (top panels), however stochastic fluctuations are anticipated to become important at criticality, as Gillespie simulations in Sec. 4 will confirm. Conversely, as soon as $\epsilon$ is lowered below zero, the system becomes unresponsive and is not able to fight a single antigen, even if replicating slowly (Fig. 8, bottom panels).

These results suggest that for T-helper/T-suppressor ratios $R>1$ (i.e. $\epsilon>0$) the host manages to remove completely the antigen, provided the affinity grows larger than an $\epsilon$-dependent value $\pi_{\nu}^{+}>2r_{\nu}/(\epsilon\kappa_{\nu})$, while for $R<1$ the immune system is impaired and does not respond to the antigen, which replicates indefinitely in the host, no matter how large and fast the affinity $\pi_{\nu}^{+}$ between BCR and antigen grows. From a more general point of view, the model captures the effectiveness and robustness of the immune system dynamics after exposure to an antigen and predicts that the most important single parameter in determining whether the system is in a healthy or in an immuno-suppressed phase is $\epsilon$, directly related to the T-helper/T-suppressor ratio $R$, and in line with recent experiments [92, 51, 53, 57]. The model predicts the ratio to affect directly the location of the peak of the time-dependent antigen concentration, occurring when the affinity $\pi_{\nu}^{+}$ becomes larger than the $\epsilon$-dependent value $2r_{\nu}/(\epsilon\kappa_{\nu})$, as well as the resting sizes of B clones, related to $\epsilon$ via $b_{\nu}=\epsilon\kappa_{\nu}/2~{}\forall~{}\nu$. This result is consistent with the equilibrium statistical mechanical analysis carried out in [93], where the B cells density was shown to be a decreasing function of the T-suppressor cells density, at equilibrium. Note that $B_{\nu}=c_{\nu}N^{1-\gamma}b_{\nu}$ so the equilibrium size of B clone $\nu$ will depend on the number $c_{\nu}N^{1-\gamma}$ of T cells signaling B clone $\nu$. B clones which can be signaled by different T clones, will have higher values of $c_{\nu}$ and will therefore reach larger sizes at stationarity.

Finally we note that, although unstable, $b_{\nu}=0$ is always a fixed point of (31), even for $\epsilon>0$ and may be selected by means of stochastic fluctuations, leading to clonotypes extinction. Stochastic fluctuations are predicted to become relevant in the regime $\gamma=1$, where the number of B clones is extensive and the number of cells per clonotype is ${\mathcal{O}}{(N^{0})}$, a scenario which has been observed in mouse [94].

4 Simulations

In this section we employ the Gillespie algorithm [95], to simulate the biochemical reactions defined in Sec. 2 at the single cell level, in the presence of a single antigen population. The underlying assumption of the algorithm is that cell populations are well-mixed and interact in a finite volume. Since we do not consider antigen mutations and cross-reactivity effects between B clones and antigen, we will only consider one B clone and the ensemble of its conjugate T cells, so that we drop clonal indeces from now on.

We have seven different species:

• •

B cells (B)

• •

Antigen presenting B cells (APB)

• •

Active and inactive suppressor T cells (S/S*)

• •

Active and inactive helper T cells (H/H*)

• •

Antigen (Ag)

and ten different reactions:

• •

Antigen replication at rate $r$ (1)

• •

Binding between antigen and B cell, at rate $\pi^{+}$ (2)

• •

Activation and deactivation of T cells, at rate $W=W_{t}(0)$ and $W^{\prime}=W_{t}(1)$, respectively (see (18))

[TABLE]

We note that the activation of T cells due to binding with APB cells, is modelled here as a single-specie reaction with APB-dependent rate. Alternatively, one could model it as a two-specie reaction (3) with APB-independent rate, and account for spontaneous activation of T cells, due to noise, separately, as a single-specie reaction with constant rates. We have checked that the two implementations are equivalent, however the above implementation has the advantage of reducing the number of reactions and connects more explicitely with the equations of Sec. 3, meaning that the same reaction rates can be used.

• •

Unbinding of B cells and antigen, at rate $\pi^{-}$ (4)

• •

Expansion and contraction of B cells

[TABLE]

• •

B cells competition, at rate $\delta$ (7)

The initial conditions of the model are chosen from a well-mixed system in equilibrium at inverse noise level $\beta=1$, in the absence of antigen. The number of cells in the T clone $cN^{1-\gamma}$ is set to $10^{4}$, of which there are active and inactive T-helper and T-suppressor cells. Taken together, the initial total T-helper and total T-suppressor populations, $N_{H}=H+{H^{\star}}$ and $N_{S}=S+{S^{\star}}$, are given by $N_{H}=cN^{1-\gamma}(1+\epsilon)/2$ and $N_{S}=cN^{1-\gamma}(1-\epsilon)/2$. To find the ratio of active to inactive T cells, we refer to the update rule (13). Given that $\beta=1$ and that the APB concentration is initially zero, we have that half of T cells are initially active and half inactive. This approach of deriving initial T cell populations is equivalent to finding the steady state solution to the ODE (27) for $m_{\mu}$ and equating it to equation (16). Finally, the steady state value of B cells is obtained from equation (25), by setting $m_{\mu}=\epsilon/2$, which gives $B=cN^{1-\gamma}\lambda\epsilon/(2\delta)$. We denote with ${\bf n}=(B,APB,S,S^{\star},H,H^{\star},{\rm Ag})$ the population vector associated to the seven different species and we let ${\bf k}=(r,\pi^{+},W,W,W^{\prime},W^{\prime},\pi^{-},\lambda,\lambda,\delta)$ be the vector of the deterministic reaction rates. Deterministic reaction rates ${\bf k}$ from the ODE approach can be converted to stochastic reaction rates ${\bf c}$, following a simple rule [95]: $c_{s}=k_{s}$ if reaction $s$ involves a single reactant, $c_{s}=k_{s}/V$ for two distinct reactants, and $c_{s}=k_{s}/V^{n-1}$ for $n$ distinct reactants, where $V$ is the volume used to convert concentrations in the ODE approach to cell numbers in the stochastic simulation. This covers most possibilities except for when there are two reactants of the same species, in which case $c_{s}=2k_{s}/V$, to account for the combinatorics of reactions involving the same species.

Since concentrations in the ODE approach were defined in terms of clonal sizes over the average number of conjugate T cells $cN^{1-\gamma}$ (as opposed to volume $V$) we divide the deterministic reaction rates for reactions involving two species by $cN^{1-\gamma}$ and for reactions with the same single reactant appearing twice by $cN^{1-\gamma}/2$. This leads to ${\bf c}=(r,\pi^{+}/(cN^{1-\gamma}),W,W,W^{\prime},W^{\prime},\pi^{-},\lambda,\lambda,2\delta/(cN^{1-\gamma}))$. At each iteration, the algorithm calculates $a_{0}=\sum_{\ell=1}^{10}a_{\ell}$, where $a_{\ell}({\bf n})=c_{\ell}({\bf n})h_{\ell}({\bf n})$ denotes the propensity and $h_{\ell}({\bf n})$ the distinct combinations of reactants in reaction $\ell$, and it determines the next reaction $s$ to execute as the one satisfying $\sum_{\ell=1}^{s-1}a_{\ell}<r_{1}a_{0}\leq\sum_{\ell=1}^{s}a_{\ell}$ for a random number $r_{1}$ generated uniformly in $(0,1)$, and the time until it occurs via $\tau=a_{0}^{-1}\ln r_{2}^{-1}$ for a random number $r_{2}$ generated uniformly in $(0,1)$, hence it updates the number of molecules [95].

The results of Gillespie simulations are shown for B cells and antigen populations, in Figures 9, 10 and 11. They are in agreement with those presented in Sec. 3, confirming the validity of the reduced (mean-field) description of the system in terms of four ODEs, involving the macroscopic variable $m_{\nu}$, presented in Sec. 3. In particular, all the simulations are seen to converge to the predicted steady states, within finite size fluctuations ${\mathcal{O}}{(N^{(1-\gamma)/2})}$. Modifying the time-dependence of $\pi^{+}$ shifts the peak and affects the maximum value of the antigen population (Fig. 9), increasing $\lambda$ or decreasing $\delta$ increases the population of B cells but it leads to the same qualitative dynamics of antigen removal following B clonal expansion (Fig. 10), and lowering $\epsilon$ below zero, the B cell population becomes extinct and the antigen population diverges (Fig. 11, bottom panels), consistently with the ODE approach. Conversely, the top panels of Fig. 11 show that at criticality, different realizations show different behaviours, consistently with the expectation that fluctuations about the mean-field solution become important. This suggests that stochasticity in the biological processes and discreteness of cells may drive immune systems operating near criticality in the suppressed phase.

5 Model extensions

In this section we look at two extensions of the model, one allowing different rates for clonal expansion and suppression, and the other inluding cross-reactivity effects, which occur when a BCR $\mu$ directed against antigen $\mu$ is also successful in binding with another antigen $\nu\neq\mu$.

5.1 The role of clonal expansion and suppression rates

In the model defined in previous sections, we made the assumption $\lambda_{\nu}^{+}=\lambda_{\nu}^{-}$. Allowing different kinetic coefficients for clonal expansion and suppression $\lambda_{\nu}^{+}\neq\lambda_{\nu}^{-}$, modifies (25) to

[TABLE]

where t_{\nu}={\langle}t_{\nu}({\mbox{\boldmath\sigma}}){\rangle} and

[TABLE]

represents the density of activated T cells, regardless of their being helper or suppressors. Under the assumption of independence between $\eta$ and $\xi$, that was used to derive (24), one has from (16) and (54) $m_{\nu}=\epsilon t_{\nu}$. Substituting in (53), we get

[TABLE]

with

[TABLE]

Since $t_{\nu}\geq 0$ at all times, from (55) we have that in the long-time limit $b_{\nu}=t_{\nu}\tilde{\lambda}_{\nu}/\delta_{\nu}$ for $\tilde{\lambda}_{\nu}>0$ and $b_{\nu}=0$ for $\tilde{\lambda}_{\nu}\leq 0$, where $t_{\nu}$ is the stationary solution of

[TABLE]

Then, from (30), it follows that the antigen is cleared for $r_{\nu}<\pi_{\nu}^{+}(\infty)\tilde{\lambda}_{\nu}/(2\delta_{\nu})$, i.e. for

[TABLE]

Under the assumption that $\pi_{\nu}^{+}$ increases to values ${\mathcal{O}}{(10^{2})}/{\rm day}$, while the other kinetic parameters are ${\mathcal{O}}{(1)}$, the critical value of $\epsilon$ above which the antigen is removed from the system is approximately given by

[TABLE]

Recent experimental results suggesting that the onset of immunosuppression is associated with lymphocyte ratios close to one, i.e. $\epsilon\simeq 0$, justify a posteriori the assumption $\lambda_{\nu}^{+}\simeq\lambda_{\nu}^{-}$.

5.2 The role of cross-reactivity between B cells and antigens

Our analysis has so far restricted to single epitope antigens, however, antigens have normally several epitopes which may be recognized by different B clones, leading to cross-reactivity effects. Here, we discuss briefly how these can be incorporated in the model. We introduce a variable $A_{\mu\nu}$ which takes value $1$ if BCR $\mu$ can bind antigen $\nu$ and [math] otherwise. The model studied ealier, with no B-Ag cross-reactivity, corresponds to $A_{\mu\nu}=\delta_{\mu\nu}$, where BCR $\mu$ can only recognize antigen $\mu$. In the presence of cross-reactivity, when an antigen $\nu$ enters the system, all the clones $\mu$ such that $A_{\mu\nu}=1$, will respond, so one has

[TABLE]

with each B clone $\mu$ evolving according to equation (25). We model the interactions $\{A_{\mu\nu}\}$ between B clones and antigen $\nu$, as random variables with distribution

[TABLE]

where $d_{\nu}={\langle}\sum_{\mu}A_{\mu\nu}{\rangle}$ is the average number of B clones reacting with antigen $\nu$, which we assume ${\mathcal{O}}{(1)}$. In the long-time limit, APB densities will approach the value $p_{\mu}=A_{\mu\nu}b_{\mu}\psi_{\nu}a_{\mu}$ and $m_{\mu}$ will evolve according to

[TABLE]

The sum on the right end side is due to clonal interference, now comprising two effects: cross-reactivity between B and T cells and cross-reactivity between B cells and antigens. If both types of cross-reactive interactions ($\xi$’s and $A$’s) are diluted, as postulated here, the sum is ${\mathcal{O}}{(N^{-\gamma})}$, as it consists of ${\mathcal{O}}{(N^{\gamma})}$ terms, each of order ${\mathcal{O}}{(N^{-2\gamma})}$. For $N\to\infty$, this vanishes and all the clones which are able to bind antigen $\nu$, of which there are $d_{\nu}$, will expand via (25) and

[TABLE]

Hence, cross-reactivity between B cells and antigens leads to the signaling of multiple B clones in parallel, but as long as cross-reactive interctions are diluted, they do not lead to signal interference. On the other hand, since multiple B clones now jointly contribute to the clearance of the antigen, the latter can be removed at smaller values of the affinities $\{\pi_{\mu}^{+}\}$, than those required in the absence of cross-reactivity.

Finally, we note that in the presence of B-Ag cross-reactivity, interclonal competition may arise, as cross-reactive B clones may compete for the same resources during clonal expansion, leading to coupled B clones dynamics $\dot{b}_{\mu}=b_{\mu}(\lambda_{\mu}m_{\mu}-\delta_{\mu}\sum_{\rho}A_{\mu\rho}b_{\rho})$. Although a full analysis of cross-reactivity effects goes beyond the scope of this work, the above discussion suggests that cross-reactivity may allow expansion of multiple clones in parallel and antigen clearence at smaller values of the affinity, on the other hand interclonal competition may arise and couple the dynamics of different clones, potentially resulting in time-variation of clones which cannot bind to the invading antigen directly.

6 Conclusions

In this work, we introduced a statistical mechanical model for the adpative immune system, which comprises B cells, T-helper cells, T-suppressor cells and antigens. When the ratio between T-helper and T-suppressor cells is above one, the model produces an immune response which is a unimodal function of time, in qualitative agreement with experimental observations, accomplishing, in particular, a complete removal of the antigen as well as lymphocyte homeostasis, for any chosen value of the control parameters (replication rates, cell death rates and noise level) and any chosen increasing function of time for the affinity between B cells and antigens, provided it increases beyond a critical value. The model correctly predicts the existence of a lag time between infection and immune response detection, and informs on the role played by the T-helper/T-suppressor ratio and by different kinetics parameters, on the steady state and relevant timescales (e.g. antigen density peak and B clonal contraction), which could be potentially useful for parameters inference.

As the ratio between T-helper and T-suppressor cells is tuned over a narrow region above one, the model exhibit a transition from a functional to an impaired phase, thus supporting the validity of the T-helper/T-suppressor ratio as an index of immuno-suppression. In the transitional state between the functional and the impaired state, noise plays an important role and B clones fail to contract after the antigen concentration reaches its peak, a pattern which is observed in ageing [96].

The identification of reliable markers of immuno-suppression is currently an active research field and several indices have been correlated to disease prognosis in recent years (e.g. CD4+/Treg, Treg/CD8+, Th17/Treg, CD4+/CD4+CD25+, CD4+/CD8+ ratios, absolute numbers, and differences in cell numbers, of different sub-populations). With more lymphocyte sub-populations being uncovered, it becomes increasingly important to have models which are able to discern relevant parameters from accidental correlations. The T-helper/T-suppressor ratio $R$ considered in this work, is given by the number of T-helper cells, i.e. CD4+ cells which are not regulatory, divided by the total number of CD8+ and T-regulatory cells, so $R=({\rm CD4}-{\rm Treg})/({\rm CD8}+{\rm Treg})$. This can be expressed in terms of the ratios $R_{1}$=CD4+/CD8+ and $R_{2}$=Treg/CD4+ commonly reported in the literature, as $R=R_{1}(1-R_{2})/(1+R_{1}R_{2})$. The latter is an increasing function of $R_{1}$ and a decreasing function of $R_{2}$, consistently with the experimental finding that high values of the CD4+/CD8+ and the CD4+/Treg ratios both correlate with positive outcomes. This puts on a firmer ground the use of the CD4+/CD8+ and CD4+/Treg ratios for prognosis monitoring. We argue, however, that given Treg cells are a small percentage of CD4+ cells, they may be subject to larger fluctuations, so the CD4+/CD8+ ratio may be a more reliable index of immuno-suppression than CD4+/Treg. Also, we propose that the combination $R$ of the two may give further insights, especially when the ratio Treg/CD4+ is abnormally large, as in cancer, autoimmunity and HIV diseases, so that $R$ may deviate significantly from $R_{1}$. We stress, however, that spatial heterogeneities have not been taken into account in the model and may play an important role. In particular, cell concentrations and kinetics may vary considerably across different tissues and organs, so the main merit of this approach lies in providing a mathematical tool to assess the relevance of the different cell ratios which have recently attracted experimental attention, rather than in determining accurately the values at which immune impairment occurs, which may be subject to spatial fluctuations.

Possible pathways for future research include the effect of fast mutating antigens like cancerous cells and retrovirus, that manage to mutate before being removed from the system, and the mechanism by which HIV infection, cancer progression and immunosenescence alter the T-helper/T-suppressor ratio. In addition, the biological realism of the model may be improved in several directions, by including interclonal competition and cross-reactivity effects, a more detailed modelling of the affinity maturation of B cells in the germinal centre, the role of self-antigens, T-T interactions and clonal expansion of active T cells. The latter process will lead to time variations of the number of T cells, which may be included in future extensions of the model by making the parameters $c_{\mu}$ dependent on time. The evolution of the size and diversity of T and B repertoires, encoded in the $c_{\mu}$ and $b_{\mu}$, in response to antigenic histories, may be an interesting pathway for future research.

We hope that even at this level of simplification, the model can offer a useful theoretical framework to understand homeostasis and impairment in the adaptive immunity, and complement recent experimental studies aimed at assessing the validity of the T-helper/T-suppressor ratio as a biological marker for immunosuppression.

Acknowledgements

AA acknowledges the stimulating research environment provided by the EPSRC Centre for Doctoral Training in Cross-Disciplinary Approaches to Non-Equilibrium Systems (CANES) (EP/L015854/1) and Dr Alexander Mozeika for many interesting discussions.