# A Fast Second-Order Explicit Predictor-Corrector Numerical Technique To
Investigating And Predicting The Dynamic Of Cytokine Levels And Human Immune
Cell Activation In
Response To Gram-Positive Bacteria: Staphylococcus Aureus

Eric Ngondiep${}^{\text{\,1}}$, Ariane Njomou Ndantouo${}^{\text{\,2,3}}$,
George Mondinde Ikomey${}^{\text{\,2,3}}$

(${}^{\text{\,1\,}}$Department of Mathematics and Statistics, College of
Science, Imam Mohammad Ibn Saud
Islamic University (IMSIU), $90950$ Riyadh $11632,$ Saudi Arabia.

${}^{\text{\,2\,}}$School of Health Sciences, Catholic University of Central
Africa, $1110$ Messa-Yaounde, Cameroon.

${}^{\text{\,3\,}}$Center for the Study and Control of Communicable Diseases
(CSCCD), Faculty of Medicine
and Biomedical Sciences, University of Yaounde 1, $8445$ Yaounde, Cameroon.

Correspondence email: ericngondiep<EMAIL_ADDRESS>
Abstract. This paper develops a second-order explicit predictor-corrector
numerical approach for solving a mathematical model on the dynamic of cytokine
expressions and human immune cell activation in response to the bacterium
staphylococcus aureus (S. aureus). The proposed algorithm is at least zero-
stable and second-order accurate. Mathematical modeling works that analyze the
human body in response to some antigens have predicted concentrations of a
broad range of cells and cytokines. This study deals with a coupled cellular-
cytokine model which predicts cytokine expressions in response to gram-
positive bacteria S. aureus. Tumor necrosis factor alpha, interleukin 6,
interleukin 8 and interleukin 10 are included to assess the relationship
between cytokine release from macrophages and the concentration of the S.
aureus antigen. Ordinary differential equations are used to model cytokine
levels while the cellular responses are modeled by partial differential
equations. Interactions between both components provide a more robust and
complete systems of immune activation. In the numerical simulations, a low
concentration of S. aureus is used to measure cellular activation and cytokine
expressions. Numerical experiments indicate how the human immune system
responds to infections from different pathogens. Furthermore, numerical
examples suggest that the new technique is faster and more efficient than a
large class of statistical and numerical schemes discussed in the literature
for systems of nonlinear equations and can serve as a robust tool for the
integration of general systems of initial-boundary value problems.

Keywords: mathematical model on dynamic of mixed cellular-cytokine, cytokine
expressions, human immune cell activation, S. aureus, a second-order explicit
predictor-corrector numerical technique, numerical simulations.

AMS Subject Classification (MSC). 65Lxx.

##  1 Introduction and motivation

Mathematical models arising from the complex biological systems are important
decision tools that allow to elucidate emergent properties of intricate
biological pathways with the human body [16, 44, 6, 48, 42]. The mathematical
modeling works based on the human immune cells in response to pathogens have
predicted the concentrations of a broad range of cells including: B-Cell,
dentritic, macrophages and plasma cells. The leukocytes of the immune system
are usually the neutrophils, monocytes, basophiles eosinophiles and
lymphocytes. They play a very important role in recognizing and fighting
infections. Their main function is to eliminate pathogens by phagocytosis and
cell debris through engulfment and chemical degradation. The interaction
between the bacterium and host during infection leads to immune cell
activation and often to a change in the relative ratio of leukocyte sub-
populations of the innate immune and adaptive immune responses [4, 7]. Raman
spectroscopy based on spectroscopic fingerprint in clinical samples are
efficient approaches to identify and differentiate these subsets of leukocytes
[69, 40, 9]. Furthermore, this technique also detects immune cells activation
and apoptosis [36, 8, 68, 3].

Adaptive immunity deals with a long-term specific response initiated to
eliminate a well-known pathogen whereas innate immunity considers many lines
of defense starting with: saliva, skin, various secretions and ending with
non-specific leukocytes [43, 33]. Additionally, the innate immune system is
capable to recognize and tag a wide set of antigens. The inability of innate
immune response to kill the pathogen leads to the activation of the adaptive
immune system, which is firstly composed with T cells and B cells. The
functions of B cells is to produce antibodies to neutralize the pathogen and
subsequently destroy them. Furthermore, T and B cells contribute in the
production of small signaling proteins released by leukocytes, so called
cytokines, which facilitate the communication between immune cells. These
cytokines include: tumor necrosis factor alpha (TNF$\alpha$), interleukin 6
(IL6), interleukin 8 (IL8) and interleukin 10 (IL10) [80, 76]. TNF$\alpha$ is
fundamental to the acute phase reaction during inflammatory response. IL6 is a
pro-inflammatory cytokine that is involved in inflammation and homeostasis
processes and can also act as an anti-inflammatory cytokine through its
inhibitory effects on TNF$\alpha$. In addition, IL6 plays a crucial role in
the recruitment of T cells and in the production of T and B cells during
inflammation and the delay in apoptosis of T cells [24]. In the literature
[78, 64, 30, 14], the authors have shown that IL8 is induced by TNF$\alpha$,
inhibited by IL10, that is an anti-inflammatory cytokine critical in the
regulation of immune response and IL8 is also involved in the recruitment of
basophils, neutrophils and T cells. Moreover, IL10 maintains homeostasis and
prevents host damage during infection while acting as an immuno-regulator.
Furthermore, this anti-inflammatory cytokine limits the production of pro-
inflammatory cytokine including IL6, while down-regulating the expression of
TNF$\alpha$, T helper type 1 cytokines, and major histocompatibility complex
class 2 molecules [24, 31]. IL6, IL8 and IL10 exert many effects on the immune
system, hematopoiesis and acute-phase response. In [46, 79], the authors
showed that monocytes, lymphocytes, cancer cells as well as macrophages have
been documented for the production and secretion of IL6, IL8 and IL10. These
cytokines in a para-crime manner induce in vitro growth of melanoma cells,
prostate, lung, ovarian, kidney and cervical cancers. Specifically, the
disease level and the activity of immune system mechanisms modulated by some
interleukins such as: IL6, IL8 and IL10, allow to evaluate the efficiency of
the treatment and prognosis in course of malignancy. For more details, we
refer the readers to [34, 66, 27, 25, 35].

Although present in all bacterial cell walls, peptidoglycan (PePG) helps to
stabilize cell structure and shape. PePG protects cells from bursting in
response to environmental stressors. Similarly, lipoteichoic acid (LTA), known
as a gram-positive bacterium, is composed with a glycolipid covalently bound
and a hydrophilic glycerophosphate polymer. Innate immune system deals with
peptidoglycan and lipoteichoic acid which are able to trigger the systemic
release of cytokines [29, 13]. Lipopolysaccharide (LPS), a hallmark of gram-
negative bacteria, commonly called endotoxin, consists of a membrane-anchor
lipid, glycan polymer and oligosaccharide core. It lies in the class of the
most potent immuno-stimulants. It has been proven [22, 72] that pro-
inflammatory and anti-inflammatory cytokine activation in gram-negative
bacterial infections are mainly driven by endotoxin. Analogously, one cause of
the pro-inflammatory cytokine activity follows from the physiological
recognition of the lipid component by the immune system [37, 1, 45, 72]. PePG
and LTA are the main sources of activation of cytokines in response to gram-
positive bacterial infections whose levels can be predicted by the use of
mathematical models [13, 11, 6, 23, 38].

Complex systems of differential equations arising in mathematical biology and
physics play a major role in various branches of science and engineering by
describing their interaction among the diffusion transport and reaction
mechanism. A large class of unsteady nonlinear differential equations such as:
the mathematical models of covid-19, Navier-Stokes problems, Stokes-Darcy
models, convection-diffusion-reaction equations, time-fractional equations,
advection-diffusion problem, mathematical models of dynamic of poverty and
corruption, etc.., are introduced to capture some features of several real-
world biological and physical phenomena [62, 17, 52, 41, 57, 58, 63, 59]. For
example: traffic flow, plasma physics, fluid mechanics, optical fibers,
chemical kinetics, mathematical immunology, population dynamics, neutron
nuclear reaction, financial derivatives pricing [15, 74]. A wide range of
applications requires the analytical solutions of such equations, which allow
to well understand the interesting features, novel phenomena and intrinsic
mechanism hidden in the dynamic systems [2]. Since the exact solutions of
these equations only exist under some restrict conditions, develop efficient
and reliable numerical techniques for such systems of unsteady nonlinear
equations is of great interest. In the literature, researchers have described
and analyzed abundant statistical and numerical approaches in an approximate
solutions [21, 70, 50, 53, 26, 28, 56, 60, 65, 55, 32, 5, 73, 77]. Most of the
existing mathematical models of complex biological systems deal with the
cytokine response to lipopolysaccharide and the basic pathway of the standard
immune response to S. aureus. Although these models represent responses to
several pathogens such as: gram-positive and gram-negative bacteria, the
cytokine expressions are computed as functions of activated macrophages. This
suggests similar interactions between cytokines and macrophages [75, 67, 5].

In this paper, the mathematical model considers a system of ordinary
differential equations (ODEs) coupled with partial differential ones (PDEs) to
emulating the relationship between the pathogen S. aureus, cytokines and cells
with the human immune system. Specifically, a second-order explicit predictor-
corrector scheme is developed to investigating and predicting the cells-based
activation and cytokine expressions induced by the bacterium S. aureus. Both
theoretical and numerical results indicate that the proposed algorithm is less
time consuming, faster and more efficient than a large class of statistical
and numerical approaches applied to such complex systems of differential
equations [28, 10, 54, 26, 51, 65]. In addition, our results show that the
developed numerical method can be observed as a robust tool to predict ex vivo
and in vivo experimental data induced by a given antigen. The highlights of
this work are the following:

i.

mathematical formulation of the dynamic of cytokine concentrations and
activation of human immune cells in response to the bacterium S. aureus,

ii.

construction of the explicit predictor-corrector scheme and analysis of
stability and convergence rate of the proposed algorithm,

iii.

numerical examples to confirm the theory.

The paper is organized as follows: we present the mathematical model on the
dynamic of cytokine levels and human immune cell activation in response to the
propagation of S. aureus in Section 2. In Section 3, we develop the two-step
explicit numerical technique for solving the considered problem and we provide
both stability analysis and convergence order of the new algorithm. Section 4
considers some numerical simulations whereas the general conclusions and
future works are described in Section 5.

##  2 Mathematical formulation of the model

This section deals with a coupled inflammatory cytokine model and cellular one
as proposed in [5, 23] (see also Figure 1). This combination is so called
mixed cellular-cytokine model and it represents an extension of some previous
works discussed in the literature [6, 48, 16, 17, 41] by considering the
inter-connectivity between the cytokine and cellular responses and captures ex
vivo and potential in vivo. The cellular model given in [67] predicts the
relationships between activated and resting macrophages, antibodies and S.
aureus whereas the cytokine model analyzes how changes in activated
macrophages and resting ones impact cytokine evolution for TNF$\alpha$, IL6,
IL8 and IL10 [5, 67]. Both models operate in the same time frame whenever the
concentrations of the initial bacteria and the activated macrophages are
scaled to match those reported in the ex vivo studies [5, 78, 75]. The coupled
cellular-cytokine model is described by a system of PDEs and ODEs and it will
be efficiently solved using a new second-order explicit predictor-corrector
method. The cellular model is represented by the PDEs while the cytokine one
considers the ODEs. It is worth noticing to mention that the gram-positive
bacteria (for instance: S. aureus) causes the wide-spread inflammation and
septic shock that is primarily due to the function LPA and PePG during an
inflammatory response [47, 22, 24]. The work carried out by LTA and PePG
induces the cytokine concentration in the host’s innate and adaptive immune
response to S. aureus pathogen. In [22, 24, 47], researchers have established
that the action mechanisms in the host which include: neutrophil flux,
phagocytosis and Sbi protein activation differ from the wall components and
the given pathogen, but indicate the same inflammatory cytokine responses.

|
---|---
|
Figure 1: Representation of: cytokine model, cellular one, mixed cellular-
cytokine model and Raman spectroscopic analysis

###  2.1 Cellular model

The cellular model discussed in [76, 73] predicts the activation of the
acquired immune response to S. aureus as a function of concentrations of
lymphocytes, antibodies and plasma cells, denoted by: (T,B), F and A,
respectively, together with activated and resting macrophages expressions,
bacteria and antibodies spatially distributed in $1$ $cm^{3}$ of lump tissue
$(x=(x_{1},x_{2},x_{3}))$ and represented by: M, A and F, respectively. The
variation of concentrations of lymphocytes, antibodies and plasma cells is
observed at the nearest of the lymph node. In this paper, we assume that the
activated macrophages behave as antigens presenting cells and they move to the
nearest lymph node where the specific response is triggered. So, these
macrophages are obtained by interacting in space with the antigen in the
tissue $(M_{A}(x,t))$ and at the lymph node $(M_{A}^{L}(t))$ by interacting
with the lymphocyte expressions varying only in time. Thus, the dependent
variables dealing with the cellular model are described as follows.

  * •

Spatial variables $(pg/mm^{3})$: S. aureus pathogen ($y_{1}(x,t)$), resting
macrophages ($y_{2}(x,t)$), activated macrophages ($y_{3}(x,t)$) and specific
antibodies $(F(x,t))$.

  * •

Temporal variables $(pg/mm^{3})$: plasma cells $(P(t))$, T-lymphocytes
$(T(t))$, antibodies $(F^{l}(t))$, B-lymphocytes $(B(t))$ and average
activated macrophages $(M_{l}^{a}(t))$.

###  2.2 Coupled cellular-cytokine model

This work considers the dynamic of mixed cellular-cytokine model given in [71]
and described by a system of ordinary and partial differential equations.
Assuming that the average tissue concentration obtained from the cellular
model is the initial condition, the expected concentrations of TNF$\alpha$,
IL6, IL8 and IL10 are functions of resting macrophages $y_{2}(x,t)$ and the
activated ones $y_{3}(x,t)$. In some previous works, the authors developed
this model to study the response of cytokines to LPS. Specifically, their
findings have shown that a combination of cytokine model and the cell-based
one is ”well-defined” because LPS and pathogen-associated molecular patterns
(PAMP) induce similar pro-inflammatory and anti-inflammatory responses [20].
The considered coupled model deals with the following unsteady variables:
$y_{1}(x,t)$ indicates the bacteria, $y_{2}(x,t)$ and $y_{3}(x,t)$ denote
resting and activated macrophages, respectively, while $y_{4}(t)$, $y_{5}(t)$,
$y_{6}(t)$ and $y_{7}(t)$ represent: tumor necrosis factor alpha, interleukin
6, interleukin 8 and interleukin 10, respectively. We remind that tumor
necrosis factor alpha, interleukin 6 and interleukin 8 are pro-inflammatory
whereas interleukin 10 is an anti-inflammatory cytokine that plays an
important role in the regulation of immune responses. Furthermore, interleukin
10 reduces the production of some pro-inflammatory cytokines (for example:
interleukin 6 and tumor necrosis factor alpha).

Let $\Omega$ be a bounded domain of $\mathbb{R}^{d}$, where $d=2$ or $3$, and
$T$ be a positive real number. The initial-boundary value problem modeled by
the S. aureus bacterium ($y_{1}(\cdot)$) is defined as

$\left\\{\begin{array}[]{ll}\frac{\partial y_{1}}{\partial
t}=\beta_{1}y_{1}(1-k_{1}^{-1}y_{1})-\mu_{1}y_{1}-\lambda_{2}y_{1}y_{2}-\lambda_{3}y_{1}y_{3},&\hbox{on
$\Omega\times[0,\text{\,}T]$}\\\ \text{\,}\\\
y_{1}(x,0)=\overline{y}_{1},\text{\,\,\,on\,\,\,}\overline{\Omega},\text{\,\,\,\,\,\,\,\,\,}\frac{\partial
y_{1}}{\partial t}(x,t)=0,&\hbox{on
$\partial\Omega\times[0,\text{\,}T]$.}\end{array}\right.$ (1)

Here: $\beta_{1}y_{1}(1-k_{1}^{-1}y_{1})$ denotes the logistic growth of the
bacteria, -$\mu_{1}y_{1}$ indicates the natural decay rate of S. aureus
without immune system processes through the natural decay coefficient
$\mu_{1}$, the terms: -$\lambda_{2}y_{1}y_{2}$ and -$\lambda_{3}y_{1}y_{3}$
describe the phagocytosis of the pathogen S. aureus through resting
macrophages and activated ones. the coefficients: $\beta_{1}$, $k_{1}$,
$\mu_{1}$, $\lambda_{2}$ and $\lambda_{3}$, are positive constants which
represent the carrying capacity, replication rate, natural decay and declined
caused by resting and activated macrophages, respectively.

Up-regulation and down-regulation for each cytokine are modeled by the use of
sigmoidal functions defined as

$H_{c}^{u}(y)=\frac{y}{\eta_{xy}+y}\text{\,\,\,\,\,or\,\,\,\,\,}H_{c}^{d}(y)=\frac{\eta_{xy}}{\eta_{xy}+y},$
(2)

where $y$ is the cytokine inducing down-regulation (superscript: d) or up-
regulation (superscript: u) of cytokine ”c”, $\eta$ denotes the half-maximum
value. We add sigmoidal functions given by $(\ref{2})$ in all the equations
dealing with cytokines to describe the relationship between cytokines (see
Table 1 below).

The mathematical model of resting macrophages ($y_{2}(\cdot)$) response to the
bacterium S. aureus is given by

$\left\\{\begin{array}[]{ll}\frac{\partial y_{2}}{\partial
t}=\mu_{2}y_{2}(1-y_{2m}^{-1}y_{2})-\left[\gamma_{3}+k_{6}H_{c}^{u}(y_{4})H_{c}^{d}(y_{7})\right]y_{1}y_{2},&\hbox{on
$\Omega\times[0,\text{\,}T]$}\\\ \text{\,}\\\
y_{2}(x,0)=\overline{y}_{2},&\hbox{on $\overline{\Omega}$}\end{array}\right.$
(3)

where $\mu_{2}y_{2}(1-y_{2m}^{-1}y_{2})$ represents the growth of the resting
macrophages and
-$\left[\gamma_{3}+k_{6}H_{c}^{u}(y_{4})H_{c}^{d}(y_{7})\right]y_{1}y_{2}$
designates macrophages activation at the rate $\gamma_{3}$ and $k_{6}$ in
response to S. aureus and taking into account the cytokines $y_{4}(t)$ and
$y_{7}(t)$ influences, respectively, $\mu_{2}$ is the influx rate associated
with $y_{2}(x,t)$, $y_{2m}$ is the maximum of $y_{2}(x,t)$ on the domain
$\Omega\times[0,\text{\,}T]$, $\overline{y}_{2}$ is the initial condition
defined as the average of resting macrophages in the tissue. This average is
an outcome of the cellular model simulation over $24h$.

Activated macrophages ($y_{3}(\cdot)$) are modeled as follows

$\left\\{\begin{array}[]{ll}\frac{\partial y_{3}}{\partial
t}=-\mu_{3}y_{3}+\left[\gamma_{3}+k_{6}H_{c}^{u}(y_{4})H_{c}^{d}(y_{7})\right]y_{1}y_{2},&\hbox{on
$\Omega\times[0,\text{\,}T]$}\\\ \text{\,}\\\
y_{3}(x,0)=\overline{y}_{3},&\hbox{on $\overline{\Omega}$.}\end{array}\right.$
(4)

In equation $(\ref{4})$, the term -$\mu_{3}y_{3}$ denotes the decay of
$y_{3}(x,t)$ at rate $\mu_{3}$,
$\left[\gamma_{3}+k_{6}H_{c}^{u}(y_{4})H_{c}^{d}(y_{7})\right]y_{1}y_{2}$ is
the macrophages activation at the rates $\gamma_{3}$ and $k_{6}$ considering
the influence of the cytokines $y_{4}(t)$ and $y_{7}(t)$. $\overline{y}_{3}$
is the initial condition, that is defined as the average of activated
macrophages in the tissue. That is, an outcome of the cell-based model
simulations during $24h$. Furthermore, the initial and boundary (of the
tissue) concentrations of both $y_{2}(\cdot)$ and $y_{3}(\cdot)$ are constant.

The dynamic of tumor necrosis factor alpha ($y_{4}(\cdot)$) is described as

$\left\\{\begin{array}[]{ll}\frac{dy_{4}}{dt}=k_{7}H_{c}^{d}(y_{5})H_{c}^{d}(y_{7})y_{3}-k_{2}(y_{4}-q_{4}),&\hbox{on
$\Omega\times[0,\text{\,}T]$}\\\ \text{\,}\\\
y_{4}(0)=\overline{y}_{4}:=0,&\hbox{}\end{array}\right.$ (5)

where, the first term: $k_{7}H_{c}^{d}(y_{5})H_{c}^{d}(y_{7})y_{3}$ denotes
the down-regulation interactions that cytokines $y_{5}(\cdot)$ and
$y_{7}(\cdot)$ have with the growth of $y_{4}(\cdot)$ (at rate $k_{7}$)
mediated by the average concentrations of $y_{3}(\cdot)$, the second term:
$-k_{2}(y_{4}-q_{4})$ represents the natural decay of $y_{4}(\cdot)$ at rate
$k_{2}$. $q_{4}>0$, is a constant so called a resting level while
$\overline{y}_{4}$ is the initial condition. Equation $(\ref{5})$ suggests
that the rate of change of $y_{4}(\cdot)$ depends on $y_{3}(\cdot)$.

The equation that models the evolution of interleukin 6 ($y_{5}(\cdot)$) is
given by

$\left\\{\begin{array}[]{ll}\frac{dy_{5}}{dt}=(k_{8}+k_{9}H_{c}^{u}(y_{4}))H_{c}^{d}(y_{5})H_{c}^{d}(y_{7})y_{3}-k_{3}(y_{5}-q_{5}),&\hbox{on
$\Omega\times[0,\text{\,}T]$}\\\ \text{\,}\\\
y_{5}(0)=\overline{y}_{5}:=0,&\hbox{}\end{array}\right.$ (6)

where, the first term:
$(k_{8}+k_{9}H_{c}^{u}(y_{4}))H_{c}^{d}(y_{5})H_{c}^{d}(y_{7})y_{3}$
designates the interaction between $y_{4}(\cdot)$ (up-regulation) and
$y_{7}(\cdot)$ (down-regulation) affecting the $y_{5}(\cdot)$ production (at
rate $k_{9}$) which also induces auto-negative feedback. The term:
$-k_{3}(y_{5}-q_{5})$ indicates the natural decay (at rate $k_{3}$) of
$y_{5}(\cdot)$ toward resting level: $q_{5}>0$.

The time-dependent equation modeled by interleukin 8 ($y_{6}(\cdot)$)
expressions is defined as

$\left\\{\begin{array}[]{ll}\frac{dy_{6}}{dt}=(k_{10}+k_{11}H_{c}^{u}(y_{4}))H_{c}^{d}(y_{7})y_{3}-k_{4}(y_{6}-q_{6}),&\hbox{on
$\Omega\times[0,\text{\,}T]$}\\\ \text{\,}\\\
y_{6}(0)=\overline{y}_{6}:=0.&\hbox{}\end{array}\right.$ (7)

Here, the term: $(k_{10}+k_{11}H_{c}^{u}(y_{4}))H_{c}^{d}(y_{7})y_{3}$ denotes
the interaction between the opposing effects of $y_{4}(\cdot)$ (up-regulation)
at rate $k_{11}$ and $y_{7}(\cdot)$ (down-regulation) at rate $k_{10}$,
simulating growth of $y_{6}(\cdot)$ at a rate proportional to the average
concentration of $y_{3}(\cdot)$ production whereas the second term:
$-k_{4}(y_{6}-q_{6})$ represents the natural decay of $y_{6}(\cdot)$ at rate
$k_{4}$, toward resting level: $q_{6}>0$.

Finally, the ODE describing the dynamic of interleukin 10 ($y_{7}(\cdot)$)
concentrations is given by

$\left\\{\begin{array}[]{ll}\frac{dy_{7}}{dt}=(k_{12}+k_{13}H_{c}^{u}(y_{5}))y_{3}-k_{5}(y_{7}-q_{7}),&\hbox{on
$\Omega\times[0,\text{\,}T]$}\\\ \text{\,}\\\
y_{7}(0)=\overline{y}_{7}:=0.&\hbox{}\end{array}\right.$ (8)

In this equation, the first term: $(k_{12}+k_{13}H_{c}^{u}(y_{5}))y_{3}$ is
the up-regulation of $y_{7}(\cdot)$ due to $y_{5}(\cdot)$ (at rate $k_{13}$)
and average concentration of $y_{3}(\cdot)$ (at rate $k_{12}$). The second
term: $-k_{5}(y_{7}-q_{7})$ denotes the natural decay of $y_{7}(\cdot)$ at
rate $k_{5}$. $q_{7}>0$ is the resting level.

Moreover, the authors [12, 67] established that the migration of both S.
aureus and macrophages should be observed as diffusion. The transfer rate of
the cells from one side to another one is represented by the diffusion
coefficients and it is proportional to the concentration gradient of the
particles (cells). For the sake of convenience, the medium is assumed to be
isotropic and it has the same diffusion coefficient in each direction [67].
That is, the diffusion term is modeled via quantity $D_{r}\Delta y_{s}$, where
$r\in\\{MA,MR\\}$ and $s\in\\{A,R\\}$. $D_{r}$ denotes the diffusion
coefficient of the macrophages in the tissue, while $\Delta$ refers to the
Laplacian operator. Following the works discussed in [76, 11, 73] on the cell-
based model, the average numbers of $y_{2}(\cdot)(t)$ and $y_{3}(\cdot)(t)$
are computed by integrating the resulting concentrations of each one over the
domain $\Omega$. That is,

$\overline{y}_{l}(t)=\frac{1}{|\Omega|}\int_{\Omega}y_{l}(x,t)dx,$ (9)

where $l=2,3$, and $|\Omega|$ denotes the volume of $\Omega$.

###  2.3 Description of parameters

This subsection deals with the parameters included in the mixed cellular-
cytokine model and described by the system of nonlinear ODEs-PDEs
$(\ref{1})$-$(\ref{8})$. Since the human immune response modeling by these
coupled equations seems to be too complex, for the sake simplicity, we assume
in this paper that the cellular parameters, cytokine half-maximum value and
hill function exponent parameters are constants.

$D_{1}:$ bacteria diffusion coefficient |  $D_{2}:$ RM diffusion coefficient
---|---
$D_{3}:$ AM diffusion coefficient |  $\beta_{1}:$ replication rate of the bacteria
$\lambda_{12}:$ destruction rate of opsonized bacteria by RM |  $\mu_{1}:$ natural decay rate of the bacteria
$\mu_{2}:$ natural decay rate of RM |  $\mu_{3}:$ natural decay rate of AM
$\gamma_{3}:$ rate of active RM |  $\lambda_{2}:$ activation rate of macrophages
$\lambda_{3}:$ destruction rate of bacteria by AM |  $k_{1}:$ carrying capacity of the bacteria
$\lambda_{13}:$ destruction rate of opsonized bacteria by AM |

where AM=$y_{3}(\cdot)$:= activated macrophages, RM=$y_{2}(\cdot)$:= resting
macrophages and TNF$\alpha$=$y_{4}(\cdot)$:= tumor necrosis factor alpha. The
above parameters determine the different rates of decay or growth for the
bacteria and cells.

$\alpha_{1}:$ migration rate of RM to site of infection |  $k_{3}:$ activation rate of interleukin 6 (per hrs)
---|---
$\alpha_{2}:$ migration rate of AM to site of infection |  $k_{4}:$ activation rate of interleukin 8 (per hrs)
$k_{2}:$ activation rate of TNF$\alpha$ (per hrs) |  $k_{5}:$ activation rate of interleukin 10 (per hrs)
$k_{6}:$ activation rate of RM influenced by $y_{4}(\cdot)$ |  $k_{7}:$ upregulation of $y_{4}(\cdot)$ by AM
$k_{8}:$ upregulation of $y_{5}(\cdot)$ by AM |  $k_{9}:$ upregulation of $y_{5}(\cdot)$ by TNF$\alpha$
$k_{10}:$ upregulation of $y_{6}(\cdot)$ by AM |  $k_{11}:$ upregulation of $y_{6}(\cdot)$ by TNF$\alpha$
$k_{12}:$ upregulation of $y_{7}(\cdot)$ by AM |  $k_{13}:$ upregulation of $y_{7}(\cdot)$ by interleukin 6

We recall that the constants $k_{sl}$ denote the rate of change in the up-
regulation rate of the cytokine secreted from activated macrophages while
$k_{s}$ represent the activation or elimination of cytokine. The initial
values of these parameters are based on the predicted conditions of the model
activated with low dose of the bacterium S. aureus.

$q_{4}:$ concentration of TNF$\alpha$ in absence of pathogem |  $q_{5}:$ concentration of IL6 in absence of antigen
---|---
$q_{6}:$ concentration of IL8 in absence of bacteria |  $q_{7}:$ concentration of IL10 in absence of microbes

where IL6=$y_{5}(\cdot)$:= interleukin 6, IL8=$y_{6}(\cdot)$:= interleukin 8
and IL10=$y_{7}(\cdot)$:= interleukin 10. it worth mentioning that the source
terms $q_{l}$ are also used to set the initial conditions for each cytokine.
The values of these parameters are obtained from the initial predicted
conditions of the model in the absence of the bacterium S. aureus. The half-
maximum value parameter $\eta_{sl}$ describes the effective cytokine
concentration at which targeted cytokine activity should reach half-maximum
with units of $pgmL^{-1}$.

Table 1

$\eta_{45}:$ HMV assoc. down-reg of TNF$\alpha$ by IL6 |  $H_{y_{7}}^{u}(y_{5}):$ HFE assoc. up-reg of IL10 by IL6
---|---
$\eta_{57}:$ HMV assoc. down-reg of IL6 by IL10 |  $H_{y_{5}}(y_{5}):$ HFE assoc. auto-neg feedback of IL6
$\eta_{54}:$ HMV assoc. up-reg of IL6 by TNF$\alpha$ |  $H_{y_{6}}^{u}(y_{4}):$ HFE assoc. up-reg of IL8 by TNF$\alpha$
$\eta_{67}:$ HMV assoc. down-reg of IL8 by IL10 |  $H_{y_{4}}^{d}(y_{7}):$ HFE assoc. down-reg of TNF$\alpha$ by IL10
$\eta_{47}:$ HMV assoc. down-reg of TNF$\alpha$ by IL10 |  $H_{y_{5}}^{u}(y_{4}):$ HFE assoc. up-reg of IL6 by TNF$\alpha$
$\eta_{55}:$ HMV assoc. auto-neg feedback of IL6 |  $H_{y_{5}}^{d}(y_{7}):$ HFE assoc. down-reg of IL6 by IL10
$\eta_{64}:$ HMV assoc. up-reg of IL8 by TNF$\alpha$ |  $H_{y_{6}}^{d}(y_{7}):$ HFE assoc. down-reg of IL8 by IL10
$\eta_{75}:$ HMV assoc. up-reg of IL10 by IL6 |  $H_{y_{4}}^{d}(y_{5}):$ HFE assoc. down-reg of TNF$\alpha$ by IL6

where HMV:=half-max value, assoc. down-reg:= associated down-regulation,
HFE:=hill function exponent, assoc. up-reg:= associated with up-regulation,
assoc. auto-neg:=associated with auto-negative.

In the literature [5, 67], the authors discussed the values of the above
parameters and the obtained results are provided in the following Table 2.

Table 2

parameter | value | unit | parameter | value | unit
---|---|---|---|---|---
$D_{1}$ | $3.7\times 10^{-15}$ | $mm^{3}/day$ | $\beta_{1}$ | $2.0$ | $day^{-1}$
$D_{2}$ | $4.32\times 10^{-2}$ | $mm^{3}/day$ | $k_{1}$ | $5\times 10^{-2}$ | $cell/mm^{3}$
$D_{3}$ | $3\times 10^{-1}$ | $mm^{3}/day$ | $\mu_{1}$ | $0.1$ | $day^{-1}$
$\mu_{2}$ | $3.3\times 10^{-2}$ | $day^{-1}$ | $q_{4}$ | $0.14$ | relative concentration
$\mu_{3}$ | $7\times 10^{-2}$ | $day^{-1}$ | $q_{5}$ | $0.6$ | relative concentration
$\gamma_{3}$ | $8.2\times 10^{-2}$ | $mm^{3}/cell.day$ | $q_{6}$ | $0.2$ | relative concentration
$\lambda_{2}$ | $5.98\times 10^{-3}$ | $mm^{3}/cell.day$ | $q_{7}$ | $0.15$ | relative concentration
$\lambda_{3}$ | $5.98\times 10^{-2}$ | $mm^{3}/cell.day$ | $\eta_{45}$ | $560$ | relative concentration
$\lambda_{13}$ | $7.14\times 10^{-2}$ | $mm^{6}/cell^{2}.day$ | $\eta_{47}$ | $17.4$ | relative concentration
$\lambda_{12}$ | $1.66\times 10^{-3}$ | $mm^{6}/cell^{2}.day$ | $\eta_{57}$ | $34.8$ | relative concentration
$\alpha_{1}$ | $4\times 10^{0}$ | $day^{-1}$ | $\eta_{55}$ | $560$ | relative concentration
$\alpha_{2}$ | $10^{-3}$ | $day^{-1}$ | $\eta_{54}$ | $185$ | relative concentration
$k_{2}$ | $2\times 10^{-1}$ | $day^{-1}$ | $\eta_{64}$ | $185$ | relative concentration
$k_{3}$ | $4.64$ | $day^{-1}$ | $\eta_{67}$ | $17.4$ | relative concentration
$k_{4}$ | $0.464$ | $day^{-1}$ | $\eta_{75}$ | $560$ | relative concentration
$k_{5}$ | $1.1$ | $day^{-1}$ | $H_{y_{4}}^{d}(y_{7})$ | $3.0$ | dimensionless
$k_{6}$ | $8.65$ | $hr^{-1}$ | $H_{y_{4}}^{d}(y_{5})$ | $2.0$ | dimensionless
$k_{7}$ | $1.5$ | $\frac{rel.cyt.conc.}{day\text{\,}of\text{\,}cell}$ | $H_{y_{5}}(y_{5})$ | $1.0$ | dimensionless
$k_{8}$ | $10^{-2}$ | $\frac{rel.cyt.conc.}{day\text{\,}of\text{\,}cell}$ | $H_{y_{7}}^{u}(y_{5})$ | $3.68$ | dimensionless
$k_{9}$ | $8.1\times 10^{-1}$ | $\frac{rel.cyt.conc.}{day\text{\,}of\text{\,}cell}$ | $H_{y_{5}}^{u}(y_{4})$ | $2.0$ | dimensionless
$k_{10}$ | $5.6\times 10^{-2}$ | $\frac{rel.cyt.conc.}{day\text{\,}of\text{\,}cell}$ | $H_{y_{5}}^{d}(y_{7})$ | $4.0$ | dimensionless
$k_{11}$ | $5.6\times 10^{-1}$ | $\frac{rel.cyt.conc.}{day\text{\,}of\text{\,}cell}$ | $H_{y_{6}}^{u}(y_{4})$ | $3.0$ | dimensionless
$k_{12}$ | $1.9\times 10^{-1}$ | $\frac{rel.cyt.conc.}{day\text{\,}of\text{\,}cell}$ | $H_{y_{6}}^{d}(y_{7})$ | $1.5$ | dimensionless
$k_{13}$ | $1.91\times 10^{-2}$ | $\frac{rel.cyt.conc.}{day\text{\,}of\text{\,}cell}$ |  |  |

where rel. cyt. conc.:=relative cytokine concentration.

###  2.4 Statistical approach

A regression model can be used to compare the numerical results and the
experimental data. The analysis should consider computed relative
concentrations of the pathogen S. aureus $y_{1}(\cdot)$, resting macrophages
$y_{2}(\cdot)$, activated macrophages $y_{3}(\cdot)$, tumor necrosis factor
alpha $y_{4}(\cdot)$, interleukin 6 $y_{5}(\cdot)$, interleukin 8
$y_{6}(\cdot)$ and interleukin 10 $y_{7}(\cdot)$, together with the
experimental concentrations for cellular and cytokines induced by either gram-
positive bacteria, gram-negative bacteria, lipoteichoic acid or peptidoglycan.
Furthermore, for either cellular or cytokine the regression model should be
linear:

$z=ay+b,$

where $a$ denotes the slope of each linear least square regression and $b$ is
the $z$-intercept.

##  3 Development of the second-order explicit predictor-corrector numerical
scheme

In this section, we develop a new second-order explicit predictor-corrector
numerical technique in a computed solution of the cellular-cytokine model
described by the nonlinear unsteady equations $(\ref{1})$-$(\ref{8})$.

Let $N$ be a positive integer and $\sigma:=\Delta t=\frac{T}{N}$, be the step
size. Set $t_{k}=k\sigma$, $t_{k+\frac{1}{2}}=\frac{t_{k}+t_{k+1}}{2}$, for
$k=0,1,2,...,N$, and
$\mathcal{T}_{\sigma}=\\{t_{k},\text{\,\,}k=0,1,...,N\\}$, be a regular
partition of $[0,\text{\,}T]$. For the convenience of writing, we set
$y_{r}(t):=y_{r}(x,t)$, for $r=4,...,7$, $y_{j}(x,t_{k})=y_{j}^{k}(x)$, for
$(x,t_{k})\in\Omega\times\mathcal{T}_{\sigma}$, $j=1,2,...,7$, and
$y^{k}(x)=[y_{1}^{k}(x),...,y_{7}^{k}(x)]\in\mathbb{R}^{7}$. The space of mesh
functions defined on $\Omega\times\mathcal{T}_{\sigma}$ is given by
$\mathcal{F}_{\sigma}=\\{y^{k}(x),\text{\,\,}0\leq k\leq
N,\text{\,}x\in\Omega\\}$. We introduce the vector functions $y$ and $F$
defined from $\Omega\times[0,\text{\,}T]$ to $\mathbb{R}^{7}$ and
$\mathbb{R}^{7}\times\Omega\times[0,\text{\,}T]$ to $\mathbb{R}^{7}$,
respectively, by

$y(x,t)=[y_{1}(x,t),...,y_{7}(x,t)]\text{\,\,\,\,\,and\,\,\,\,\,}F(y(x,t))=[F_{1}(y(x,t)),...,F_{7}(y(x,t))],$
(10)

where

$F_{1}(y(x,t))=\beta_{1}(1-k_{1}^{-1}y_{1})y_{1}-\mu_{1}y_{1}-\lambda_{2}y_{1}y_{2}-\lambda_{3}y_{1}y_{3}$,
---
$F_{2}(y(x,t))=\mu_{2}(1-y_{2m}^{-1}y_{2})y_{2}-[\gamma_{3}+k_{6}H_{c}^{u}(y_{4})H_{c}^{d}(y_{7})]y_{1}y_{2}$,
$F_{3}(y(x,t))=[\gamma_{3}+k_{6}H_{c}^{u}(y_{4})H_{c}^{d}(y_{7})]y_{1}y_{2}-\mu_{3}y_{3}$,
$F_{4}(y(x,t))=k_{7}H_{c}^{d}(y_{5})H_{c}^{d}(y_{7})y_{3}-k_{2}(y_{4}-q_{4})$,
$F_{5}(y(x,t))=[k_{8}+k_{9}H_{c}^{u}(y_{4})]H_{c}^{d}(y_{5})H_{c}^{d}(y_{7})]y_{3}-k_{3}(y_{5}-q_{5})$,
$F_{6}(y(x,t))=[k_{10}+k_{11}H_{c}^{u}(y_{4})]H_{c}^{d}(y_{7})y_{3}-k_{4}(y_{6}-q_{6})$,
$F_{7}(y(x,t))=[k_{12}+k_{13}H_{c}^{u}(y_{5})]y_{3}-k_{5}(y_{7}-q_{7})$.
(11)

Utilizing equations $(\ref{10})$-$(\ref{12})$, it is not hard to observe that
the system of nonlinear differential equations $(\ref{1})$-$(\ref{8})$ can be
expressed as

$\frac{\partial y}{\partial t}(x,t)=F(y(x,t)),$ (12)

where $\frac{\partial y}{\partial t}=[\frac{\partial y_{1}}{\partial
t},...,\frac{\partial y_{7}}{\partial t}]$. Subjects to initial-boundary
conditions

$y_{j}(x,0)=\overline{y}_{j},\text{\,\,\,for\,\,\,}j=1,...,7,\text{\,\,\,on\,\,\,\,}\overline{\Omega},$
(13) $\frac{\partial y_{1}}{\partial
t}(x,t)=0,\text{\,\,\,on\,\,\,\,}\partial\Omega\times[0,\text{\,}T].$ (14)

The integration of $(\ref{13})$ over the interval
$(t_{k},\text{\,}t_{k+\frac{1}{2}})$ provides

$y(x,t_{k+\frac{1}{2}})=y(x,t_{k})+\int_{t_{k}}^{t_{k+\frac{1}{2}}}F(y(x,t))dt,$

which is equivalent to

$y^{k+\frac{1}{2}}(x)=y^{k}(x)+\int_{t_{k}}^{t_{k+\frac{1}{2}}}F(y(x,t))dt.$
(15)

Let $P_{j}^{(1)}(x,t)$ be the first-order polynomial approximating the
functions $F_{j}(y(x,t))$ at the mesh points $(t_{m},F_{j}(y^{m}(x)))$, for
$m\in\\{k,k+\frac{1}{2}\\}$, thus

$F_{j}(y(x,t))=P_{j}^{(1)}(x,t)+E_{j}^{(1)}(y(x,t)),$ (16)

where

$P_{j}^{(1)}(x,t)=\frac{t-t_{k+\frac{1}{2}}}{t_{k}-t_{k+\frac{1}{2}}}F_{j}(y^{k}(x))+\frac{t-t_{k}}{t_{k+\frac{1}{2}}-t_{k}}F_{j}(y^{k+\frac{1}{2}}(x))=\frac{\sigma}{2}[F_{j}(y^{k+\frac{1}{2}}(x))-F_{j}(y^{k}(x))]+$
$t_{k+\frac{1}{2}}F_{j}(y^{k}(x))-t_{k}F_{j}(y^{k+\frac{1}{2}}(x)),$ (17)

and the associated error $E_{j}$ is defined as

$E_{j}^{(1)}(y(x,t))=\frac{1}{2}(t-t_{k})(t-t_{k+\frac{1}{2}})\frac{\partial^{2}F_{j}}{\partial
t^{2}}(y(x,t_{\epsilon}^{1})),$ (18)

where $t_{\epsilon}^{1}$ is between the minimum and maximum of
$t_{k+\frac{1}{2}}$, $t_{k}$ and $t$.

Integrating both sides of equation $(\ref{18})$ over the interval
$[t_{k},\text{\,}t_{k+\frac{1}{2}}]$, results in

$\int_{t_{k}}^{t_{k+\frac{1}{2}}}P_{j}^{(1)}(x,t)dt=\frac{1}{\sigma}[F_{j}(y^{k+\frac{1}{2}}(x))-F_{j}(y^{k}(x))](t_{k+\frac{1}{2}}^{2}-t_{k}^{2})+\frac{2}{\sigma}[t_{k+\frac{1}{2}}F_{j}(y^{k}(x))-t_{k}F_{j}(y^{k+\frac{1}{2}}(x))](t_{k+\frac{1}{2}}-t_{k}).$
(19)

Since $t_{m}=m\sigma$, so
$(t_{k+\frac{1}{2}}^{2}-t_{k}^{2})=(t_{k+\frac{1}{2}}-t_{k})(t_{k+\frac{1}{2}}+t_{k})=(2k+\frac{1}{2})\frac{\sigma^{2}}{2}$
and
$[t_{k+\frac{1}{2}}F_{j}(y^{k}(x))-t_{k}F_{j}(y^{k+\frac{1}{2}}(x))](t_{k+\frac{1}{2}}-t_{k})=\frac{\sigma^{2}}{2}[(k+\frac{1}{2})F_{j}(y^{k}(x))-kF_{j}(y^{k+\frac{1}{2}}(x))]$.
These facts along with equation $(\ref{20})$ give

$\int_{t_{k}}^{t_{k+\frac{1}{2}}}P_{j}^{(1)}(x,t)dt=\frac{\sigma}{4}[F_{j}(y^{k}(x))+F_{j}(y^{k+\frac{1}{2}}(x))].$
(20)

In addition, the integration of equation $(\ref{19})$ on the interval
$[t_{k},\text{\,}t_{k+\frac{1}{2}}]$, yields

$\int_{t_{k}}^{t_{k+\frac{1}{2}}}E_{j}^{(1)}(y(x,t))dt=\frac{1}{2}\int_{t_{k}}^{t_{k+\frac{1}{2}}}(t-t_{k})(t-t_{k+\frac{1}{2}})\frac{\partial^{2}F_{j}}{\partial
t^{2}}(y(x,t))dt.$

The absolute value in both sides of this equation implies

$\left|\int_{t_{k}}^{t_{k+\frac{1}{2}}}E_{j}^{(1)}(y(x,t))dt\right|\leq\frac{1}{2}\int_{t_{k}}^{t_{k+\frac{1}{2}}}|t-t_{k}||t-t_{k+\frac{1}{2}}|\left|\frac{\partial^{2}F_{j}}{\partial
t^{2}}(y(x,t))\right|dt\leq\frac{\sigma^{3}}{16}\underset{0\leq t\leq
T}{\sup}\left|\frac{\partial^{2}F_{j}}{\partial t^{2}}(y(x,t))\right|.$ (21)

But, utilizing equations $(\ref{1})$-$(\ref{8})$, it follows that the
functions $y_{j}(\cdot)$ have continuous first-order partial derivatives on
the bounded domain $\overline{\Omega}\times[0,\text{\,}T]$. This fact together
with equations $(\ref{2})$ and $(\ref{12})$ show that the functions $F_{j}$
have continuous second-order partial derivatives. Thus, $\underset{0\leq t\leq
T}{\sup}\left|\frac{\partial^{2}F_{j}}{\partial t^{2}}(y(x,t))\right|\leq
C_{0}$, where $C_{0}>0$, is a constant independent of $\sigma$. Using this and
estimate $(\ref{22})$, one can write

$\int_{t_{k}}^{t_{k+\frac{1}{2}}}E_{j}^{(1)}(y(x,t))dt=O(\sigma^{3}),\text{\,\,\,for\,\,\,\,}j=1,2,...,7.$
(22)

Integrating equation $(\ref{17})$ on $[t_{k},\text{\,}t_{k+\frac{1}{2}}]$ and
utilizing equations $(\ref{21})$ and $(\ref{23})$, to get

$\int_{t_{k}}^{t_{k+\frac{1}{2}}}F_{j}(y(x,t))dt=\frac{\sigma}{4}[F_{j}(y^{k}(x))+F_{j}(y^{k+\frac{1}{2}}(x))]+O(\sigma^{3}),\text{\,\,\,for\,\,\,\,}j=1,2,...,7.$

Substituting this into equation $(\ref{16})$, to obtain

$y^{k+\frac{1}{2}}(x)=y^{k}(x)+\frac{\sigma}{4}[F(y^{k}(x))+F(y^{k+\frac{1}{2}}(x))]+\mathcal{O}(\sigma^{3}),$
(23)

where $\mathcal{O}(k^{3})=(O(\sigma^{3}),...,O(\sigma^{3}))$. We should
approximate the term $F_{j}(y^{k}(x))+F_{j}(y^{k+\frac{1}{2}}(x))$, by the sum
$c_{1}F_{j}(y^{k}(x))+c_{2}F_{j}[y^{k}(x)+\sigma pF(y^{k}(x))]$, in which the
scalars $c_{1}$, $c_{2}$ and $p$ are chosen so that the following equation
holds

$\frac{2}{\sigma}\left(y^{k+\frac{1}{2}}(x)-y^{k}(x)\right)-\frac{1}{2}[F(y^{k}(x))+F(y^{k+\frac{1}{2}}(x))]=\mathcal{O}(\sigma^{2}).$

Expanding the Taylor series for $y_{j}(x,\cdot)$ and $F_{j}$ about the points
$t_{k}$ and $y^{k}(x)$, respectively, with steplength $\frac{\sigma}{2}$
utilizing forward difference formulation, simple computations result in

$\frac{2}{\sigma}\left(y_{j}^{k+\frac{1}{2}}(x)-y_{j}^{k}(x)\right)-\frac{1}{2}[F_{j}(y^{k}(x))+F_{j}(y^{k+\frac{1}{2}}(x))]=(1-\frac{c_{1}+c_{2}}{2})F_{j}(y^{k}(x))+\frac{\sigma}{4}(1-2c_{2}p)\frac{\partial
F_{j}}{\partial t}(y^{k}(x))+O(\sigma^{2}).$

The right side of this equations equals $O(\sigma^{2})$, if and only if,
$c_{1}+c_{2}=2$ and $c_{2}p=\frac{1}{2}$. For instance, we can take
$c_{1}=\frac{3}{2}$, $c_{2}=\frac{1}{2}$ and $p=1$. So,

$F_{j}(y^{k}(x))+F_{j}(y^{k+\frac{1}{2}}(x))=\frac{3}{2}F_{j}(y^{k}(x))+\frac{1}{2}F_{j}[y^{k}(x)+\sigma
F(y^{k}(x))],\text{\,\,\,for\,\,\,}j=1,2,...,7.$

This is equivalent to

$F(y^{k}(x))+F(y^{k+\frac{1}{2}}(x))=\frac{3}{2}F(y^{k}(x))+\frac{1}{2}F[y^{k}(x)+\sigma
F(y^{k}(x))].$ (24)

Plugging equations $(\ref{24})$ and $(\ref{24a})$, we obtain

$y^{k+\frac{1}{2}}(x)=y^{k}(x)+\frac{\sigma}{8}\left[3F(y^{k}(x))+F\left(y^{k}(x)+\sigma
F(y^{k}(x))\right)\right]+\mathcal{O}(\sigma^{3}).$

Truncating the error term $\mathcal{O}(\sigma^{3})$ and replacing the
analytical solution $y(x,\cdot)$ with the computed one $Y(x,\cdot)$, this
gives

$Y^{k+\frac{1}{2}}(x)=Y^{k}(x)+\frac{\sigma}{8}\left[3F(Y^{k}(x))+F\left(Y^{k}(x)+\sigma
F(Y^{k}(x))\right)\right].$ (25)

Equation $(\ref{25})$ represents the predictor phase of the desired algorithm.

In a similar manner, denoting by $P_{j}^{(2)}(x,t)$ be the linear polynomial
in $t$ interpolating the function $F_{j}(y(x,t))$ at the mesh points
$(t_{k+\frac{1}{2}},F_{j}(y^{k+\frac{1}{2}}(x)))$ and
$(t_{k+1},F_{j}(y^{k+1}(x)))$, and $E_{j}^{(2)}(y(x,t))$ be the corresponding
error. Replacing $F(y(x,t))$ with $P^{(2)}(x,t)+E^{(2)}(y(x,t))$ into equation
$(\ref{13})$ and integration over the interval
$[\text{\,}t_{k+\frac{1}{2}},t_{k+1}]$, where
$P^{(2)}(x,t)=[P_{1}^{(2)}(x,t),...,P_{7}^{(2)}(x,t)]$ and
$E^{(2)}(y(x,t))=[E_{1}^{(2)}(y(x,t)),...,E_{7}^{(2)}(y(x,t))]$,
straightforward computations provide

$y^{k+1}(x)=y^{k+\frac{1}{2}}(x)+\frac{\sigma}{4}[F(y^{k+\frac{1}{2}}(x))+F(y^{k+1}(x))]+\mathcal{O}(\sigma^{3}).$
(26)

By the application of Taylor series, the quantity
$F(y^{k+\frac{1}{2}}(x))+F(y^{k+1}(x))$ can be approximated as

$F_{j}(y^{k+\frac{1}{2}}(x))+F_{j}(y^{k+1}(x))=F_{j}(y^{k+\frac{1}{2}}(x))+F_{j}[y^{k+\frac{1}{2}}(x)+\frac{\sigma}{2}F(y^{k+\frac{1}{2}}(x))],\text{\,\,\,for\,\,\,}j=1,2,...,7.$

Thus,

$F(y^{k+\frac{1}{2}}(x))+F(y^{k+1}(x))=F(y^{k+\frac{1}{2}}(x))+F[y^{k+\frac{1}{2}}(x)+\frac{\sigma}{2}F(y^{k+\frac{1}{2}}(x))].$
(27)

Combining equations $(\ref{26})$ and $(\ref{27})$, this yields

$y^{k+1}(x)=y^{k+\frac{1}{2}}(x)+\frac{\sigma}{4}\left[F(y^{k+\frac{1}{2}}(x))+F\left(y^{k+\frac{1}{2}}(x)+\frac{\sigma}{2}F(y^{k+\frac{1}{2}}(x))\right)\right]+\mathcal{O}(\sigma^{3}).$
(28)

Tracking the infinitesimal term $\mathcal{O}(\sigma^{3})$ and replacing the
exact solution $y(x,\cdot)$ with the approximate one $Y(x,\cdot)$, to get the
corrector step of the proposed approach

$Y^{k+1}(x)=Y^{k+\frac{1}{2}}(x)+\frac{\sigma}{4}\left[F(Y^{k+\frac{1}{2}}(x))+F\left(Y^{k+\frac{1}{2}}(x)+\frac{\sigma}{2}F(Y^{k+\frac{1}{2}}(x))\right)\right].$
(29)

But equation $(\ref{15})$ suggests that $\frac{\partial y_{1}}{\partial
t}(x,t)=0$, on $\partial\Omega\times[0,T]$. This indicates that $y_{1}(x,t)$
is constant on $\partial\Omega\times[0,T]$. Since,
$y_{1}(x,0)=\overline{y}_{1}$, so $y_{1}^{k}(x)=\overline{y}_{1}$, for
$k=0,1,...,N$, $x\in\partial\Omega$. Plugging this equation together with
equations $(\ref{14})$, $(\ref{25})$ and $(\ref{29})$, to obtain the developed
second-order explicit predictor-corrector numerical technique for solving the
mathematical model of the dynamic of coupled cellular-cytokine problem
$(\ref{1})$-$(\ref{8})$. That is, for $k=0,1,...,N-1$, and every $x\in\Omega$,

$Y^{k+\frac{1}{2}}(x)=Y^{k}(x)+\frac{\sigma}{8}\left[3F(Y^{k}(x))+F\left(Y^{k}(x)+\sigma
F(Y^{k}(x))\right)\right],$ (30)
$Y^{k+1}(x)=Y^{k+\frac{1}{2}}(x)+\frac{\sigma}{4}\left[F(Y^{k+\frac{1}{2}}(x))+F\left(Y^{k+\frac{1}{2}}(x)+\frac{\sigma}{2}F(Y^{k+\frac{1}{2}}(x))\right)\right],$
(31)

with initial and boundary conditions

$Y_{j}^{0}=\overline{y}_{j},\text{\,\,\,for\,\,\,}j=1,2,...,7,\text{\,\,\,on\,\,\,}\Omega,\text{\,\,\,\,and\,\,\,\,}Y^{k}_{1}(x)=\overline{y}_{1},\text{\,\,\,for\,\,\,}k=0,1,2,...,N,\text{\,\,\,}x\in\partial\Omega.$
(32)

We remind that $Y=[Y_{1},...,Y_{7}]$ and $F=[F_{1},...,F_{7}]$.

######  Theorem 3.1.

The numerical approach $(\ref{s1})$-$(\ref{s2})$ is zero-stable and second-
order accurate for any values of the initial-boundary conditions given by
$(\ref{s3})$.

###### Proof.

A combination of approximations $(\ref{s1})$ and$(\ref{s2})$ results in

$Y^{k+1}(x)-Y^{k}(x)=\frac{\sigma}{8}\left[3F(Y^{k}(x))+F\left(Y^{k}(x)+\sigma
F(Y^{k}(x))\right)\right]+$
$\frac{\sigma}{4}\left[F(Y^{k+\frac{1}{2}}(x))+F\left(Y^{k+\frac{1}{2}}(x)+\frac{\sigma}{2}F(Y^{k+\frac{1}{2}}(x))\right)\right].$

It follows from this equation that the first characteristic polynomial in term
of $\lambda^{\frac{1}{2}}$ of the developed explicit predictor-corrector
numerical scheme is defined as

$P_{2}(\lambda^{\frac{1}{2}})=(\lambda^{\frac{1}{2}})^{2}-1=(\lambda^{\frac{1}{2}}-1)(\lambda^{\frac{1}{1}}+1).$

It is not difficult to see that the roots of this polynomial are:
$\lambda_{1}^{\frac{1}{2}}=1$ and $\lambda_{2}^{\frac{1}{2}}=-1$. Since the
two roots are simple and lie in the closed unit disc, it comes from the
definition of zero-stability [19, 18] that the proposed scheme
$(\ref{s1})$-$(\ref{s3})$ is zero-stable. Furthermore, since the truncation
errors in both approximations $(\ref{25})$ and $(\ref{29})$ equal
$\mathcal{O}(\sigma^{3})$, the proof of the second-order accuracy for the
constructed technique $(\ref{s1})$-$(\ref{s3})$ is similar to that established
in [49]. ∎

##  4 Numerical experiments and Discussions

This section presents some numerical simulations to demonstrate the utility
and effectiveness of the new second-order explicit predictor-corrector scheme
$(\ref{s1})$-$(\ref{s3})$ in approximate solutions of the mathematical model
on the dynamic of cytokine levels and human immune cell activation in response
to the pathogen S. aureus $(\ref{1})$-$(\ref{8})$. Since the proposed
technique is an explicit predictor-corrector approach and second-order
accuracy, it is faster and more efficient than a wide set of statistical and
numerical methods discussed in the literature for solving general systems of
mixed ODEs/PDEs modeling real-world problems [26, 28, 65] and references
therein. Some data used in the simulations are taken from [71]. Both tables
and graphs (Tables $3$-$4$ and Figures 2-3) suggest that the approximate
solutions (macrophages and cytokine expressions) increase until a maximum
concentration, so called ”peaks”, after approximately $t$ hours, and become
constant when the antigens are eliminated by the immune response. The
activated macrophages behave as the pathogens presenting cells and they move
to the nearest lymph node where the specific response is triggered. The peaks
of TNF$\alpha$, IL6, IL8 and IL10 (see Tables $3$-$4$), agree with the model
results (Figure 3). Figure 2 considers the relationship between S. aureus
$(Y_{1}(\cdot))$, resting macrophages $(Y_{2}(\cdot))$ and activated
macrophages $(Y_{3}(\cdot))$ whereas Figures 3 indicates how the changes in
activated macrophages and resting macrophages impact TNF$\alpha$
$(Y_{4}(\cdot))$, IL6 $(Y_{5}(\cdot))$, IL8 $(Y_{6}(\cdot))$ and IL10
$(Y_{7}(\cdot))$. Additionally, the simulated results do not include some
factors such as: humoral immune response effects, neutrophil flux or
complement response. The initial concentrations of the pathogens, macrophages
and cytokines are chosen to match well those reported in the ex vivo studies
[78, 75]. Furthermore, the numerical analysis shows that the constructed
predictor-corrector explicit method $(\ref{s1})$-$(\ref{s3})$ accurately
predicts the dynamic of human immune cell activation and cytokine levels in
response to S. aureus. Thus, the new algorithm can be considered as a robust
tool to predict ex vivo and in vivo experimental data induced by a given
antigen (Tables $3$-$4$). The graphs (Figures 2-3) also indicate that the body
response varies based on the damage inflicted by the bacteria S. aureus.
Because of the activation of cytokines in the human immune response, the
numerical experiments are performed assuming low initial concentration of
gram-positive bacteria S. aureus as discussed in [71]. However, high
expressions of the pathogen S. aureus due to the tissue damage associated with
endotoxicity of gram-positive bacteria provoke a fast increase in cellular and
cytokine responses (see Figure 4). This case is not analyzed in this work and
should be considered as the topic of our future investigations. Furthermore,
high concentrations of the antigens are not sufficiently discussed in the
literature because of the lack of data to validate the increased expressions
[71].

Table 3. Concentrations of S. aureus, macrophages and cytokines (in
cells$/mm^{3}$) at point $x=(2\times 10^{-3},\text{\,}10^{-3},\text{\,}2\times
10^{-3})$.

time (hrs) | S. aureus | rest. macroph | act. macroph |  TNF$\alpha$ | IL6 | IL8 | IL10
---|---|---|---|---|---|---|---
0 | $2.0\times 10^{-1}$ | $10^{-2}$ | $5.0\times 10^{-3}$ | $0.0\times 10^{0}$ | $0.0\times 10^{0}$ | $0.0\times 10^{0}$ | $0.0\times 10^{0}$
1 | $2.22\times 10^{-2}$ | $1.29\times 10^{-2}$ | $4.80\times 10^{-3}$ | $2.47\times 10^{-2}$ | $5.93\times 10^{-1}$ | $7.20\times 10^{-2}$ | $9.84\times 10^{-2}$
3 | $1.91\times 10^{-2}$ | $2.18\times 10^{-2}$ | $4.40\times 10^{-3}$ | $6.33\times 10^{-2}$ | $6.00\times 10^{-1}$ | $1.48\times 10^{-1}$ | $1.45\times 10^{-1}$
6 | $1.90\times 10^{-2}$ | $4.19\times 10^{-2}$ | $4.20\times 10^{-3}$ | $9.88\times 10^{-2}$ | $6.00\times 10^{-1}$ | $1.88\times 10^{-1}$ | $1.505\times 10^{-1}$
9 | $1.90\times 10^{-2}$ | $6.42\times 10^{-2}$ | $4.40\times 10^{-3}$ | $1.18\times 10^{-1}$ | $6.00\times 10^{-1}$ | $1.97\times 10^{-1}$ | $1.508\times 10^{-1}$
12 | $1.90\times 10^{-2}$ | $8.09\times 10^{-2}$ | $5.00\times 10^{-3}$ | $1.29\times 10^{-1}$ | $6.00\times 10^{-1}$ | $1.97\times 10^{-1}$ | $1.509\times 10^{-1}$
15 | $1.90\times 10^{-2}$ | $8.98\times 10^{-2}$ | $5.70\times 10^{-3}$ | $1.35\times 10^{-1}$ | $6.00\times 10^{-1}$ | $2.00\times 10^{-1}$ | $1.510\times 10^{-1}$
18 | $1.90\times 10^{-2}$ | $9.39\times 10^{-2}$ | $6.40\times 10^{-3}$ | $1.38\times 10^{-1}$ | $6.00\times 10^{-1}$ | $2.004\times 10^{-1}$ | $1.511\times 10^{-1}$
21 | $1.90\times 10^{-2}$ | $9.55\times 10^{-2}$ | $7.00\times 10^{-3}$ | $1.40\times 10^{-1}$ | $6.00\times 10^{-1}$ | $2.005\times 10^{-1}$ | $1.512\times 10^{-1}$
24 | $1.90\times 10^{-2}$ | $9.61\times 10^{-2}$ | $7.50\times 10^{-3}$ | $1.41\times 10^{-1}$ | $6.00\times 10^{-1}$ | $2.005\times 10^{-1}$ | $1.513\times 10^{-1}$
peak | — | $9.61\times 10^{-2}$ | $7.50\times 10^{-3}$ | $1.411\times 10^{-1}$ | $6.00\times 10^{-1}$ | $2.005\times 10^{-1}$ | $1.513\times 10^{-1}$

Table 4. Concentrations of S. aureus, macrophages and cytokines (in
cells$/mm^{3}$) at point $x=(2\times 10^{-3},\text{\,}3\times
10^{-3},\text{\,}6\times 10^{-3})$

time (hrs) | S. aureus | rest. macroph | act. macroph |  TNF$\alpha$ | IL6 | IL8 | IL10
---|---|---|---|---|---|---|---
0 | $4.68\times 10^{-1}$ | $2.33\times 10^{-2}$ | $1.17\times 10^{-2}$ | $0.0\times 10^{0}$ | $0.0\times 10^{0}$ | $0.0\times 10^{0}$ | $0.0\times 10^{0}$
1 | $2.24\times 10^{-2}$ | $2.87\times 10^{-2}$ | $1.13\times 10^{-2}$ | $2.51\times 10^{-2}$ | $5.93\times 10^{-1}$ | $7.22\times 10^{-2}$ | $9.92\times 10^{-2}$
3 | $1.91\times 10^{-2}$ | $4.29\times 10^{-2}$ | $1.03\times 10^{-2}$ | $6.43\times 10^{-2}$ | $6.00\times 10^{-1}$ | $1.50\times 10^{-1}$ | $1.46\times 10^{-1}$
6 | $1.90\times 10^{-2}$ | $6.54\times 10^{-2}$ | $9.40\times 10^{-3}$ | $1.00\times 10^{-1}$ | $6.00\times 10^{-1}$ | $1.88\times 10^{-1}$ | $1.515\times 10^{-1}$
9 | $1.90\times 10^{-2}$ | $8.15\times 10^{-2}$ | $9.10\times 10^{-3}$ | $1.20\times 10^{-1}$ | $6.00\times 10^{-1}$ | $1.98\times 10^{-1}$ | $1.516\times 10^{-1}$
12 | $1.90\times 10^{-2}$ | $9.02\times 10^{-2}$ | $9.00\times 10^{-3}$ | $1.30\times 10^{-1}$ | $6.00\times 10^{-1}$ | $2.00\times 10^{-1}$ | $1.516\times 10^{-1}$
15 | $1.90\times 10^{-2}$ | $9.40\times 10^{-2}$ | $9.10\times 10^{-3}$ | $1.36\times 10^{-1}$ | $6.00\times 10^{-1}$ | $2.01\times 10^{-1}$ | $1.516\times 10^{-1}$
18 | $1.90\times 10^{-2}$ | $9.55\times 10^{-2}$ | $9.20\times 10^{-3}$ | $1.39\times 10^{-1}$ | $6.00\times 10^{-1}$ | $2.01\times 10^{-1}$ | $1.516\times 10^{-1}$
21 | $1.90\times 10^{-2}$ | $9.62\times 10^{-2}$ | $9.30\times 10^{-3}$ | $1.41\times 10^{-1}$ | $6.00\times 10^{-1}$ | $2.01\times 10^{-1}$ | $1.516\times 10^{-1}$
24 | $1.90\times 10^{-2}$ | $9.64\times 10^{-2}$ | $9.40\times 10^{-3}$ | $1.42\times 10^{-1}$ | $6.00\times 10^{-1}$ | $2.01\times 10^{-1}$ | $1.516\times 10^{-1}$
peak | — | $9.64\times 10^{-2}$ | $1.17\times 10^{-2}$ | $1.419\times 10^{-1}$ | $6.00\times 10^{-1}$ | $2.01\times 10^{-1}$ | $1.516\times 10^{-1}$

##  5 General conclusion and future works

In this paper, we have developed a second-order explicit predictor-corrector
numerical approach in simulated results of a coupled cellular-cytokine model
defined by the system of nonlinear equations $(\ref{1})$-$(\ref{8})$. The
theory has suggested that the proposed numerical technique is second-order
convergent and stable for any value of the initial datum (Theorem 3.1).
Additionally, the study indicates that the new algorithm is faster and more
efficient than a broad range of statistical techniques and numerical schemes
discussed in the literature for solving such systems of nonlinear differential
equations [54, 28, 26, 65, 10]. Furthermore, the graphs show that the cellular
model predicts the relationship between the S. aureus bacteria ($Y_{1}$),
resting macrophages ($Y_{2}$) and activated ones ($Y_{3}$) whereas the
cytokine model analyzes how the changes in the resting and activated
macrophages influence the propagation of tumor necrosis factor alpha
($Y_{4}$), interleukin 6 ($Y_{5}$), interleukin 8 ($Y_{6}$) and interleukin 10
($Y_{7}$). The developed numerical approach is used to investigating and
predicting the dynamic of cytokine levels and human immune cell activation in
response to the pathogen S. aureus. This suggests that the proposed technique
$(\ref{s1})$-$(\ref{s3})$ can be observed as a robust tool to predicting in
vivo and ex vivo experimental data induced by a specific bacterium (see Tables
$3$-$4$). The considered bacteria neutrophil interactions are human-specific
and may impact the predicted cytokine average expressions. More accurate
simulated results should be obtained by incorporating the effect of
neutrophils in the ODE modeling the interleukin 8. Finally, the model defined
by the initial-boundary value problem $(\ref{1})$-$(\ref{8})$ does not
incorporate relevant complement proteins. Complement response of human immune
system to gram-positive bacterium (for instance, S. aureus) is fundamental in
the activation of chemoattractants for phagocytosis of the antigens.
Specifically, cytokine and complement responses have overlapping biological
effects on the human body under septic conditions. However, solve a coupling
of ordinary and partial differential equations by fast and efficient numerical
approaches provides a more realistic representation of the complex
relationships within the immune system and may serve as a good test of drugs
in silico. Our future works will develop a third-order explicit numerical
technique to predicting the dynamic of mixed cellular-cytokine model with
complement proteins.

### Ethical Approval

Not applicable.

### Availability of supporting data

Not applicable.

### Declaration of Interest Statement

The authors declare that they have no conflict of interests.

### Funding

Not applicable.

### Authors’ contributions

The whole work has been carried out by the three authors.

## References

  * [1] C. Alexander, E. T. Rietschel. ”Bacterial Lipopolysaccharides and Innate Immunity”, J. Endotoxin Res., $7(3)$ $(2001)$, $167$-$202$.
  * [2] L. M. B. Assas. ”New exact solutions for the Kawahara equation using Exp-function method”, J. Comput. Appl. Math., $233(2009)$, $97$-$102$.
  * [3] C. G. Atkins, K. Buckley, M. W. Blades, R. F. B. Turner. ”Raman Spectroscopy of Blood and Blood Components”, Appl. Spectrosc., $71(2017)$, $767$-$793$.
  * [4] J. A. Badwey, M. L. Karnovsky. ” Active Oxygen Species and the Functions of Phagocytic Leukocytes”, Annu. Rev. Biochem., $49(1)$ $(1980)$, $695$-$726$
  * [5] R. Brady, F. -I. Do, H. T. Tran, S. Janum, K. M$\phi$ller, S. Brix, et al. ”Personalized Mathematical Model Predicting Endotoxin-Induced Inflammatory Responses in Young Men”, Math. Model. Nat. Phenom., $13(42)$ $(2016)$, $1$-$20$.
  * [6] L. Caudill, F. Lynch. ”A Mathematical Model of the Inflammatory Response to Pathogen Challenge”, Bull. Math. Biol.”, $80(8)$ $(2018)$, $2242$-$2271$.
  * [7] D. S. Chabot-Richards, T. I. George. ”White blood cell counts: Reference methodology”, Clin. Lab. Med., $35(2015)$, $11$-$24$.
  * [8] N. Chaudhary, T. N. Q. Nguyen, D. Cullen, A. D. Meade, C. Wynne. ”Discrimination of immune cell activation using Raman micro-spectroscopy in an in-vitro $\&$ ex-vivo model”, Spectrochim. Acta Part A Mol. Biomol. Spectrosc., $248(2021)$, $119118$.
  * [9] M. Chen, N. McReynolds, E. C. Campbell, M. Mazilu, J. Barbosa, K. Dholakia, S. J. Powis. ”The Use of Wavelength Modulated Raman Spectroscopy in Label-Free Identification of T Lymphocyte Subsets, Natural Killer Cells and Dendritic Cells”, PLoS ONE, $10(2015)$, $e0125158$.
  * [10] X. Chen, W. Wang, G. Xie, R. Hontecillas, M. Verma, A. Leber, et al. ”Multi-Resolution Sensitivity Analysis of Model of Immune Response to Helicobacter Pylori Infection via Spatio-Temporal Metamodeling”, Front. Appl. Math. Stat. 5, $4(4)$ $(2019)$.
  * [11] J. Cole, J. Aberdein, J. Jubrail, D. H. Dockrell. ”The Role of Macrophages in the Innate Immune Response to Streptococcus Pneumoniae and Staphylococcus Aureus: Mechanisms and Contrasts”, Adv. Microb. Physiol., $65(2014)$, $125$-$202$.
  * [12] J. Crank. ”Mathematics of Diffusion”, (London: Oxford $\&$ Clarendon Press), $(1975)$.
  * [13] S. J. De Kimpe, M. Kengatharan, C. Thiemermann, J. R. Vane. ”The Cell Wall Components Peptidoglycan and Lipoteichoic Acid From Staphylococcus Aureus Act in Synergy to Cause Shock and Multiple Organ Failure”, Proc. Natl. Acad. Sci. U.S.A., $92(22)$ $(1995)$, $10359$-$10363$.
  * [14] R. de Waal Malefyt, J. Haanen, H. Spits, M. G. Roncarolo, A. te Velde, C. Figdor, et al. ”Interleukin 10 (IL-10) and Viral IL-10 Strongly Reduce Antigen-Specific Human T Cell Proliferation by Diminishing the Antigen-Presenting Capacity of Monocytes via Downregulation of Class II Major Histocompatibility Complex Expression”, J. Exp. Med., $174(4)$ $(1991)$, $915$-$924$.
  * [15] P. J. Dellar, R. Salmon. ”Shallow water equations with a complete coriolis force and topography”, Phys. Fluids, $17(10)$ $(2005)$, $106601$.
  * [16] M. V. DiLeo, J. A. Kellum, W. J. Federspiel. ”A Simple Mathematical Model of Cytokine Capture Using a Hemoadsorption Device”, Ann. BioMed. Eng., $37(1)$ $(2009)$, $222$-$229$.
  * [17] R. Eftimie, J. J. Gillard, D. A. Cantrell. ”Mathematical Models for Immunology: Current State of the Art and Future Research Directions”, Bull. Math. Biol., $78(10)$ $(2016)$, $2091$-$2134$.
  * [18] G. Dahlquist. ”A special stability problem for linear multistep methods”, BIT $3(1963)$, $27$-$43$.
  * [19] G. Dahlquist. ”Cnvergence and stability in the numerical integration of ordinary differential equations”, Math. Scand., $4(1956)$, $33$-$56$.
  * [20] H. Fan, D. L. Williams, B. Zingarelli, K. F. Breuel, G. Teti, G. E. Tempel, et al. ”Differential Regulation of Lipopolysaccharide and Gram-Positive Bacteria Induced Cytokine and Chemokine Production in Macrophages by Gai Proteins”, Immunology, $122(1)$ $(2007)$, $116$-$123$.
  * [21] W. A. Farahat, L. B. Wood, I. K. Zervantonakis, A. Schor, S. Ong, D. Neal, et al. ”Ensemble analysis of angiogenic growth in three-dimensional microfluidic cell cultures”, PLoS ONE, $7(2012)$, $e37333$.
  * [22] R. J. Feezor, C. Oberholzer, H. V. Baker, D. Novick, M. Rubinstein, L.L. Moldawer, et al. ”Molecular Characterization of the Acute Inflammatory Response to Infections With Gram-Negative Versus Gram-Positive Bacteria”, Infection Immun., $71(10)$ $(2003)$, $5803$-$5813$.
  * [23] G. Z. Feuerstein, T. Liu, F. C. Barone. ”Cytokines, Inflammation, and Brain Injury: Role of Tumor Necrosis Factor-Alpha”, Cerebrovasc Brain Metab. Rev., $6(4)$ $(1994)$, $341$-$360$.
  * [24] B. Fournier, D. J. Philpott. ”Recognition of Staphylococcus Aureus by the Innate Immune System”, Clin. Microbiol. Rev., $18(3)$ $(2005)$, $521$-$540$.
  * [25] J. Fu, J. Zheng, W. Fang, B. Wu. ”Effect of interleukin-6 on the growth of human lung cancer cell line”, Chin Med J., $111(1998)$, $265$-$8$.
  * [26] D. Gamerman, H. Lopes. ”Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference”, New York, NY: Chapman and Hall/CRC $(2006)$.
  * [27] D. Giri, M. Ozen, M. Ittmann. ”Interleukin-6 is an autocrine growth factor in human prostate cancer”, Am. J. Pathol., $159(2001)$, $2159$-$65$.
  * [28] W. Hastings. ”Monte carlo sampling methods using markov chains and their applications”, Biometrika, $57(1970)$, $97$-$109$.
  * [29] C. B. Hergott, A. M. Roche, E. Tamashiro, T. B. Clarke, A. G. Bailey, A. Laughlin, et al. ”Peptidoglycan From the Gut Microbiota Governs the Lifespan of Circulating Phagocytes at Homeostasis”, Blood, $127(20)$ $(2016)$, $2460$-$2471$.
  * [30] W. E. Holmes, J. Lee, W. J. Kuang, G. C. Rice, W. I Wood. ”Structure and Functional Expression of a Human Interleukin-8 Receptor”, Science, $253(5025)$ $(1991)$, $1278$-$1280$.
  * [31] W. K. E. Ip, N. Hoshi, D. S. Shouval, S. Snapper, R. Medzhitov. ”Anti-Inflammatory Effect of IL-10 Mediated by Metabolic Reprogramming of Macrophages”, Science, $356(6337)$ $(2017)$, $513$-$519$.
  * [32] M. Jambunathan. ”Some properties of beta and gamma distributions”, Ann. Math. Stat., $25(1954)$, $401$-$5$.
  * [33] J. Kim, P. Bentley, U. Aiklelin, J. Greesmith, G. Tedesco, J. Twycross. ”Immune System Approaches to Intrusion Detection-A Review”, Natural Computing. Natural Comput., $6(2007)$, $413$-$466$.
  * [34] A. E. Koch, P. J. Polverini, S. L. Kunkel SL, L. A. Harlow, L. A. DiPietro, V. M. Elner, S. G. Elner, R. M. Strieter. ”Interleukin-8 as a macrophage-derived mediator of angiogenesis”, Science, $258(1992)$, $1798$-$801$.
  * [35] E. Kovacs. ”Investigation of interleukin-6 (IL-6), soluble IL-6 receptor (sIL-6R) and soluble gp130 (sgp130) in sera of cancer patients”, Biomed Pharmacother, $55(2001)$, $391$-$6$.
  * [36] T. Kulahava, L. Golubewa, N. Vasilieu, O. Paddubskaya, M. Shuba. ”P-73-Neutrophils activation by carbon nanotubes functionalized with COOH-groups and conjugated with DNA”, Free Radic. Biol. Med., $120(2018)$, $S66$-$S67$.
  * [37] X. Liu, S. Yin, Y. Chen, Y. Wu, W. Zheng, H. Dong, et al. ”LPS-induced Proinflammatory Cytokine Expression in Human Airway Epithelial Cells and Macrophages via NF-kb, $STAT3$ or $AP$-$1$ Activation”, Mol. Med. Rep., $17(4)$ $(2018)$, $5484$-$5491$.
  * [38] K. M. Lowings, M. Christodoulides, J. A. Warner. ”Cytokine Response to Staphylococcus Aureus (SA) in Human Lung, Skin and Nasal Tissue”, J. Allergy Clin. Immunol., $123(2)$ $(2009)$,$S51$.
  * [39] M. Luca, S. Huang, J. E. Gershenwald, R. K. Singh, R. Reich, M. Bar-Eli. ”Expression of interleukin-8 by human melanoma cells up-regulates MMP-$2$ activity and increases tumor growth and metastasis”, Am. J. Pathol., $151(1997)$, $1105$-$13$.
  * [40] S. Managò, P. Mirabelli, M. Napolitano, G. Zito, A. C. De Luca. ”Raman detection and identification of normal and leukemic hematopoietic cells”, J. Biophotonics, $11(2018)$, $e201700265$.
  * [41] G. I. Marchuk. ”Mathematical Modelling of Immune Response in Infectious Diseases”, Netherlands: Springer, $(1997)$.
  * [42] M. Meier-Schellersheim, R. Varma, B. R. Angermann. ”Mechanistic Models of Cellular Signaling, Cytokine Crosstalk, and Cell-Cell Communication in Immunology”, Front. Immunol. 10, $2268(2268)$ $(2019)$.
  * [43] R. Mirzaei, R. Ranjbar, S. Karampoor, R. Goodarzi, H. Hasanvand. ”The Human Immune System Toward Staphylococcus Aureus”, Open Microbiol. J., $14(2020)$, $164$-$170$.
  * [44] W. Möbius, L. Laan. ”Physical and Mathematical Modeling in Experimental Papers”, Cell, $163(7)$ $(2015)$, $1577$-$1583$.
  * [45] M. Monguió-Tortajada, M. Franquesa, M. R. Sarrias, F. E. Borràs. ”Low Doses of LPS Exacerbate the Inflammatory Response and Trigger Death on TLR3-Primed Human Monocytes,” Cell Death Dis., $9(5)$ $(2018)$, $499$.
  * [46] K. W. Moore, R. de Waal Malefyt, R. L. Coffman, A. O’Garra. ”Interleukin-10 and the interleukin-10 receptor”, Annu. Rev. Immunol., $19(2001)$, $683$-$765$.
  * [47] P. Moreillon, P. A. Majcherczyk. ”Proinflammatory Activity of Cell-Wall Constituents From Gram-Positive Bacteria”, Scand. J. Infect. Dis., $35(9)$ $(2003)$, $632$-$641$.
  * [48] P. A. Morel, R. E. C. Lee, J. R. Faeder. ”Demystifying the Cytokine Network: Mathematical Models Point the Way”, Cytokine $98(2017)$, $115$-$123$.
  * [49] E. Ngondiep. ”A robust numerical two-level second-order explicit approach to predicting the spread of Covid-19 pandemic with undetected infectious cases”, J. Comput. Appl. Math., $403$ $(2022)$, $113852$.
  * [50] E. Ngondiep. ”Unconditional stability over long time intervals of a two-level coupled MacCormack/Crank-Nicolson method for evolutionary mixed Stokes-Darcy model”, J. Comput. Appl. Math., $409(2022)$, $114148$, Doi: 10.1016/j.cam.2022.114148.
  * [51] E. Ngondiep. ”Stability analysis of MacCormack rapid solver method for evolutionary Stokes-Darcy problem”, J. Comput. Appl. Math. $345(2019)$, $269$-$285$.
  * [52] E. Ngondiep. ”A two-level fourth-order approach for time-fractional convection-diffusion-reaction equation with variable coefficients”, Commun. Nonlinear Sci. Numer. Simul., $111(2022)$, $106444$, Doi: 10.1016/j.cnsns.2022.106444.
  * [53] E. Ngondiep. ”A novel three-level time-split MacCormack scheme for two-dimensional evolutionary linear convection-diffusion-reaction equation with source term”, Int. J. Comput. Math., $98(1)$ $(2021)$, $47$-$74$.
  * [54] E. Ngondiep. ”A novel three-level time-split approach for solving two-dimensional nonlinear unsteady convection-diffusion-reaction equation”, J. Math. Computer Sci., $26(3)$ $(2022)$, $222$-$248$.
  * [55] E. Ngondiep. ”An efficient three-level explicit time-split approach for solving two-dimensional heat conduction equation”, Appl. Math. Inf. Sci., $14(6)$, $(2020),$ $1075$-$1092$.
  * [56] E. Ngondiep. ”Long time stability and convergence rate of MacCormack rapid solver method for nonstationary Stokes-Darcy problem”, Comput. Math. Appl., $75$ $(2018)$, $3663$-$3684$.
  * [57] E. Ngondiep. ”An efficient three-level explicit time-split scheme for solving two-dimensional unsteady nonlinear coupled Burgers equations”, Int. J. Numer. Methods Fluids, $92(4)$ $(2020)$, $266$-$284$.
  * [58] E. Ngondiep. ”A robust three-level time-split MacCormack scheme for solving two-dimensional unsteady convection-diffusion equation”, J. Appl. Comput. Mech., $7(2)$ $(2021)$, $559$-$577$.
  * [59] E. Ngondiep. ”Long time unconditional stability of a two-level hybrid method for nonstationary incompressible Navier-Stokes equations”, J. Comput. Appl. Math., $345(2019)$, $501$-$514$.
  * [60] E. Ngondiep. ”A two-level factored Crank-Nicolson method for two-dimensional nonstationary advection-diffusion equation with time dependent dispersion coefficients and source/sink term”, Adv. Appl. Math. Mech., $13(5)$ $(2021)$, $1005$-$1026$.
  * [61] E. Ngondiep, N. Kerdid, M. A. M. Abaoud, I. A. I. Aldayel. ”A three-level time-split MacCormack method for two-dimensional nonlinear reaction-diffusion equations”, Int. J. Numer. Meth. Fluids, $92(12)$ $(2020)$, $1681$-$1706$.
  * [62] E. Ngondiep. ”A six-level time-split Leap-Frog/Crank-Nicolson approach for two-dimensional nonlinear time-dependent convection-diffusion-reaction equation”, Int. J. Comput. Meth., $(2023)$, Doi: 10.1142/S0219876222500645.
  * [63] H. F. Nijhout, J. A. Best, M. C. Reed. ”Using Mathematical Models to Understand Metabolism, Genes, and Disease”, BMC Biol., $13(1)$ $(2015)$, $79$.
  * [64] Y. Osawa, M. Nagaki, Y. Banno, D. A. Brenner, T. Asano, Y. Nozawa, et al. ”Tumor Necrosis Factor Alpha-Induced Interleukin-8 Production via NF-kappaB and Phosphatidylinositol 3-Kinase/Akt Pathways Inhibits Cell Apoptosis in Human Hepatocytes”, Infection Immun., $70(11)$ $(2002)$, $6294$-$6301$.
  * [65] A. Patil, D. Huard, C. Fonnesbeck. ”PyMC: Bayesian stochastic modelling in Python”, J. Stat. Softw., $(2010)$, $35:i04$.
  * [66] B. Petrini, B. Andersson, O. Strannegard, J. Wasserman, H. Blomgren, U. Glas. ”Monocyte release and plasma levels of inter-leukin-6 in patients irradiated for cancer”, In Vivo, $6(1992)$, $531$-$4$.
  * [67] B. D. M. Quintela, R. W. dos Santos, M. Lobosco. ”On the Coupling of Two Models of the Human Immune Response to an Antigen”, BioMed. Res. Int., $(2014)$, $410$-$457$.
  * [68] A. Ramoji, O. Ryabchykov, K. Galler, A. Tannert, R. Markwart, R. P. Requardt, I. Rubio, M. Bauer, T. Bocklitz, J. Popp, et al. ”Raman Spectroscopy Follows Time-Dependent Changes in T Lymphocytes Isolated from Spleen of Endotoxemic Mice”, ImmunoHorizons, $3(2019)$, $45$-$60$.
  * [69] A. Ramoji, U. Neugebauer, T. Bocklitz, M. Foerster, M. Kiehntopf, M. Bauer, J. Popp. ”Toward a Spectroscopic Hemogram: Raman Spectroscopic Differentiation of the Two Most Abundant Leukocytes from Peripheral Blood”, Anal. Chem., $84(2012)$, $5335$-$5342$.
  * [70] R. Smith. ”Uncertainty Quantification: Theory, Implementation, and Applications”, Philadelphia, PA: SIAM-Society for Industrial and Applied Mathematics, $(2013)$.
  * [71] K. Talaei, S. A. Garan, B. D. M. Quintela, M. S. Olufsen, J. Cho, J. R. Jahansooz, P. K. Bhullar, et al. ”A Mathematical Model of the Dynamics of Cytokine Expression and Human Immune Cell Activation in Response to the Pathogen Staphylococcus aureus”, Frontier Cell. Infect. Microb., $11(2021)$.
  * [72] D. M. Tawfik, J. M. Lankelma, L. Vachot, E. Cerrato, A. Pachot, W. J. Wiersinga, et al. ”Comparison of Host Immune Responses to LPS in Human Using an Immune Profiling Panel, In Vivo Endotoxemia Versus Ex Vivo Stimulation”, Sci. Rep., $10(1)$ $(2020)$, $9918$.
  * [73] M. Torres, J. Wang, P. J. Yannie, S. Ghosh, R. A. Segal, A. M. Reynolds. ”Identifying Important Parameters in the Inflammatory Process With a Mathematical Model of Immune Cell Influx and Macrophage Polarization”, PloS Comput. Biol., $15(7)$ $(2019)$, $e1007172$.
  * [74] N. K. Vitanov. ”Modified mothod of simplest equation: powful tool for obtaining exact and approximate traveling-wave solutions of nonlinear PDEs”, Commun. Nonlinear Sc. Numer. Simulat., $16(2011)$, $1176$-$85$.
  * [75] J. E. Wang, P. F. J$\phi$rgensen, M. Almlöf, C. Thiemermann, S. J. Foster, A. O. Aasen, et al. ”Peptidoglycan and Lipoteichoic Acid From Staphylococcus Aureus Induce Tumor Necrosis Factor Alpha, Interleukin 6 (IL-6), and IL-10 Production in Both T Cells and Monocytes in a Human Whole Blood Model”, Infect. Immun., $68(7)$ $(2000)$, $3965$-$3970$.
  * [76] D. A. Winer, S. Winer, L. Shen, P. P. Wadia, J. Yantha, G. Paltser, et al. ”B Cells Promote Insulin Resistance Through Modulation of T Cells and Production of Pathogenic IgG Antibodies”, Nat. Med., $17(5)$ $(2011)$, $610$-$617$.
  * [77] X. Xue, D. M. Falcon. ”The Role of Immune Cells and Cytokines in Intestinal Wound Healing”, Int. J. Mol. Sci., $20(23)$ $(2019)$, $6097$.
  * [78] L. Yao, F. D. Lowy, J. W. Berman. ”Interleukin-8 Gene Expression in Staphylococcus Aureus-Infected Endothelial Cells”, Infection Immun., $64(8)$ $(1996)$, $3407$-$3409$.
  * [79] H. Yssel, R. De Waal Malefyt, M. G. Roncarolo, J. S. Abrams, R. Lahesmaa, H. Spits, J. E. de Vries. ”IL-10 is produced by subsets of human CD4+ T cell clones and peripheral blood T cells”, J. Immunol., $149(1992)$, $2378$-$84$.
  * [80] J. M. Zhang, J. An. ”Cytokines, Inflammation and Pain”, Int. Anesthesiol Clin., $45(2)$ $(2007)$, $27$-$37$.

Approximate solutions at different spatial points $x$.  $x=(2\times
10^{-3},\text{\,}10^{-3},\text{\,}2\times 10^{-3})$ $x=(3\times
10^{-3},\text{\,}1.5\times 10^{-3},\text{\,}3\times 10^{-3})$ $x=(2\times
10^{-2},\text{\,}5\times 10^{-3},\text{\,}10^{-3})$ $x=(2\times
10^{-3},\text{\,}3\times 10^{-3},\text{\,}6\times 10^{-3})$

Figure 2:  Concentrations: S. aureus, resting macroph., activated macroph.,
TNF alpha, IL6, IL8 and IL10

Simulated results at point $x=(2\times
10^{-3},\text{\,}10^{-3},\text{\,}2\times 10^{-3})$

Figure 3:  Concentrations: S. aureus, resting macroph., activated macroph.,
TNF alpha, IL6, IL8 and IL10

Diffusion of S. aureus, macrophages and cytokines as provided in [71]
Parameter adjust.: TNF$\alpha$(A), IL6 (B), IL8 (C), IL10 (D) S. aureus
diffusion at different periods rest. macroph. diff.: oh(A), 3h(B), 12h(C),
24h(D) act. macroph. diff.: oh(A), 3h(B), 12h(C), 24h(D)

Figure 4:  TNF alpha, IL6, IL8, IL10, S. aureus, resting macroph. and
activated macroph.