Mathematical Model to Predict Polyclonal T-Cell-Dependent Antibody Synthesis Responses

Abstract

Mathematical models are becoming indispensable tools to explore the complexities of biological systems at cellular levels. We present a model to explore the baseline immune cell interactions for in vitro polyclonal antibody synthesis via B-cells regulated by helper and regulatory T-cells. The model incorporates interactions of antigen-presenting cells, T-cells, regulatory T-cells, and B-cells with each other and predicts time-dependent trajectories of these cells and antibody synthesis stimulated by pokeweed mitogen. We used an ordinary differential equation-based approach to simulate the dynamic changes in the cells and cytokines numbers due to the cellular and humoral response to pokeweed mitogen stimulation. The parameters of the ordinary differential equations model are determined to yield a normal immune response as observed in the pokeweed mitogen-stimulated in vitro antibody synthesis via normal T, B, and antigen-presenting cells. The dose effects of antigen load and basal values of regulatory T-cells on the profiles of various immune response variables are also evaluated.

1. Introduction

Mathematical modeling is a powerful tool to explore complex interactions and mechanisms in immunoregulatory networks of cells and cytokines. The development of models for immune responses can provide insights into the immune system’s multi-layered, organized, interactive, and regulated network of cellular components and cytokines/chemokines. Such a model can potentially improve biologically and clinically relevant predictions for the diagnosis and treatment of human diseases as well as vaccine responses. Here, we present a model that functionally exhibits optimal T-cell-dependent antibody synthesis response via B-cells upon polyclonal stimulation, while incorporating the upstream complex interactions involved in the T-cell-dependent B-cell antibody response using pokeweed mitogen (PWM) as a nominal antigen.

Upon antigen recognition, cognate B-cells become activated and initiate proliferation and differentiation processes to produce antibody-secreting plasma cells (PCs) crucial in the initial response, and long-lived PCs and memory cells that contribute to sustained immunity. Earlier in vitro studies have shown polyclonal antibody synthesis via PWM-stimulated B-cells using unfractionated peripheral blood mononuclear cells (PBMC) or defined proportions of T- and B lymphocytes or T-cell subpopulations in the context of stem cell transplantation, infections, and vaccinations [1,2,3]. PWM-induced immunoglobulins (Ig) synthesis via B-cells [4,5,6,7] requires the presence of both monocytes acting as an antigen-presenting cell (APC) and helper T-cells acting as regulators of B-cell activation [8,9,10]. Proliferation and differentiation of B-cells, development of memory B-cells, and antibody affinity maturation are regulated by cytokines produced by T helper cells (T_h) and regulatory T-cells (T_reg) [7,11,12].

Given that PWM stimulation of B-cells provides an immunologically well-defined platform to investigate immune responses, the PWM system supported by historically acknowledged and robust experimental data [13,14,15,16,17] in normal PBMC was used to develop the model. This model incorporates the functional dynamics of the selected immune variables considered to be critical to dissect the process of T-cell receptor (TCR) activation leading to antibody production by B-cells. The model could eventually be applied to predict optimal immune response that can prevent infectious diseases such as SARS-CoV2 and infections and can be expanded to study a variety of liquid and solid tumors [18] and autoimmune diseases. Our simulations provide new insights into how the immune system responds to cell population dynamics under balanced and imbalanced immunologic states.

2. Methods

2.1. Process Description

The in vitro method for measuring I_g production uses T-cells, non-T-cells (comprising of B-cells, monocytes, natural killer cells, immune cell-derived cytokines, and PWM, while the numbers and proportions of T_h and T_regs are adjusted to regulate in vitro I_g production [19]. PWM was used as a polyclonal antigen (polyclonal Ag) trigger to stimulate, proliferate, and differentiate B-cells into antibody-producing PCs with the help of T-cells.

Immune response to Ag stimulation: Our objective is to define the proliferation and differentiation of PC under simulated immunologic conditions and quantify the biological processes. Figure 1 describes the immune components and their interactions leading to antibody production. The antigen is recognized by a naïve dendritic cell (DC) or APC and B-cells via their surface receptors and presented to naïve T-cells (T_n) that differentiate into T_h cells and T_reg. The activated T-cells and APC, in turn, produce IL-2, IL-4, and IL-6 cytokines [20] (T-cell growth, B-cell proliferation, and B-cell differentiation factors, respectively). The IL-4 and IL-6 enriched microenvironment promotes B-cell proliferation and differentiation into antibody-producing PC. To keep the model simple, we chose dominant cytokines that participate in B-cell proliferation, differentiation, and regulatory activity. The mathematical framework was developed taking into consideration the T-B-cells interaction, B-cell proliferation and differentiation, and I_g synthesis by PC.

2.2. Development of Mathematical Model

With the advances in data storage and computational power, the applications of mathematical modeling [21] to complex biological processes are gradually gaining popularity. The quantitative techniques can provide deep insights into cellular biology and other biological processes [22]. Over the last decades, mathematical modeling has also been seen as an alternative tool to dive deeper into the investigation of HIV infection [23] and to investigate spatio-temporal dynamics of tumor growth [24]. Tumor models have also been developed [25] to differentiate normal and cancer cell growth regimes using radially symmetric reaction–diffusion equations. Mathematical models are useful to generate time-dependent quantitative data of the cell populations, growth/decay rates, and maximum/minimum values. These models allow us to explore the dynamics of complex diseases in terms of their latent molecular mechanisms, help in treatment strategy optimization, and new drug development and discovery. Several modeling approaches have been applied for biological systems to simulate dynamic changes in the regulatory networks, disease progression, and cellular system. One of the widely used methods is based on setting up ordinary differential equations (ODEs) [26,27] for the biological variables. Within this approach, a first-order rate of change in a biological variable [28] is expressed in terms of the variables that affect its dynamics. Since the biological variables are dependent on each other, we obtain a set of coupled equations that are solved numerically to obtain the time-dependent continuous evolution of the variables. The ODE-based methods can be put into three classes: (1) Michaelis Menten Kinetics [29], (2) the law of mass action [30], and (3) the Hill function [31]. Each of these methods is suited for a specific biological process. In our study, we used the Michaelis Menten Kinetics approach to model continuous changes in populations of intercellular networks of various immune cells and cytokines.

Another mathematical approach is based on Petri nets formalism, which allows different abstraction levels of the biological system, from purely qualitative to more complex quantitative models. Mapping a biological system into an intuitive network representation, Petri nets [32], offers a convenient quantitative formalism to investigate biological systems. Agent-based modeling [33] is another quantitative technique that is being used to simulate a large set of active components represented by agents of biological systems. Within this approach, a biological system is represented via agents, which are the model components performing independent biological processes and interacting with other agents and the microenvironment [34].

In addition to the investigations of the dynamic evolution of the cytokines and chemokines as discussed in the above methods, their diffusive motion within a restricted spatial biological environment is also of vital interest [35,36] and can be discussed by space-time spectral order sinc-collocation method using a nonlinear fractional differential equation. The predictor–corrector method offers greater numerical stability and high accuracy thus making it suitable for complex fluidic dynamic systems.

We used a deterministic mathematical model [37] to quantify the dynamics of the immune response by developing first-order differential equations for each of the immune variables. The procedure of setting up the equations for the immune variables is similar to what has been used earlier in other studies [38], i.e., the rate of change in a variable is expressed as a linear combination of a minimum of two or more immunological processes that include basal production rate (π_nn), proliferation rate (α_nn), transformation/differentiation rate (θ_nn), suppression rate (δ_nn), uptake rate (µ_nn), and turnover rate (λ_nn). The subscript nn in the immune response process denotes the variable name.

Although the immune response to Ag stimulation in vivo is a result of coordinated interactions of the hundreds of variables, here, we have limited the evaluation to the major immune components known as the primary drivers of the immune response. We selected 12 immunological variables whose time-dependent evolution is determined by various immunologic parameters. The evolution of dynamic interactions was stretched to 21 days because most of the immune responses occur within this time frame and nearly all the variables display their optimal changes, growth, saturation, or steady-state behaviors within 21 days. Since we used in vitro antibody synthesis by PCs and I_g production profiles in response to PWM stimulation, the framework of our mathematical model is based on polyclonal activation of the T-cell-dependent I_g synthesis. In order to describe the behavior of B-cells differentiating into antibody-producing PCs, we considered various cells including T_h-cells, T_reg cells, B-cells, and APC in the presence of IL-2, IL-4, and IL-6 produced by T-cells. While T_h provides helper activity for B-cell activation, proliferation, and differentiation, T_reg and suppressor macrophages (M2) can inhibit the B-cell activation process or block T-cell helper activity. Interactions of these cellular components and their proliferative and suppressive functions are depicted in Figure 1.

2.3. Detailed Description of the Model Components

Antigen (Ag). Ag (PWM) triggers adaptive cellular and humoral immune responses, leading to the development of Ag-specific T_h cells, T_regs cells, and B-cell responses. PWM triggers a T-cell-dependent antibody response by B-cells. The rate of change in Ag due to its interaction [39] with B-cells and dendritic cells is described by

$$\frac{dAg}{dt}={\pi }_{Ag}-\left({\mu }_{AgB}B+{\mu }_{AgN}{DC}_N+{\mu }_{AgA}{DC}_A\right)Ag-{\lambda }_{Ag}Ag.$$

(1)

Here, Ag denotes the concentration of antigen in µg/mL (used in in vitro experiments), the constant ${\pi }_{Ag}$ represents basal production rates, the terms ${\mu }_{AgB}BAg$, ${\mu }_{AgN}{DC}_{N}Ag$, and ${\mu }_{AgA}{DC}_{A}Ag$ represent uptake of Ags by B-cells, naïve dendritic cells (DC_N) and activated dendritic cells (DC_A), respectively. The parameters, ${\mu }_{AgB},$ ${\mu }_{AgN}, \mathrm{a}\mathrm{n}\mathrm{d} {\mu }_{AgA}$ are the rates of elimination of Ag (Ag processing) by B-cells, naïve dendritic cells (DC_N) and activated dendritic cells (DC_A), respectively. The natural degradation of Ag is represented by ${\lambda }_{Ag}Ag$. The interaction of Ag with APC (B-cells, DC_N, DC_A) initiates the activation of the immune system that propagates through the immune network, leading to the elimination of antigens.

Naïve Dendritic Cells (DC_N). The main role of DC_N is to recognize, capture, and process Ags. After the processing of the Ag sequences into small pieces called peptides, they are presented to the TCR on T_h-cells in the context of class II MHC molecules on the cell surface. The interactions of Ag with DCs and their activation are expressed by the following equation [40]:

$$\frac{d{DC}_N}{dt}={\pi }_{DCN}-{\theta }_{DCA}{DC}_{N}Ag-{\theta }_{M2}{DC}_{N}Ag-{\lambda }_{DCN}{DC}_N.$$

(2)

where in the first term, ${\pi }_{DCN}$, is the basal production rate of DC_N; the second term, ${\theta }_{DCA}{DC}_{N}Ag$, is the rate of transformation of DC_N into activated dendritic cells (DC_A); ${\theta }_{M2}{DC}_{N}Ag$ represents the rate of transformation of DC_N cells/monocytes into M2 macrophages; and the last term,${\lambda }_{DCN}{DC}_N,$ is natural turnover rate of DC_N. It should be noted that high Ag doses can induce tolerance or suppression in both in vitro and in vivo models [41].

Activated Dendritic Cells (DC_A). Dendritic cells are activated by inflammatory cytokines to express a costimulatory molecule (CD80, CD86), which is recognized by specific receptors (CD28) on the T-cells. DC_A is required to trigger T-cell activation, differentiation into effector T-cells, and clonal expansion of naïve T-cells. We determine the dynamics of DC_A using Equation (3) below. The first term of the equation,${\theta }_{DCA}{DC}_{N}Ag,$ is the rate of generation of activated dendritic cells [42], which is equal to the rate of activation of the naïve DC into DC_A; and the last term,${\lambda }_{DCA}{DC}_A,$ is the degradation rate of the DC_A.

$$\frac{d{DC}_A}{dt}={\theta }_{DCA}{DC}_{N}Ag-{\lambda }_{DCA}{DC}_A.$$

(3)

Suppressor Macrophages (M2). Macrophages activated via an alternative pathway are referred to as suppressive macrophages or M2. Under non-infectious conditions, the presence of IL-4, IL-10, and TGF-β favors the polarization of monocytes or DC_N [43] to an M2 subpopulation that produces IL-10 in high concentrations along with TGF-β that are able to suppress immune effector cells and antibody synthesis. Equation (4) describes the rate of production of M2, which is equivalent to the differentiation of naïve dendritic cells into M2 via antigens interaction and their decay, ${\lambda }_{M2}{M}_2.$

$$\frac{d{M}_2}{dt}={\theta }_{M2}{DC}_{N}Ag-{\lambda }_{M2}{M}_2,$$

(4)

M2 macrophages suppress the T-cell activation process leading to suppression of antibody synthesis.

Naïve T-cell (T_n). Naïve T-cells are derived from cell division and thymic export. T_n cells circulate between the blood and peripheral lymphoid organs (lymph nodes, spleen, and mucosal lymphoid tissues) until they encounter a foreign Ag. This encounter activates T_n cells and initiates the differentiation of Tn into T_h and T_reg cells. Since our focus is on antibody responses, the development of cytotoxic T lymphocytes was excluded from the model. We have considered here that T_n cells [44] are produced at a constant rate ${\pi }_{Tn}$ and their dynamics in the periphery are determined by homeostatic mechanisms driven by cell death, differentiation, and rate of generation [44,45]. These processes are summarized in the following expression Equation (8).

$$\frac{d{T}_n}{dt}={\pi }_{Tn}-\left({\theta }_{Th}+{\theta }_{Treg}\right)T_n\frac{{DC}_A}{{DC}_A+{DC}_{A}50}-{\delta }_{Tn}{T}_n{M}_2-{\lambda }_{Tn}{T}_n,$$

(5)

The ${\theta }_{Th}{T}_n\frac{{DC}_A}{{DC}_A+{DC}_{A}50}$ and ${\theta }_{Treg}{T}_n\frac{{DC}_A}{{DC}_A+{DC}_{A}50}$ are transformation rates of T_n cells into T_h-cells and T_reg-cells, respectively. Since DCs concentration floats between 1% to 3% of the total cell population, we set ${DC}_{A}50=0.5\ast DCA0\frac{cells}{ml}~2800 \mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\mathrm{s}/\mathrm{m}\mathrm{L}$. DCA0 varies from 0.3% to 0.6 % of the total PBMC. The ${\delta }_{Tn}{T}_n{M}_2$ denotes the suppression of T_n cells by M2 macrophages, and the last term, ${\lambda }_{Tn}{T}_n$, is the degradation rate of T_n cells due to all possible mechanisms.

Helper T-cells (T_h). The dynamic behavior of T_h cells was evaluated in the stimulation of the cellular arm of immune responses. CD4+ T helper lymphocytes play a key role in the adaptive immune system exerting a wide spectrum of biological functions required for effective control of infections and avoidance of autoimmune diseases. Tight regulation of immune responses by CD4+ T-cells is exerted by various subpopulations of CD4+ T-cells including classical Th₁ cells, Th₂ cells, and Treg cells predominantly governed by the cytokines in the milieu and, the intensity of stimulation. When T_h cells are activated by APCs, they not only secrete multiple cytokines and chemokines but also express specific stimulatory co-receptors on their surface. The behavior of T_h cells is given by

$$\frac{d{T}_h}{dt}={\theta }_{Th}{T}_n\frac{{DC}_A}{{DC}_A+{DC}_{A}50}+{\alpha }_{Th}{T}_h\frac{IL2}{IL2+IL250}\frac{Ag}{Ag+Ag50}-{\theta }_{Tr}{T}_{reg}{T}_h-{\lambda }_{Th}{T}_h.$$

(6)

The first term, ${\theta }_{Th}{T}_n\frac{{DC}_A}{{DC}_A+{DC}_{A}50}$, represents a generation of T_h from naïve T-cells [46,47,48]; the second term, ${\alpha }_{Th}{T}_h\frac{IL2}{IL2+IL250}\frac{Ag}{Ag+Ag50}$ represents a proliferation of antigen-bounded T_h by IL-2; ${\theta }_{Tr}{T}_{reg}{T}_h$ represents the depletion of T_h cells due to their differentiating into T_reg cells; and the last term is their turnover rate. In our calculations, we have set $\mathrm{I}\mathrm{L}250=5.0 \mathrm{n}\mathrm{g}/\mathrm{m}\mathrm{L}.$

Regulatory T-cells (T_reg). The primary function of T_reg cells is to prevent autoimmune diseases by maintaining self-tolerance. However, the main suppressive function of T_reg cells in the context of antibody production is exerted by immunosuppressive cytokine TGF-β, which can inhibit the secretion of immunoglobulins [49]. Treg cells have essential roles in the maintenance of peripheral tolerance, immune homeostasis, and prevention of autoimmunity by regulating effector T-cell responses and thus preventing their potentially pathogenic effects via various mechanisms [50]. Features such as the plasticity of T_reg cells have changed the understanding of T_reg-cell biology in terms of their interaction with other immune and non-immune cells, their functions in specific tissues, and the implications of this for the pathogenesis of autoimmune diseases. In this paper, we have restricted the role of T_reg cells immune cells, and antigens. We determine the dynamics of T_reg cells [51] using the following equation.

$$\frac{d{T}_{reg}}{dt}={\theta }_{Treg}{T}_n\frac{{DC}_A}{{DC}_A+{DC}_{A}50}+{\alpha }_{Treg}{T}_{reg}\frac{IL2}{IL2+IL250}\frac{Ag}{Ag+Ag50}+{\theta }_{Tr}{T}_{reg}{T}_h-{\lambda }_{Treg}{T}_{reg},$$

(7)

The first term, ${\theta }_{Treg}{T}_n\frac{{DC}_A}{{DC}_A+{DC}_{A}50}$, represents the generation of T_reg cells from the naïve T-cells; ${\alpha }_{Treg}{T}_{reg}\frac{IL2}{IL2+IL250}\frac{Ag}{Ag+Ag50}$ is the proliferation of T_reg by the IL-2, ${\theta }_{Tr}{T}_{reg}{T}_h$ represents rate of generation of T_reg cells due to differentiation of T_h cells; and the last term is their natural decay rate, ${\lambda }_{Treg}{T}_{reg}.$

B-cells. B-cells when activated differentiate into PCs and secrete antibodies [7]. While some Ags can trigger a direct B-cell response, most B-cell responses are derived after interaction with the T_h cells. Activated T_h cells produce IL-4 which helps in the proliferation of B-cells and IL-6, which supports differentiation of B-cells. B-cell behavior is determined by the equation below, which incorporates these processes.

$$\frac{dB}{dt}={\pi }_B+{\alpha }_{BTh}B{T}_h+{\alpha }_{IL6}\frac{IL6}{IL6+IL650}B-{\theta }_{BTh}B{DC}_A-{\lambda }_{B}B,$$

(8)

${\pi }_B$ is the basal level generation rate of B-cells; and ${\alpha }_{BP}B{T}_h$ represents a positive feedback by T_h via interaction with B-cells. Cytokine proliferation of B-cells is given by ${\alpha }_{IL6}\frac{IL6}{IL6+IL650}B$ term. For simplicity, we have taken IL650 = $\mathrm{I}\mathrm{L}250$. The term ${\theta }_{BTh}B{DC}_A$ represents the differentiation of B-cells into PCs mediated by T_h cells and ${\lambda }_{B}B$ represents B-cell turnover rates.

Plasma cells (PCs). The long-lived PCs represent the terminal differentiation step of B-cells that secrete antibodies (Abs). The dynamics of PC [52,53] is given by

$$\frac{dPC}{dt}={\theta }_{BTh}B{DC}_A-{\lambda }_{PC}PC.$$

(9)

T-cell-dependent antibody responses require the activation of B-cells by T_h cells that respond to the same antigen. The first term, ${\theta }_{BTh}B{DC}_A$, describes a generation of PCs due to B-cells differentiation by activated dendritic cells and ${\delta }_{PC}PC$ describes their natural turnover rates.

Immunoglobulins (Ig). Ig are glycoproteins produced by B-cells that play a crucial role in protective immunity. Structurally, Ig contains two components: the Fc region (tail part) and the Fab region; the Fc tail of the Ig is responsible for triggering effector functions of the innate immune system while the Fab arms are responsible for antigen binding. Among the subtypes of Ig, the immunoglobulin G (IgG) class of antibodies is a major effector molecule of the humoral immune response in humans and accounts for about 75% of the total immunoglobulins in plasma. Antibodies of the IgG class have a relatively high antigen affinity and long serum persistence. Antibody production via PCs is described below as follows:

$$\frac{dIg}{dt}={\alpha }_{Ig}PC-{\lambda }_{Ig}Ig.$$

(10)

The first term, ${\alpha }_{Ig}PC$, is the rate of antibody production by PCs, and the last term, ${\lambda }_{Ig}Ig$, is Ab degradation rate.

Cytokine IL-2. IL-2 is a T-cell growth factor that supports the expansion of T-cells after engaging with Ag. IL-2 is essential for the development of T-cells and T_reg. T_regs function to prevent other T-cells from recognizing and reacting against “self-Ags”. Not only do T-cells produce IL-2 [54] but they consume IL-2 for their proliferation [12]. The dynamics of IL-2 is given by

$$\frac{dIL2}{dt}=\left({\alpha }_{IL1}+{\alpha }_{IL2}\frac{{DC}_A}{{DC}_A+{DC}_{A}50}\right)T_{h}f\left(t\right)-\left({\mu }_{ThIL2}{T}_h+{\mu }_{TrIL2}{T}_{reg}\right)IL2-{\lambda }_{IL2}IL2.$$

(11)

The first term, $\left({\alpha }_{IL1}+{\alpha }_{IL2}\frac{{DC}_A}{{DC}_A+{DC}_{A}50}\right)T_{h}f\left(t\right)$ describes the production of IL-2 by T_h [46,55]; $\left({\mu }_{ThIL2}{T}_h+{\mu }_{TrIL2}{T}_{reg}\right)IL2$ term accounts for IL-2 consumption by T_h and T_reg; and the last term, ${\lambda }_{IL2 }IL2$, accounts for their normal degradation. The dynamic offset between two opposing forces (proliferation and consumption terms of IL-2) determines the trajectory of the immune response. IL-2 production is transient due to a network of autocrine and paracrine signals present in normal and pathological responses. We have introduced a factor $f\left(t\right)=\left(\frac{t}{{\tau }_0}\right).\exp \left(-{\left(\frac{t}{{\tau }_0}\right)}^2\right)$ to account for this transient behavior by setting τ₀ to 5 days.

Cytokine IL-4 and IL-6. IL-4 and IL-6 are critical in the activation process. While IL-6 helps in the differentiation, growth, and activation of T and B lymphocytes and stimulates the production of antibodies, IL-4 induces the differentiation of naive helper T-cells [56,57]. IL-4 and IL-6 are essential to transition from innate to acquired immunity and drive B-cell proliferation and differentiation, respectively. In this model, for brevity, we have combined [58,59] the effects of both IL-4 and IL-6 into just one variable, IL-6. The rate of change in IL-6 is given by

$$\frac{dIL6}{dt}=\left({\alpha }_{IL1}+{\alpha }_{IL2}\frac{{DC}_A}{{DC}_A+{DC}_{A}50}\right)T_{h}f\left(t\right)-{\mu }_{BIL6}BIL6-{\lambda }_{IL6}IL6.$$

(12)

The second term, ${\mu }_{BIL6}B.IL6,$ is the rate of consumption of IL-6 by B-cells for their differentiation and the natural degradations of IL-6 is given by ${\lambda }_{IL6}IL6$.

2.4. Parameters Determination

The profiles of the immune response variables are strongly dependent on the parameter values used in the equations and the initial (time = 0) values of the variables are shown in

Table 1. Initial values of the variables.

No	Variables	Values (per mL)
1	Ag (Ag)	0.5 µg, 1.0 µg and 2.0 µg
2	Naive Dendritic cells (DCN)	16,000
3	Activated Dendritic cells (DCA)	8000
4	Regulatory Macrophages (M2)	0
5	Naïve T-cells (Tn)	520,000
6	B-cells (B)	20,000
7	Antibody, Ig	0
8	Plasma cells (PC)	0
9	Treg cells	(1%, 5%, and 10%) of Th
10	T Helper cells (Th)	150,000
11	Interleukin 2 (IL-2)	10 ng
12	Interleukin 4 (IL-4)	10 ng

. The parameter and basal values are strongly dependent on each other. For a normal healthy person, the values of the variables lie in a certain range as shown in Figure 2. The range dictates the cell populations and their distributions in 1 mL of PBMC (1 × 10⁶ cells) as we move from left to right in the figure. Based on these values, we have taken basal values of the variables in 1 mL of PBMC as shown in

Table 2. Description of parameters, mean, and range of values.

Parameter (Units)	Description	Mean Value	Min Value	Max Value
${\pi }_{Ag}$ (ng/day−1)	Basal production rate of Ag	1.621 × 10−1	1.093 × 10−1	1.990 × 10−1
${\mu }_{AgB}$ (1/cell.day)	Rate of Ag binding with B-cells	9.03 × 10−7	2.18 × 10−7	1.69 × 10−6
${\mu }_{AgN}$ (1/cell.day)	Rate of Ag binding with naïve dendritic cell	8.27 × 10−7	1.06 × 10−7	1.61 × 10−6
${\mu }_{AgA}$ (1/cell.day)	Rate of Ag binding with activated dendritic cells	8.74 × 10−7	2.09 × 10−7	1.67 × 10−6
${\lambda }_{Ag}$ (day−1)	Turnover rates of antigen	2.50 × 10−1	1.01 × 10−1	4.94 × 10−1
${\pi }_{DCN}$ (cells/mL.day)	Basal production rate of dendritic cells	4.54 × 102	2.04 × 101	1.00 × 103
${\theta }_{DCA}$ (1/ng.day)	Rate of DC binding with Ags, which leads to their differentiation into activated DC	2.53 × 10−2	4.29 × 10−4	5.22 × 10−2
${\theta }_{M2}$ (1/ng.day)	Rate of DC binding with Ag, which leads to their differentiation into M2 cells	1.85 × 10−2	3.20 × 10−3	3.44 × 10−2
${\lambda }_{DCN}$ (1/day)	Decay rates of DC	1.77 × 10−1	9.49 × 10−2	2.98 × 10−1
${\lambda }_{DCA} (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Decay rates of DCA	3.90 × 10−1	2.05 × 10−1	5.94 × 10−1
${\lambda }_{M2} (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Decay rates of M2	3.75 × 10−1	3.03 × 10−1	4.86 × 10−1
${\pi }_{Tn}$ (cells/day)	Basal production rates of naïve T-cells	9.03 × 100	8.11 × 100	9.90 × 100
${\theta }_{T\mathrm{h}} (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Rate of Tn cells differentiating into Th cells	2.61 × 10−5	9.65 × 10−6	4.99 × 10−5
${\theta }_{Treg} (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Rate of Tn cells differentiating into Treg cells	2.95 × 10−4	5.12 × 10−5	4.78 × 10−4
${\delta }_{Tn}$ (1/cell.day)	Rate of suppression of Tn cells by M2	5.49 × 10−6	1.03 × 10−6	9.35 × 10−6
${\lambda }_{Tn} (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Turnover rate of Tn cells	3.27 × 10−2	1.07 × 10−2	5.56 × 10−2
${\pi }_B$ (cells/day)	Basal rate of production of B-cells	3.48 × 102	3.03 × 102	3.91 × 102
${\alpha }_{BT\mathrm{h}}$ (1/cell.day)	Rate of proliferation of B-cells by Th cells	3.25 × 10−7	5.55 × 10−7	4.42 × 10−7
${\alpha }_{IL6} (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Growth rate of B-cells by IL-6	3.81 × 10−4	2.83 × 10−4	4.63 × 10−4
${\theta }_{BT\mathrm{h}}$ (1/cell.day)	Rate differentiation of B-cells into B2-cells induced by Th cells	1.48 × 10−10	9.07 × 10−12	2.83 × 10−10
${\lambda }_B (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Degradation rates of B-cells	2.05 × 10−2	1.62 × 10−2	2.49 × 10−2
${\lambda }_{PC} (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Degradation rates of PC- cells	1.64 × 10−2	4.00 × 10−3	2.63 × 10−2
${\theta }_{Tr}$ (1/cell.day)	Rate of Th cells differentiation into Treg cells	1.20 × 10−6	6.18 × 10−7	2.13 × 10−6
${\alpha }_{T\mathrm{h}} (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Proliferation rate of Th cells by IL2	2.70 × 10−1	1.10 × 10−1	3.47 × 10−1
${\lambda }_{T\mathrm{h}} (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Turnover rate of Th cells	1.39 × 10−1	3.86 × 10−2	2.96 × 10−1
${\alpha }_{Treg} (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Proliferation rate of Treg by IL2	1.19 × 10−4	8.76 × 10−5	1.57 × 10−4
${\lambda }_{Treg} (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Turnover rates of Treg cells	2.13 × 10−1	9.95 × 10−2	2.94 × 10−1
${\alpha }_{IL1}$ (ng/cell.day)	Expansion rates of IL-2 by Th cells	8.67 × 10−1	1.21 × 10−1	1.59 × 100
${\alpha }_{IL2}$ (ng/cell.day)	Expansion rates of IL-2 by activated Th cells	1.01 × 100	2.42 × 10−1	1.46 × 100
${\mu }_{T\mathrm{h}IL2}$ (1/cell.day)	Rate of consumption of IL-2 by Th cells	2.80 × 10−2	6.75 × 10−3	5.98 × 10−2
${\mu }_{TrIL2}$ (1/cell.day)	Rate of consumption of IL-2 by Treg cells	2.35 × 10−2	1.39 × 10−3	4.98 × 10−2
${\lambda }_{IL2} (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Rate of degradation of IL-2	3.12 × 10−2	1.05 × 10−2	5.00 × 10−2
${\mu }_{BIL6}$ (1/cell.day)	Rate of consumption of IL-6 by B-cells	8.18 × 10−2	4.35 × 10−3	1.54 × 10−1
${\lambda }_{IL6} (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Rate of degradation of IL-6	2.79 × 10−2	1.13 × 10−2	4.94 × 10−2
${\alpha }_{Ig}$ (µg/cell.day)	Rate of production of Ig by plasma cells	3.44 × 10−1	1.77 × 10−1	4.85 × 10−1
${\lambda }_{Ig} (1/\mathrm{d}\mathrm{a}\mathrm{y})$	Turnover rate of Ig	2.61 × 10−2	2.32 × 10−2	2.95 × 10−2

. We solve the system of differential Equations (1)–(12) using the R programming language. For the integration of these equations, we used the ODE function of the package “deSolve”. Using Equations (1)–(12), we determined parameter values that generate the normal immune response of a healthy person. The computational algorithm follows the following steps: we start with the basal values of the variables and use a uniform distribution function bounded by minimum and maximum values to generate very large values of the parameter for each parameter. Starting with the first value of each parameter, the time-dependent equations are solved for the time window of 21 days retaining only the solutions that are convergent. This process is repeated for all the values of the parameters and thus generates a very large set of output solutions. We analyzed the output profiles of the immune variables and selected only those parameter sets that produced the normal immune response of the cells as shown in Figure 3 where each cell and cytokine show characteristic features in its profile. To obtain a large collection of the time-dependent profiles of the variables, we varied the minimum and maximum values of the parameters and generated a new set of parameters, and repeated the calculations. In the end, we collected only those parameters that produced the desired profiles of the immune variables. From this collection of the parameters, we determined their mean, minimum, and maximum values as reported in Table 2.

3. Results

Overview of the Mathematical Model. A theoretical model was developed to assess the immune response to PWM by PBMC from healthy donors incorporating cell–cell interactions of APC-T-cells, T-cells-B-cells, B-cells-APC, and contributions from cytokines that promote proliferation and differentiation of T and B-cells. The system of equations was solved for a set of parameter values presented in Table 2 using the initial values for the variables, in Table 1, a range of normal immune responses was determined for the immune response variables—Ag, naïve (DCN) and activated DC (DCA), naïve (T_n), helper (T_h) and regulatory (T_reg) T-cells, B-cells, plasma cells (PC), M2, I_g, IL-2, IL-4 and IL-6, and I_g synthesis over a period of three weeks. The three-week interval was selected to capture the time course in which the immune responses and induction-specific antibodies are essentially complete when an individual is vaccinated or challenged with a pathogen.

In Vitro Experimental Data Supporting the Model. The in vitro model was used as a guide for establishing the range of immune response variables [18,60,61]. Individual I_gG-secreting B-cells were measured by B-cell ELISpots in cultures of PBMC, CD4⁺ T-cells + non-T-cells (B-cells, monocytes, and NK cells), or CD8⁺ T-cells + non-T-cells in co-cultures stimulated with PWM (Figure 2B). The data from individuals published by Lum et al. [18,60,61] were reanalyzed and presented in the context of this study. Similarly, purified CD4⁺ or CD8⁺ T-cell subsets were evaluated for enhancement or suppression of control culture containing a fixed number of T + B-cells (Figure 2C). The effect of monocytes and subpopulations of monocytes was estimated based on the effects of adding purified monocytes to purified T and B-cells (monocytes were depleted) in the same assay using normal and patienT-cells after HLA-identical allogeneic sibling transplants [10].

Recapitulating Experimental Data in Mathematical Model. The model’s results are based on experimental data we used in our equations via the initial time = 0 values of the variables (Table 1). In Figure 3, we show the simulation profiles for antigen concentration, IL-2, IL-4/IL-6, naïve DC, activated DC, M2 macrophages, naïve T-cells, Th₁ cells, Treg cells, B-cells, PC, and I_g concentration over a period of three weeks using in vitro experimental results of in vitro polyclonal antibody synthesis PWM-stimulation.

3.1. Dynamics of the Variables

Dose responses to PWM. The dynamic and kinetics of the immune responses to three different doses (0.5 µg, 1.0 µg, and 2.0 µg/mL) of Ag were evaluated. As a function of time, the Ag dose declines from the initial value due to its binding with APCs and the natural degradation process (Figure 3a). Increasing initial values of Ag increases the overall profile of Ag but its decay dynamics remain the same.

Dynamics of Dendritic Cells. DC_N decline with time as the DC_N differentiate into activated DCs (Figure 3d), while the DC_A population initially increases due to a higher rate of a generation that peaked at days 1–2 but declines later due to turnover rates (Figure 3e); M2 macrophages increase initially but later they begin to decline in numbers due to the dominance of turnover rate to their generation rate (Figure 3f). This behavior is consistent with the recent work by Lichtnekert and Kawakami [62].

Dynamics of T-cells. Naïve T-cells decline continuously from their initial value soon after activation and the beginning of differentiation into effector T-cells. This is due to their higher rates of differentiation into T_h and T_reg cells after stimulation and activation (Figure 3g). Soon after the presentation of Ag, the number of T_h cells rose from the basal state as expected in the normal immune response. T_h cells continue to rise due to additional positive feedback driven by IL-2 produced by T_h cells themselves until the number of T_h cells reaches a peak between 5 to 15 days. Later, T_h cells start declining due to a decrease in IL-2 levels as well as due to their natural turnover rate (Figure 3h) and eventually reach a homeostasis value. The proliferation of the T_h cells and T_reg cells is determined by the IL-2, Ag, and activated DC whose functional potency declines after a week from the start of immune activation. The time-dependent profile of T_reg cells is a reciprocal of the T_h cells pattern. The number of T_regs is small during the initial expansion phase of T_h cells, and once T_reg starts increasing (Figure 3i), the T_h cells begin to decline due to suppressor activity mediated by T_reg cells.

Dynamics of B-cells. B-cells show a slower increase due to their differentiation into PC and turnover rate (Figure 3j). Around 3 weeks after stimulation, the number of B-cells peaks and begins to decline whereas PC increases, peaks, and declines (Figure 3k). Antibody synthesis increases steadily during the timeframe (Figure 3l) and declines around the 4th week after stimulation.

Dynamics of Cytokines. The IL-2 (Figure 3b) and IL-6 patterns are functionally similar for T and B-cells, respectively. The main source of generation for IL-2 and IL-6 are T_h cells. The IL-6 has similar kinetics except it is being utilized by B-cells (Figure 3c).

3.2. Effects of Varying Ag Dose on the Dynamics of Immune Response

In addition to the input parameter values (Table 2) required for evaluation of the coupled differential equations, the outcomes of the equations also depend on the initial (time = 0 day) values of the variables as given in Table 1. Since Ag is the only external element that induces immune responses, the effects of varying Ag doses were investigated. Given a highly non-linear relationship between the immune response components, it is not trivial to predict the behavior of each element to the changing Ag doses. Figure 4A shows the variations in the immune profiles for all the 12 variables at 0.5 µg/mL, 1.0 µg/mL, and 2.0 µg/mL values of Ag. Apart from T_n, IL-2, and M2, all variables exhibited pronounced increases in their concentrations on increasing Ag loads, whereas DC_N decreased with increasing Ag loads. The reciprocal response pattern of DC_N can be understood from Equation (2) where Ag drives the differentiation of DC_N into DC_A and M2 leading to decreases in DC_N with increasing Ag. This explains why M2 does not change with varying doses of Ag see Equation (4). Naive T-cell profiles are nearly independent of Ag concentration because of M2 behavior and saturation effects of DC_A at their higher concentrations. Increasing doses of Ag increases the overall concentration profile of Ag but its decaying behavior remains the same. Since the kinetics of IL-2 is determined by the time-dependent behavior of T_h and T_reg, changes in absolute numbers of T_h and T_reg cells in the presence of Ag are compensated by rates of IL-2 production and consumption rates. For this reason, the IL-2 profile does not change much with Ag. There is a steep increase in the slope of T_h at higher Ag load (2.0 g/mL) before it peaks. The behavior of T_reg cells as a function of increasing Ag is similar to T_h cells, however, T_reg showed a much stronger increase at its peak. The B-cell population also increases with Ag dose but only at later times after activation. Increasing B-cell population increases the PC profile and, in turn, increases I_g production, which continues to increase even beyond the time frame of 3 weeks shown here. Ag dose clearly affects the peak height and peak position of immune cell profiles.

The changes in the values of the peak heights of some of the selected variables—T_h, T_reg, B-cell, PC, and I_g—are shown in Figure 4B (left panel). All these variables show a linear increase in the peak heights with increasing Ag dose. The concentration of T_h cells at the peak increases by a factor of 1.5 while that of T_reg at the peak increases by a factor of 2.3. The B-cell and plasma cell numbers increase by a factor of 2 and 2.5, respectively. Figure 4B right panel shows shifts in the profile patterns of T_h, T_reg, and PC with increasing Ag loads. Other variables did not show changes in their peak patterns and hence are not shown here. The most remarkable increase occurred in the T_reg pattern that increased to day 11 from day 4 with increasing Ag dose, while the T_h peak position shifted to day 6 from day 4. PC peak position showed a gradual increase with increasing Ag load.

3.3. Effects of Varying Initial Value of T_regs on the Dynamics of Immune Response

Given that T_regs in healthy subjects vary from 1 to 10% of T_h, simulations were conducted at 1%, 5%, and 10% of T_h (Figure 5A). T_h, T_reg, B, PC, and I_g were the components that apparently changed with varying initial values of T_regs. The shifts for each of the above-mentioned components were statistically different by <0.05. Except for T_regs, all other variables decreased by ~10%. The most remarkable was a 5× increase in T_reg cells while T_h showed a ~16% decline when the initial value of T_reg was increased from 1 to 10%. Similarly, B-cells, PC, and I_g showed a decline in their overall profiles when T_reg initial values increased to 10%. Changes in peak heights of these variables are shown in Figure 5B (left panel). Unlike the T_reg population, B-cells, PC, T_h, and I_g decrease in magnitude. This behavior is consistent with a healthy immune system response. B-cells declined by 32% while PC declined by 22% when T_regs increased from 1 to 10%. There was a concomitant decline in I_g synthesis while the PC numbers declined. In Figure 5B (right panel), we show the shifts in the peak positions of immune variables. T_h, T_reg, and PC showed decreases in the peak position with increasing T_reg initial numbers. The largest shift is in the T_reg maximum peak position, which shifted to day 6 from day 10 on increasing proportions of T_regs from 1% to 10%.

4. Discussion

This study reports a mathematical model for normal T-cell-dependent B-cell antibody synthesis using PWM. This system reflects changes in tightly regulated responses by helper and regulatory T-cells that could be used to simulate dysregulated immune responses [6,63]. The kinetics of the temporal interactions for T_h, T_regs, B, and dendritic cells were incorporated in the model for optimum production of PWM-stimulated antibody production, which is consistent with in vitro experimental data [13,16,19,64,65]. The key question for this model was “what are best proportions of T-cells, T-cell subsets, B-cells, and APC required to generate optimal antibody synthesis in a 3-week time frame? To address this question, various conditions were simulated using different parameter values to optimize the differentiation of B-cells to I_g-producing PC (Figure 3). The peak proliferative response of T_h was on day 6, which was associated with a concurrent rise in PC numbers that peaked on day 8. Unlike the rapid increase in PC, the increase in the accumulated I_g production continued up to day 21 (Figure 3k,l). Data in Figure 3 were consistent with our previous in vitro studies [13,16,19,64,65] that were considered optimal for all immune parameters.

Increasing or decreasing the Ag dose altered I_g production when compared to a standard dose of PWM (1.0 µg/mL) shown in Figure 3. As the Ag dose is increased, both qualitative and quantitative changes occur in immune response profiles and their peak responses (Figure 4A,B). The number of T_h cells increased by 32% and the number of B-cells increased by 54%, which led to a ~70% increase in PC. As a result of the increase in PCs, there is a corresponding increase in I_g production, which continues even after three weeks as seen in Figure 4A. Similar to T_h cells, T_reg cells also increased by 90%. Interestingly, T_reg cells did not negatively impact I_g synthesis at higher numbers.

Decreasing Ag load by half (0.5 µg/mL) leads to a 20–40% decrease in the profiles of T_h cells, PC, and I_g production as well as a reduction in B-cells differentiating into PCs. On the other hand, T_regs decreased by 20% at their peak value on day 12 (Figure 4A,B). The weaker stimulatory effect of reduced Ag dose on T_h and B-cells resulted in reduced differentiation of B-cells to PC and thus reduced I_g production regardless of T_reg numbers. The features of other cells and IL-2 remain almost the same on reducing the Ag dose. Our model shows that Ag concentration levels have noteworthy effects on the immune cell profiles.

The effects of increasing or decreasing the proportion of Tregs on I_g production are shown in Figure 5A,B. Similar to the simulation results of varying the Ag dose, changing the proportion of T_reg initial value qualitatively and quantitatively affects profiles of immune response. As the value of T_regs increased from 5% to 10% (2×), there was a ~10% decrease in T_h cells accompanied by a slight decrease in their peak position while the peak value of T_regs increased by 60% compared to standard T_regs numbers at 5%. The number of PCs diminished by 10%, resulting in a corresponding decrease in I_g produced.

Increasing T_regs was associated with changes in the immune response profiles; however, decreasing T_regs to 1% showed a marked decrease by 70% in T_regs followed by a delay in the peak performance on day 10 at 1% compared to day 8 at 5%. Decreasing T_regs cells showed positive effects on T_h cells and B-cell numbers while decreasing to 1% showed much less of an effect on PC numbers and slightly enhanced I_g synthesis (Figure 5B). However, a number of antibodies producing PCs remain more or less unchanged, suggesting that decreasing or increasing T_regs may enhance or suppress antibody synthesis, respectively, by augmenting or suppressing T_h cells so that B-cells and PCs receive less help to drive B-cell proliferation and differentiation.

The development of an in vitro model that simulates an in vivo systemic immune response to the “vaccination” of normal subjects or patients provides a model for dissecting immunopathologic mechanisms, as well as helps estimate and/or predict the effects of immunomodulatory drugs on patients ongoing immunotherapy [66]. In this model, exposure to pathogenic Ags or tumor-associated antigens would trigger network responses via cell–cell interactions leading to a concerted immune response to secrete antibodies and cytokines [4,67]. This model shows how network interactions between Ag dose, type, and dose of various cell types, and cytokines are dependent on initial or basal states of these variables and the strength or potency of the immune component. The model permits exploration of how manipulating specific or multiple components can provide insight into the type and magnitudes of immune responses or outcomes. Over the past decade, mathematical models describing T-cell activation [45], specific immune response to infections [68], cancer cell growth [69,70], and autoimmune diseases [71,72] were successful in explaining experimental observations, providing new insights, and making predictions of complex immunological processes [40,73]. In addition to predicting immune response for various triggers, another major application of mathematical modeling has been on the dose estimation and schedule to deliver immunotherapeutic or anti-cancer agents. Recently, Carla et al. [66] developed a mathematical model to determine the minimum dose for the yellow fever vaccine to achieve optimal immunity. The parameter values of such models can provide insights into the clearance of the infection, maintenance of a chronic infection, or recurring infections [74,75,76]. Earlier, some mathematical models [77,78,79] have investigated the role of tumor-immune interaction and provided some insight into tumor progression, therapeutic dose estimation, therapeutic efficacy, and treatment regimens [77,78,79].

In summary, the model we presented here provides a real-time relationship among some of the mutually competing immune response variables in normal physiological conditions when challenged by some external antigen. It provides a predictive framework to address immune responses not only to preclinical models for immune therapeutics but also for clinical interventions such as vaccines for infectious diseases, immune-modulating agents in cancer, or immune deficiency disorders. In the future direction, we would like to investigate preventive modeling for autoimmune diseases and treatment response of infections and inflammation. We expect, in the future, with the advancement of new techniques in computer science and mathematics, mathematical modeling may provide novel insights into the cellular mechanisms of system biology.