The major key players and biological processes in the proposed MMQE (Figure 2) are: modeling the variations of the antigen (pathogen) agent, the joint interaction of cytokine agents IL-2 and IL-4, the dynamics of activation, proliferation and differentiation of B cell agent, the dynamics of activation and proliferation of T cell agent, the proliferation of plasma B cell agents and the activation of DC agent for overall immune system response. We also account for the life span of various lymphocytes during an immune response. For reasons of consistencies, we assign variables to agents involved in MMQE (Table 1).

One of the most effective and widespread strategies employed by pathogens in order to escape the human immune system is antigenic variation (Vink et al., 2012). This process can be described as a system that allows a pathogen to change its (surface) antigens that are presented to, and targeted by, the host’s adaptive immune system (Vink et al., 2012; Mosa, 2020). As a consequence of antigenic variation, subpopulations of antigenically distinct organisms arise within a population that are temporarily not recognized by the primary adaptive immune response (Barbour et al., 2006; Schmid-Siegert et al., 2017). This variability provides the pathogen with an extended window of opportunity to persist within a host or to infect a previously colonized host (Vink et al., 2012).

We define the dynamics of the pathogen agent distribution in the immune system as a random variable. Often, the distribution of activity per cell was found to be well described by a log normal distribution function (Koch, 1966) which is a type of probability distribution in which the logarithm of a variable is normally distributed (Crow and Shimizu, 1987; Balakrishnan and Chen, 1999). Correspondingly, theoretic modelling of cell survival as a function of mean activity per cell showed that survival curves differed substantially when the activity per cell was log normally distributed versus normally distributed (Limpert et al., 2001; Neti and Howell, 2006). Moreover, Sartwell’s model has been used to demonstrate that the incubation period for infectious diseases and genetic diseases fits a lognormal distribution (Horowitz et al., 2001). We therefore consider the variation of pathogens in this work as variations of the microscopic state of antigens ( A ag ) , described by the log normalized density function (Mayer et al., 2015; Carignano and Dalchau, 2018) as: ln ( A ag ) ∼ N ( μ ag , V ag ) (1) where μ ag and V ag stands for the mean and variance respectively of the variation.

We next illustrate with ODEs the complex interactions between the dynamics of cytokines IL-2 and IL-4. In general, cytokines are secreted proteins that carry intra and intercellular signals which regulate the immune response to pathogens or other pathological conditions (Chen et al., 2018; Muñoz-Carrillo et al., 2019). To effect their functions, cytokines must bind to specific extracellular receptors (Schreiber and Walter, 2010). Specifically, IL-2 and IL-4 are growth factors secreted by both T and B cells and demonstrate antagonistic marginal effects during immune response (Schreiber and Walter, 2010). The proliferation response to IL-2 stimulation is controlled by three factors: 1) the initial concentration of IL-2, 2) the density of IL-2R (receptor specific to IL-2) on the surface of the cells and 3) the duration of the cell exposure to IL-2 (Burke et al., 1997). On the other hand, studies on B cells have reported that IL-4 represents a growth and differentiation factor for the activated B cells and the IL-4 binding may have inhibitory effects on IL-2 induction (Castro et al., 1999). Therefore, a mathematical model based on the synthesis of all the known interactions is necessary to predict with high accuracy, the dynamics of T cell and B cell immune response. Let N and n be the bounded and unbounded number of cytokines such that N 2 and N 4 represent the bounded numbers of IL-2 and IL-4 to a specific cell with IL-2R and IL-4R receptors respectively. Likewise, n 2 and n 4 represent the unbounded numbers of IL-2 and IL-4 from IL-2R and IL-4R respectively. Based on the law of mass action, the quantitative terms for IL-2 and IL-4 binding can be derived and described by the following (Eqs 2–5). d n 4 dt = ϕ 4 − θ 4 I 4 n 4 + φ 4 N 4 − γ 4 n 4 (2) d N 4 dt = θ 4 I 4 n 4 − φ 4 N 4 − δ 4 N 4 (3) d n 2 dt = ϕ 2 − θ 2 I 2 n 2 + φ 2 N 2 − γ 2 n 2 − κ N 4 n 2 (4) d N 2 dt = θ 2 I 2 n 2 − φ 2 N 2 − δ 2 N 2 (5) where, I 2 and I 4 represent the concentration of IL-2 and IL-4 respectively. The parameters ϕ 4 and ϕ 2 represent the production rates of unbound receptors IL-4R and IL-2R (Burke et al., 1997). φ 4 N 4 performs the dissociation of IL-4 from IL-4R. The factor γ 4 n 4 serves as the internalization of unbounded receptor; the binding of IL-4 to unbounded receptor IL-4R is resumed by the term θ 4 I 4 n 4 . The terms in Eq. (5) are identical to the terms expressed in Eq. (3) which is modeled with respect to IL-2 and IL-2R. Likewise, except the additional term κ N 4 n 2 in Eq. (4), all the terms involved in Eq. (4) are identical to the terms in Eq. (2) with respect to IL-2 and IL-2R where κ N 4 n 2 models the phenomenon of down-regulation of unbounded IL-2R by IL4 binding. The parameter κ defines the downregulation of IL-2R by IL-4.

CD4+ T cells are a key part of the adaptive immune system. They differentiate into different phenotypes to carry out different functions. With the purpose of understanding the immune system as an integrated whole, Wertheim et al. present a new multi-approach and multi-scale modeling framework that can be used to model diverse immune responses at molecular, cellular, and systemic scales (Wertheim et al., 2021). In this study, cell population dynamics and systemic cytokine concentrations were described by ODEs (Wertheim et al., 2021). Excitingly, their model demonstrates the dynamics of CD4+ T cells in response to influenza infections in different cytokine milieus, considering heterogeneous populations of T helper cells and 11 cytokines in three spatial compartments (an infection site or target organ, lymphoid tissues, and a circulatory system). Furthermore, studies on the interplay among extracellular cytokines and the differentiation of T cells provide new insights into the plasticity of CD+ T cell differentiation and the generation of T cell phenotypes (Puniya et al., 2018). They model cellular dynamics in terms of concentration of cytokines which in turn control the differentiation of CD4+T cell. Similarly, MMQE uses a multiscale and multicellular model to assess immune response of T cells in terms of different activation protocols of the immune system in addition to the downregulation effect of IL-4 on IL-2 in the presence of antigens.

By default, T-cell antigen receptor (TCR) signaling is essential for activation, proliferation, and effector function of T cells. To this end, much attention has been focused on determining the minimum signaling requirements for TCR driven proliferation (Au-Yeung et al., 2017). Produced primarily by activated T cells, IL-2 enhances activation of T cells as well as promotes the differentiation of T cells into effector T cells and memory T cells (Liao et al., 2011). IL-2 serves dual opposing functions; it potently amplifies proliferative responses of effector T (Teff) and natural killer (NK) cells, while regulating immune homeostasis by driving regulatory T (Treg) cell proliferation, differentiation, and function (Abbas et al., 2018) and both axes have been leveraged to treat human diseases (Ghelani et al., 2020). Consequently, neither naive CD4 nor naive CD8 T cells are responsive to low doses of IL-2 (Au-Yeung et al., 2017), they both require a threshold of IL-2 signal that would selectively induce their proliferation after activation (Au-Yeung et al., 2014).

We propose a non-spatial and continuous time ODE to describe T-cell responses with regards to 1) the activation with the APCs, 2) the disappearance, and 3) the expansion of activated T cells during the presence of pathogens. The differentiation of helper T cell is controlled by the APCs Agents, which translates information about the microbial threat to the T cells. However, the activation of these helper T cells is performed by the classical APCs (DC, and B cells) which present MHC II peptide complex (Amsen et al., 2004). The main role of helper T cells is to mediate both the humoral and cellular immunity response by direct receptor binding or the releasing of specific cytokines that boost the immune response (Bianca et al., 2012). Since there is possibility that the helper T cell get activated independently by DC and B cells or both simultaneously, we assigned a variable ( Λ ) to the classic APCs defining the activation protocols controlled by Agents DC and B cells in the MMQE.

The mathematical model describing the dynamics of helper T cells would be a system of nonlinear first order ODEs that capture the additional time dependent change of the different levels of T cells representing linear independent factors depending on: 1) the distribution of the agent of activation, 2) the level of the agent of proliferation (cytokines), and 3) the death of T cells. Thus, it is expressed as: d T t dt = ϕ Th Λ t + σ 2 ( N 2 N 2 + τ 2 ) T t − δ Th T t (6) where the term ϕ Th Λ t performs the estimation of the activated helper T cell with rate parameter ϕ Th and the dependent activation variable Λ t which varies in terms of the activation protocol controlled by DC and B cells. The term σ 2 ( N 2 N 2 + τ 2 ) T t in Eq. (6) describes the stimulation of T helper cells by the cytokine IL-2 with level N 2 , threshold τ 2 and the stimulation rate σ 2 . It is worth noting that the threshold represents the level of immunogenic stimulation required to elicit an immune response, determined by the sum of the negative regulatory mechanisms that work at all levels of the immune system (Guram et al., 2019). Stimulation of the immune system to any level below this threshold may still activate individual cells of the immune system but is insufficient to overcome various regulatory mechanisms and mount an effective systemic immune response. The T cell death factor is modeled by the term − δ Th T t with death rate δ Th .

Germinal centers (GCs) are specialized compartments within the secondary lymphoid organs, where B cells proliferate, differentiate, and mutate their antibody genes (Martínez et al., 2012). Upon exit from the GC, B cells terminally differentiate into plasma cells or memory B cells (Thomas et al., 2019). The GC is divided into two functionally distinct zones, a dark zone, where B cells proliferate, and a light zone, characterized by the presence of follicular dendritic cells and bound antigen and which is associated with B cell selection (Camacho et al., 1998; Victora and Nussenzweig, 2012). Investigating the release of B cells from the GC (Meyer-Hermann et al., 2012) has identified that antigen-retaining B cells differentiate to plasma cells and leave the GC through the dark zone. This is of great importance in Synthetic Biology since it will allow a better quantification of B cells exiting from GC which will provide a better prediction on its differentiation into plasma cells. Additionally, the modeling of B cell maturation enables the characterization of the evolutionary process and competition at the heart of the GC dynamics (Pélissier et al., 2020). It shows that during an immune response one clone (or a small number of clones) of B cells is preferentially selected and proliferated, apparently at random, from a heterogeneous population of cells capable of responding to the given antigen (Williamson et al., 1976; Pélissier et al., 2020). Based on T cell-dependent activation of B cells (Liao et al., 2017; Cox and Brokstad, 2020), we aim to evaluate the dynamic proliferation of B cells in terms of 1) different APCs activating the immune systems and 2) the effect of the interaction between IL-2 and IL-4.

We model the dynamics of B cell population immune response as a function of B-cell activation, disappearance, or differentiation into more antibodies. With regards to their activation, B lymphocytes use B cell receptors (BCRs) to sense the physical features of the antigens so B cell receptor (BCR) signaling is pivotal for optimal B cell activation and differentiation (Guo and Rothstein, 2013; Wan et al., 2015). Hence, BCR interaction with an antigen should overcome the “survival” threshold to achieve full B cell activation (Avalos et al., 2014). In this regards, IL-4 was shown to provide survival or threshold signals to BCR for the activation and proliferation of B cells (Illera et al., 1993; Chung et al., 2002). This makes IL-4 to play an essential role in the activation of mature B cells as a cofactor for CD40L, and antigen stimulation to induce B cell differentiation, proliferation, and antibody secretion (Granato et al., 2014). Furthermore, the B cells begin to proliferate rapidly in the presence of the proper cytokines (IL-2 and IL-4), at their threshold levels (Feldkamp and Carey, 1996).

B cells are well known for their ability to differentiate into antibody secreting plasma cells in response to foreign antigens or pathogens (Lund, 2008). However, recent studies suggested that B cells also appear to regulate both protective and pathologic immune responses by an antibody independent mechanism (Hurdayal et al., 2017). Hence, after B cells get activated by IL-4, in a first stage, they can act as specialized antigen presenting cells, by recognizing pathogens through their specialized receptors and can then present peptide sequences to helper T cells. Consequently, a successful interaction with helper T cells, results in differentiation of B cells into plasma B cells, which release more of the same stimulating antibodies. Therefore, modeling the dynamics of B cell immune response should include frameworks that investigate on 1) the activation of B cells by the cytokine IL-4, 2) the proliferation or stimulation of activated B cells, and 3) the differentiation and death of mature B cells (Bianca et al., 2012; Figueredo et al., 2013). Considering these three points, the model’s equation illustrate the dynamics between B cell ( B t ), T cells ( T t ) and the cytokines IL-2 ( N 2 ), and IL-4 ( N 4 ) described by the following differential equation: d B t dt = ϕ B ( N 4 N 4 + τ 4 ) T t + σ 2 ( N 2 N 2 + τ 2 ) B t − ( ϕ p + δ B ) B t (7) where the term ϕ B ( N 4 N 4 + τ 4 ) performs the activation and proliferation of B cells with respect to IL-4 where N 4 stands for the bounded level of IL-4 to IL-4R. ϕ B and τ 4 represent the activation rate of IL-4 and the binding threshold of IL-4 respectively. The term σ 2 ( N 2 N 2 + τ 2 ) describes the proliferation of B cells by duplication due to IL-2 with N 2 defining the bounded level of IL-2 to IL-2R and, σ 2 and τ 2 represents the duplication rate and the binding threshold of IL-2 respectively. The term − ( ϕ p + δ B ) emphasizes the differentiation and the death of the mature B cells ( B t ) where ϕ p and δ B represent the differentiation and death rate respectively of mature B cells.

Besides being able to recognize antigens such as APC, B cells can be differentiated into a plasma B-cell, whose specialized job is to mass-produce the antibodies that match the activating invader (antigen) and provide immune protection for the human body (Hoffman et al., 2016). Thus, each plasma B-cell makes antibodies unique to only one antigen. Hence, B cells leave the germinal centers at rate ϕ p to become plasma B cells ( P t ) through the ODE (Erwin and Ciupe, 2017): d P t dt = ϕ p B t (8)

B cells maintain immune tolerance and suppress a wrong or excessive immune response to avoid damage to the body by controlling their rate of differentiation, proliferation, and death.

Classical APCs including DC and B cells are a heterogeneous group of immune cells that mediate the cellular immune response by processing and presenting antigens for recognition by certain lymphocytes such as T cells. Therefore, we illustrate the dynamic of the dependent activation variable Λ t in Eq. (6) in terms of the dynamic proliferation of agents DC and B cells as follow.

• In the case where DC Agent makes helper T cells recognize the antigens for its activation (Figure 1A) then, the distribution of the activation variable ( Λ t ) will be defined in respect of the differentiation of agent DC as:

We first obtain estimates for the parameter values and initial conditions for our model from published in-vitro or in-vivo experimental research. All the parameters involved in (Eqs 2–10) are summarized in Table 2. However, by referring to the initial conditions, we provide more details on some parameters such as: the initial number of naïve T cells and B cells, the death rates of T and B cells, the production rate of the plasma B cells, the mean and variance of the lognormal distribution function of the pathogens. The microscopic states of antigens are set to the lognormal distribution with mean μ a g = 0.075 and variance V a g = 1 . This will present a continuous skewed variation of pathogens and will enable the progressive evaluation of the proliferation of lymphocytes as immune responses. The population of T cells in the body is naturally divided into clonotypes which represents the set of cells that have identical T-cell receptors (Lythe and Molina-París, 2018). The total number of cells, such as naïve CD4+ T cells, are large enough to be modeled by ordinary differential equations using appropriate starting points. All naïve CD4+ T cells are about 10 7 in an adult mouse or about 10 11 in an adult human (den Braber et al., 2012). Accordingly, we set the initial number of naive CD4+T cells ( T h 0 ) which is recruited by the antigens ( A a g ) as T h 0 = 805 cells per ml. The death rates of all CD4+T cells are equal to δ T h = 0.01   ( s − 1 ) such that the activation rate of helper T cells is defined as: ϕ T h = T h 0 × δ T h . The adaptive immune response is slow to develop on first exposure to a new pathogen, as specific clones of B cells have to become activated and expand; it can therefore take time before the response is effective (Roberts et al., 2002). Here, we assume that at the initial time when the body is exposed to the pathogens, T cells receive the first signal before communicating with B cells. Therefore, at the initial time, we consider a negligible number of B cells (approximately 100 cells per ml) comparing it to the initial number of T cells. The death rate of mature B cells ( δ B ), the activation rate of B cells ( ϕ 4 ) and the production rate of plasma B cells are defined as: δ B = 0.01   ( s − 1 ) , ϕ B = 1   ( s − 1 ) and ϕ p = 3   ( s − 1 ) respectively. Moreover, due to the fact that the regulation of cytokines is necessary in mediating a pathogen invasion response and only limited information is available on the timing of cytokine production and decay (Thakar et al., 2007), we assume that the regulatory dynamics of all the three lymphocytes (helper T cells, B cells and plasma B cells) are independent to each other.

The model contains 19 parameters which have specific biological meanings. These parameters have been derived from literature, from measurements in different in-vivo experiments which is similar to other immune system studies (Bianca et al., 2012; Gao et al., 2016; Nickaeen et al., 2019). Therefore, to account for uncertainty in computational predictions, we simulated our model by taking account of error margins allowing parameters to vary in the range [ Parameter − σ p ,  Parameter + σ p ] where σ p represents the estimated standard error corresponding to the 95% confidence interval for the model parameters. Varying parameter values in a suitable range while fixing the rest, represents the classical way to do sensitivity analysis (Bianca et al., 2012).

We first evaluate our model immune response resulting from the interactions between IL-2 and IL-4 ligand binding and their receptors. A rapid association rate of IL-2 binding to IL2-R occurs followed by a slower dissociation rate (Figure 3A) for both bounded and unbounded ligands. Specifically, the maximal rate of association for bounded ligands seems to occur after 15,000 s. A rapid maximal rate of dissociation of IL-2 from IL-2R occurs after 6,000 s for unbounded ligands. On the other hand, IL-4R binding shows a much slower association rate of IL-4 ligand as well as a longer delay of dissociation rate of IL-4 from IL-4R (Figure 3B). Quantitatively, the maximum rate of binding of IL-4 to IL-4R or unbinding rate of IL-4 from IL-4R can take a long time [defined in the range (36e+04, 48e+04) seconds] to complete. Note that the number of IL-4 bounded to the IL-4R is double the number of IL-4 unbounded from the IL-4R. Additionally, since the binding or unbinding mechanisms of IL-2 and IL-4 to/from IL-2R and IL-4R occurs simultaneously in response to the presence of the pathogens in the immune system, it is worth noting that the more the bounded number of IL-4 to IL-4R increases, the more the bounded number of IL-2 to IL-2R decreases (Figures 3A,B).

Upon recognition of specific antigens via their T-cell receptor (TCR) in the context of appropriate co-stimulatory signals, T cells clonally expand and traffic to tissues, where they perform effector functions including direct killing of infected cells and secretion of cytokines coordinating local immune responses (Brenchley et al., 2002). Sustained B cell activation also leads to increase in B cell proliferation and differentiation. In response to activation signals, naive mature B cells proliferate and differentiate into effector cells. This is driven by integration of several infection-related signals, including binding of specific antigens to the B cell receptor (BCR) and pattern recognition receptor ligands (Ruprecht and Lanzavecchia, 2006; Pone et al., 2010).

The proliferation of T and B cells increases to a maximal level followed by degradation regardless of the immune system activation agent (Figures 4A1, B1). However, the maximum number of T and B cells produced during their proliferation lifecycle depends on the activation agent. For example, B cell activating agents lead to a higher maximal level of T and B cells proliferation compared to DC agents. In addition, the proliferation of T cells and B cells immune response with respect to the B cell agent activator decreases from its maximum level to zero after about 390000s of proliferation. On the other hand, the proliferation of T and B cells based on DC agent activators reaches steady state during the same period. Conversely, the proliferation levels of plasma B cells (Figure 4C1) reach steady state earlier at 270000s when the activation mechanism is regulated by B cell agent. However, their proliferation rate keeps increasing with time when the immune response is activated by DC agents. In addition to regulating B cell growth and immunoglobulin secretion, IL-4 also affects the proliferation of lymphocytes. The proliferation level of T, B and Plasma B cells tightly depend on the interaction between cytokines IL-2 and IL-4. Hence, T cells, B cells, and Plasma B cells together display a high level of their immune response which is associated with a low concentration of IL-4 versus high concentration of IL-2 (Figures 4A2, B2, C2). In other words, the concentration of IL-4 decreases proportionally with increase in proliferation levels of T and B cells which automatically influences the growth of the Plasma B cells. However, the level of proliferation of T and B cells is highly influenced by the activation protocol defined by the joint activation with B cell and DC. Thus, T and B cells present a higher level of proliferation comparing it to the case where the immune system was activated independently by B cell or DC (Figures 5A,B). Similarly, the proliferation of plasma B cells reached a high level of steady state after 270000s of proliferation time (Figure 5C). Supplementary Figure S1–S3 provides additional plots of similar dynamic profiles as in Figure 3 and 4 of lymphocyte proliferation from different proportions of IL2-IL4 concentration levels.

We next calibrate the model by investigating the effects of simultaneously varying IL-2 and IL-4 concentration levels that are likely to influence the productivity growth of the lymphocytes during an immune response. Due to the permanent exposure of lymphocytes to pathogens in the immune system, our aim is to determine the rate of proliferation of lymphocytes during the immune response in relation to the abundance of cytokines IL-4 and IL-2 acting jointly on individual cells. We investigate the marginal effects of eight different intervals of concentrations of IL-2 and IL-4 in pM as follows: C 0 = [ 0 − 0.75 ] , C 1 = [ 0.76 − 1.5 ] , C 2 = [ 1.51 −   2.25 ] , C 3 = [ 2.26 − 3 ] , C 4 = [ 3.01 − 3.75 ] , C 5 = [ 3.76 − 4.5 ] , C 6 = [ 4.51 − 5.25 ] , C 7 = [ 5.26 − 6 ] . From each interval, we 1) generate uniformly 50 samples of concentration of IL-2 and IL-4 to compute the proliferation of lymphocytes during immune response based on (Eqs 6–10, Eq. 2) and identify the maximal level of each lymphocyte produced during their proliferation. Accordingly, we collect in total 50 data points of maximal proliferation levels of each of the lymphocytes associated with each interval.

Denoting y 1 , ... , y k for k = 50 as the maximal proliferation levels of lymphocyte collected in each range of cytokine concentration range R   = [ [ c y t o k i n e ] m i n , [ c y t o k i n e ] m a x   ] , we define the rate of proliferation growth ( γ p ) as: γ p = mean ( y 1 , ⋯ , y k ) sum ( y 1 , ⋯ , y k ) × 100

Figure 6 displays the bar plots representing both lymphocyte proliferation growth in terms of number of cells as well as the associated rates of proliferation growth across various cytokine concentration intervals. The proliferation of T cells, B cells, and plasma B cells is left skewed with decrease in the concentration levels of IL-2 and right skewed with increase in IL-4 concentration indicating opposing trends with plasma B cells showing the largest difference.

Our computational model aims to explain biological phenomenon of T cells activation followed by the proliferation of T, B and plasma B cells, as well as the antagonistic effects of IL-2 and IL-4 on lymphocyte proliferation in the presence of antigens. CD4+ T cell functions include activating other immune cells, releasing cytokines, and helping B cells to produce antibodies as well as differentiate into plasma B cells. Due to CD4+ T cell’s ability to shape, activate and regulate the adaptive immune response (Henninger et al., 2010; Swain et al., 2012), we validate our model by targeting T cells as a variable of interest based on their relative importance in immune response (Henninger et al., 2010). We analyze data from two different experimental studies (Morris et al., 2000; Ahmadi et al., 2008) to help validate 1) the comparison of CD4+ T cell proliferation with the activation of the naïve T cells initiated by different APCs: activated B cells, and mature DC (Supplementary Figure S4); 2) the influence of the antagonistic effect of cytokines IL-2 and IL-4 on the proliferation of CD4+ T cells (Supplementary Figure S5).

Figure 7A summarizes a direct head-to-head comparison of proliferation of CD4+ T cells between CD40L-activated B cells and mature DCs with regards to presenting antigens to the naïve T cells for activation (Supplementary Figure S2 for more details on the first study). It is worth noting that CD40 ligand (CD40L) is the critical membrane-expressed molecule responsible for T cell dependent B-cell activation denoted as CD40L-activated B cells serving as effective antigen APC. These data indicate that primary CD40L expanded B cells are efficient presenters of antigens to T cells and characterized CD40L-activated B cells as more influential APC compared to mature DCs. Figures 7B,C (Supplementary Figure S3), show differences in percentage change of CD4+ T cells and IL-2 signaling expression respectively between joint IL-4C signaling and antigen presence, and only antigen production. We observe a higher CD4+ T cell production in the absence of IL-4C stimulation (Figure 7B). The higher production of antibody by B-cells recognized by the increased presence of mouse antigen (GαMδ) resulted in the largest increase in the level of CD4+ T cells compared to the untreated case (case with no GαMδ). However, the treatment with IL-4C inhibited GαMδ-induced IL-2 response (Figure 7C) systematically shows a reduction of CD4+ T cells and proliferation levels since IL-2 is the cytokine growth factor of CD4+ T cells (Figure 7B). The above results provide evidence of the antagonistic effect of IL-4 on IL-2 by decreasing not only CD4+ T cell activation but also the T cell production of IL-2.

We use a hybrid ODE models with a stochastic component called MMQE to capture the variability in dynamic interaction effects resulting from communication between extracellular and intracellular dependent agents (Figure 2) within an immune response environment. Our model combines agents using network-based rules within an ODE system with specific initial conditions presented in the subsection of Numerical Evaluation. We implemented the given ODEs and estimate parameters analytically using MATLAB (Macal and North, 2005; Cilfone et al., 2015; Martin and Schlüter, 2015; Hunter et al., 2018). The codes used in this study were written using MATLAB version MATLAB_R2021a and is available at https://github.com/atiteyk2/Downregulation_Immunology.git.