Mathematical Model of Infection Propagation Mediated by Circulating Macrophage

1. Introduction

Viral infections remain a major global health concern due to their rapid transmissibility and potential to cause significant morbidity and mortality. Epidemics such as seasonal influenza and pandemics like COVID-19 have underscored the urgency of understanding the biological mechanisms underlying virus-host interactions, particularly during the early stages of infection [26,30]. Once inside host tissue, viruses infect susceptible cells, hijack their replication machinery, and produce progeny virions. These viral particles propagate through local diffusion and direct intercellular transmission, forming spatial infection waves [2,30].

The immune system responds with a layered defense strategy involving both innate and adaptive responses. The innate immune system acts rapidly, within hours, and provides a generalized defense before the adaptive immune system becomes fully engaged. Central to this early response are inflammatory cytokines and immune cells such as macrophages, which shape the tissue environment and limit viral replication [13,17].

Inflammatory cytokines and programmed cell death.

Upon detecting viral components, infected cells initiate signaling cascades leading to the production of cytokines such as tumor necrosis factor (TNF), interleukin-6 (IL-6), and interferons. These molecules orchestrate immune responses and trigger programmed cell death (PCD) in infected cells, thereby halting viral replication [18]. Several forms of PCD have been implicated in viral infections, including apoptosis, pyroptosis, necroptosis, and PANoptosis—a combined form integrating features of all three. Although essential, dysregulated or excessive cell death may result in tissue damage and systemic inflammation [10,17].

In our model, inflammatory cytokines are secreted by both infected cells and activated macrophages. These cytokines enhance immune-mediated destruction of infected cells, establishing a feedback loop that can contain infection, depending on parameter values.

Macrophage recruitment and activation.

Macrophages are key innate immune effectors involved in pathogen detection, phagocytosis, cytokine production, and immune regulation. During infection, blood monocytes are recruited to the tissue and differentiate into macrophages. These cells exhibit phenotypic plasticity, polarizing into M1 (pro-inflammatory, microbicidal) or M2 (anti-inflammatory, tissue-repairing) states based on environmental signals [13,22].

Here, we model circulating macrophages that become activated through interactions with virus particles and pro-inflammatory cytokines. Once activated, they amplify the cytokine pool and facilitate the clearance of infected cells. Their population dynamics are shaped by production, activation, and degradation processes.

Modeling infection propagation in tissue.

The spatial spread of infection arises from the interaction between viral replication, diffusion, and immune activity. Reaction-diffusion models provide a powerful framework to study such phenomena, enabling the analysis of wavefront propagation and spatiotemporal dynamics [2,7,23]. The wave speed reflects the efficiency of viral expansion, while the total viral load (integrated spatial virus concentration) quantifies infection burden [23].

We extend classical models by integrating immune regulation through cytokine signaling and macrophage activation. Our system consists of five key components: uninfected cells (E), infected cells (F), virus particles (V), inflammatory cytokines (S), and activated macrophages (N). Virus and infected cells are modeled as spatially distributed quantities, while S and N are treated as spatially homogeneous due their redistribution by blood circulation [7,19].

This study presents a novel spatial model of early viral infection incorporating immune regulatory feedback. We incorporate time delay representing intracellular viral replication and couple nonlinear dynamics of cytokine production and immune-mediated cell death. Next, we consider activation of macrophages by virus and cytokines, contributing to immune feedback. Analyzing the model, we determine the conditions of infection spreading and find the wave speed and viral load.

Numerical simulations complement the theoretical analysis and demonstrate how macrophage degradation rates and cytokine feedback strength can critically affect infection dynamics. Our results provide insights into the thresholds that separate controlled from uncontrolled viral propagation.

The remainder of the paper is structured as follows. Section 2 presents the mathematical formulation of the model. Section 3 analyzes the traveling wave solutions and conditions for their existence. Section 4 provides numerical simulations and investigates the impact of immune parameters. Section 5 discusses biological implications and outlines future research directions.

Figure 1. Schematic representation of the model. Healthy cells E are infected by the virus V and become infected cells F, which produce more virus and stimulate the immune response. The virus activates circulating macrophages M into activated macrophages N, which in turn release inflammatory cytokines S. These cytokines amplify macrophage activation and promote the elimination of infected cells.

Figure 1. Schematic representation of the model. Healthy cells E are infected by the virus V and become infected cells F, which produce more virus and stimulate the immune response. The virus activates circulating macrophages M into activated macrophages N, which in turn release inflammatory cytokines S. These cytokines amplify macrophage activation and promote the elimination of infected cells.

2. Reaction-Diffusion Model of Viral Infection

2.1. Model Formulation

The propagation of viral infection in an infected tissue occurs at the first stage of the innate immune response. The replication of the virus within infected cells promotes the production of inflammatory cytokines by both infected cells and circulating macrophages. Inflammatory cytokines, in turn, trigger the programmed cell death of infected cells. We consider the following system of equations:

$${\displaystyle \frac{\partial U}{\partial t}}=-a_{1}UV,$$

(2.1)

$${\displaystyle \frac{\partial F}{\partial t}}=a_{1}UV-k_{1}FS-k_{3}FN-{\sigma }_{1}F,$$

(2.2)

$${\displaystyle \frac{\partial V}{\partial t}}=D{\displaystyle \frac{{\partial }^{2}V}{\partial x^2}}+b{F}_{\tau }-k_{4}NV-{\sigma }_{3}V,$$

(2.3)

$${\displaystyle \frac{\mathrm{d}N}{\mathrm{d}t}}=h(J\left(V\right)+\beta S)(m_0-N)-{\sigma }_{2}N,$$

(2.4)

$${\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=a_{2}J\left(F\right)-k_{2}J\left(F\right)S+dN-{\sigma }_{4}S,$$

(2.5)

where $F_{\tau }(x,t)=F(x,t-\tau )$,

$${\displaystyle J\left(F\right)={\int }_{-\infty }^{+\infty }F(x,t)\mathrm{d}x,J\left(V\right)={\int }_{-\infty }^{+\infty }V(x,t)\mathrm{d}x.}$$

Here, $E(x,t)$ represents the concentration of uninfected target cells, such as epithelial cells, while $F(x,t)$ denotes the concentration of infected cells, and $N\left(t\right)$ corresponds to activated circulating macrophages. The virus concentration is represented by $V(x,t),$ and the inflammatory cytokines concentration is denoted by $S\left(t\right).$ Unlike $E,F,$ and $V,$ which depend on both time and space, the variables S and N depend only on time due to their distribution within the tissue via blood circulation and diffusion.

The right-hand side of equation (2.1) represents the rate of infection of uninfected cells by virus. A similar term enters the next equation with sign plus. The other terms in the right-hand side of this equation describe the rate of infected cell death due to inflammatory cytokines, their elimination by activated macrophages, and their death independent of the immune response. Equation (2.3) for the virus concentration in the infected tissue describes, respectively, its diffusion, production by the infected cells, its elimination by macrophages and its death independent of the immune response. The rate of macrophage activation by virus antigens and by inflammatory cytokines is given in the right-hand side of equation (2.4). The activation rate is limited by the total level of circulating macrophages $m_0$. The last term in this equation represents deactivation or death of macrophages. Note that the level of virus antigens in the body, which determines the activation of macrophages, is proportional to the total virus concentration $J\left(V\right)$ in the tissue. The last equation for inflammatory cytokines includes the rate of their production by infected cells proportional to their total quantity $J\left(F\right)$ in the tissue, their consumption in the interaction with infected cells, production by activated macrophages, and their degradation.

2.2. Infection Spreading as a Wave

We are looking for a solution of system (2.1)–(2.5) in the form of a travelling wave: $U(x,t)=u\left(\xi \right),F(x,t)=f\left(\xi \right),$ and $V(x,t)=v\left(\xi \right),$ where $\xi =x-ct$ and c is the wave speed. Since the two last equations of the system do not depend on the space variable, then ${\displaystyle \frac{\mathrm{d}N}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=0.$ We obtain:

$$c{u}^{\prime }=a_{1}uv,$$

(2.6)

$$-c{f}^{\prime }=a_{1}uv-k_{1}fS-k_{3}fN-{\sigma }_{1}f,$$

(2.7)

$$-c{v}^{\prime }=D{v}^{\prime \prime }+b{f}_{\tau }-{\sigma }_{3}v-k_{4}Nv,$$

(2.8)

$$N={\displaystyle \frac{h(J\left(v\right)+\beta S)m_0}{{\sigma }_2+h(J\left(v\right)+\beta S)}},$$

(2.9)

$$S={\displaystyle \frac{a_{2}J\left(f\right)+dN}{k_{2}J\left(f\right)+{\sigma }_4}},$$

(2.10)

where $f_{\tau }=f(\xi +c\tau ),$${\displaystyle J\left(f\right)={\int }_{-\infty }^{+\infty }f\left(\xi \right)\mathrm{d}\xi }$ and ${\displaystyle J\left(v\right)={\int }_{-\infty }^{+\infty }v\left(\xi \right)\mathrm{d}\xi .}$ We consider the following limits at infinity:

$$u(-\infty )=u_f,u(+\infty )=u_0,f(\pm \infty )=v(\pm \infty )=0,$$

(2.11)

where $u_0$ is the initial concentration of the uninfected cells, and $u_f$ is their unknown final concentration. Along with this final concentration, we will determine the integral the viral load $J\left(v\right)$ and the wave speed c.

3. System Without Degradation of Activated Macrophages $({\sigma }_2=0)$

We begin the study of system (2.6)-(2.10) with the case ${\sigma }_2=0$. This assumption simplifies the analysis and will allow us to determine $u_f$ from a single algebraic equation. Under this assumption, we have $N=m_0$. Therefore, we consider the following system:

$$c{u}^{\prime }=a_{1}uv,$$

(3.1)

$$-c{f}^{\prime }=a_{1}uv-k_{1}fS-k_{3}f{m}_0-{\sigma }_{1}f,$$

(3.2)

$$-c{v}^{\prime }=D{v}^{\prime \prime }+b{f}_{\tau }-({\sigma }_3+k_4{m}_0)v,$$

3.3)

$$S={\displaystyle \frac{a_{2}J\left(f\right)+d{m}_0}{k_{2}J\left(f\right)+{\sigma }_4}}.$$

(3.4)

Integrating (3.1), we obtain:

$$cln\left(\eta \right)=-a_{1}J\left(v\right).$$

Integration of the sum of equations (3.1) and (3.2) gives

$$c(u_0-u_f)=(k_{1}S+{\sigma }_1+k_3{m}_0)J\left(f\right).$$

(3.5)

From equation (3.3), we express $J\left(f\right)$ as follows:

$$J\left(f\right)={\displaystyle \frac{{\sigma }_3+k_4{m}_0}{b}}J\left(v\right).$$

Thus, we obtain the following system of equations with respect to the variables $\eta ={\displaystyle \frac{u_f}{u_0}}$ ($0<\eta <1$), $J\left(f\right)$ and $J\left(v\right)$:

$$\begin{array}{cc}\hfill cln\left(\eta \right)& =-a_{1}J\left(v\right),\hfill \end{array}$$

(3.6)

$$\begin{array}{cc}\hfill c{u}_0(1-\eta )& =\left(k_1{\displaystyle \frac{a_{2}J\left(f\right)+d{m}_0}{k_{2}J\left(f\right)+{\sigma }_4}}+{\sigma }_1+k_3{m}_0\right)J\left(f\right),\hfill \end{array}$$

(3.7)

$$\begin{array}{cc}\hfill J\left(f\right)& ={\displaystyle \frac{{\sigma }_3+k_4{m}_0}{b}}J\left(v\right).\hfill \end{array}$$

(3.8)

From (3.6) and (3.8), we obtain:

$$J\left(f\right)=-{\displaystyle \frac{c({\sigma }_3+k_4{m}_0)}{a_{1}b}}ln\left(\eta \right).$$

Substituting this expression into (3.7), we derive the equation with respect to $\eta$:

$$R(\eta -1)=ln\left(\eta \right)(1+g(ln\left(\eta \right)\left)\right),$$

(3.9)

where

$$R={\displaystyle \frac{a_{1}b{u}_0}{({\sigma }_1+k_3{m}_0)({\sigma }_3+k_4{m}_0)}},g\left(y\right)=\alpha {\displaystyle \frac{\beta -y}{\gamma -y}},$$

$$\alpha ={\displaystyle \frac{a_2{k}_1}{k_2({\sigma }_1+k_3{m}_0)}},\beta ={\displaystyle \frac{a_{1}bd{m}_0}{a_{2}c({\sigma }_3+k_4{m}_0)}},\gamma ={\displaystyle \frac{a_{1}b{\sigma }_4}{k_{2}c({\sigma }_3+k_4{m}_0)}}.$$

3.1. Existence of Solution

We will now study the existence of solutions of equation (3.9) for $\eta \in (0,1)$.

Proposition 3.1.  If

$$R>1+g\left(0\right)=1+{\displaystyle \frac{k_{1}d{m}_0}{{\sigma }_4({\sigma }_1+k_3{m}_0)}},$$

(3.10)

then equation (3.9) has a solution $\eta \in (0,1)$.

Proof. Set $f\left(y\right)=y(1+g(y\left)\right)$. Then $f(ln(\eta \left)\right)=0$ for $\eta =1$, and $f(ln(\eta \left)\right)\to -\infty$ as $\eta \to 0$. Therefore, if

$${\displaystyle \frac{d\left(R\right(\eta -1\left)\right)}{d\eta }}{{|}_{\eta =0}>{\displaystyle \frac{df(ln(\eta \left)\right)}{d\eta }}|}_{\eta =0},$$

then equation (3.9) has a solution $\eta \in (0,1)$. The last inequality is equivalent to (3.10). The proposition is proved.

□

In order to determine the number of solutions of equation (3.9), we find the second derivative $d^{2}f(ln\left(\eta \right))/d{\eta }^2$. We have

$${\displaystyle \frac{df(ln(\eta \left)\right)}{d\eta }}={\displaystyle \frac{1}{\eta }}{f}^{\prime }\left(y\right),{\displaystyle \frac{d^{2}f(ln\left(\eta \right))}{d{\eta }^2}}={\displaystyle \frac{1}{{\eta }^2}}\left(f^{\prime \prime }\left(y\right)-f^{\prime }\left(y\right)\right),$$

where $y=ln\left(\eta \right)$. Next,

$$f^{\prime \prime }\left(y\right)-f^{\prime }\left(y\right)={\displaystyle \frac{p_1{y}^3+p_2{y}^2+p_{3}y+p_4}{{(\gamma -y)}^3}}$$

where

$$p_1=\alpha +1,p_2=-\alpha \gamma -3\gamma ,p_3=\alpha \beta \gamma +2\alpha {\gamma }^2+3{\gamma }^2,p_4=\gamma (2\alpha \beta -2\alpha \gamma -\alpha \beta \gamma -{\gamma }^2).$$

Proposition 3.2.  If

$$2\alpha \beta \le 2\alpha \gamma +\alpha \beta \gamma +{\gamma }^2,$$

(3.11)

then

$${\displaystyle \frac{d^{2}f(ln\left(\eta \right))}{d{\eta }^2}}<0,\eta \in (0,1).$$

(3.12)

If

$$2\alpha \beta >2\alpha \gamma +\alpha \beta \gamma +{\gamma }^2,$$

(3.13)

then there is such ${\eta }_0\in (0,1)$ that

$${\displaystyle \frac{d^{2}f(ln\left(\eta \right))}{d{\eta }^2}}<0,\eta \in (0,{\eta }_0),{\displaystyle \frac{d^{2}f(ln\left(\eta \right))}{d{\eta }^2}}>0,\eta \in ({\eta }_0,1).$$

(3.14)

Proof. The sign of the second derivative of $f(ln(\eta \left)\right)$ for $\eta \in (0,1)$ coincides with the sign of the difference $f^{\prime \prime }\left(y\right)-f^{\prime }\left(y\right)$ for $y<0$, or with the sign of the polynomial

$$P\left(y\right)=p_1{y}^3+p_2{y}^2+p_{3}y+p_4,y<0.$$

Note that $p_1>0,p_2<0,p_3>0$. Therefore,

$$P_0\left(y\right)\equiv p_1{y}^3+p_2{y}^2+p_{3}y<0,y<0.$$

If condition (3.11) is satisfied, then $p_4\le 0$ and $P\left(y\right)<0$ for all $y<0$. If condition (3.11) is satisfied, then $p_4>0$. Since $P_0\left(y\right)$ is a monotonically increasing function for $y<0$, then there is some $y_0<0$ such that $P\left(y\right)>0$ for $y_0<y<0$ and $P\left(y\right)<0$ for $y<y_0$. This property of the function $P\left(y\right)$ provides (3.14). The proposition is proved.

□

Theorem 3.3.  If condition (3.11) is satisfied, then equation (3.9) has a solution $\eta \in (0,1)$ if and only if inequality (3.10) holds. Moreover, this solution is unique. If condition (3.13) is satisfied, then there exists such $R^{*}\in (0,1+g\left(0\right))$ that equation (3.9) has a solution $\eta \in (0,1)$ in and only if $R\ge R^{*}$. This solution is unique for $R\ge 1+g\left(0\right)$, and there are two solutions for $R^{*}<R<1+g\left(0\right)$.

Proof. Let (3.11) be satisfied. According to Proposition 3.2, the function $h\left(\eta \right)=f(ln(\eta \left)\right)$ is concave for $\eta \in (0,1)$. Since equation (3.9), that is,

$$R(\eta -1)=h\left(\eta \right)$$

has solution $\eta =0$, then it has a unique solution $\eta \in (0,1)$ if $R>h^{\prime }\left(1\right)$, and no such solution otherwise.

Let now (3.13) be satisfied. Then existence of solution for $R>1+g\left(0\right)$ follows from Proposition 3.1. Suppose that the solution is not unique and, for certainty, there are two solutions ${\eta }_1,{\eta }_2\in (0,1)$, ${\eta }_1<{\eta }_2<0$. The derivatives of the left-hand side and of the right-hand side of equation (3.9) are equal to each other between solutions ${\eta }_1$ and ${\eta }_2$, and between solutions ${\eta }_2$ and $\eta =0$. Therefore, $h^{\prime }\left(\eta \right)=R$ at least for two values of $\eta \in (0,1)$. According to Proposition 3.2, the $h^{\prime }\left(\eta \right)$ is a positive function in the interval $(0,1)$ with a single minimum. Since $h^{\prime }\left(1\right)<R$, then equation $h^{\prime }\left(\eta \right)=R$ can have at most one solution. This contradiction proves uniqueness of solution of equation (3.9) in $(0,1)$.

Similar arguments show that equation (3.9) does not have solution in $(0,1)$ for R sufficiently small, in particular, if $R<{min}_{\eta \in (0,1)}{h}^{\prime }\left(\eta \right)$. Hence, there is such $R^{*}<1+g\left(0\right)$ that equation (3.9) does not have solutions for $R<R^{*}$.

Finally, there are exactly two solutions for $R^{*}<R<h^{\prime }\left(1\right)=1+g\left(0\right)$ (Figure 2). This follows from a similar analysis of the derivatives $h^{\prime }\left(\eta \right)$ and $h^{\prime \prime }\left(\eta \right)$. The details are left to the reader. The theorem is proved.

□

Approximate solution.

Analysis of equation (3.9) shows that its solution converges to 0 as R increases. Therefore, for R sufficiently large, $\eta \ll 1$, and $\eta$ can be neglected in the left-hand side of this equation. On the other hand, $g(ln(\eta \left)\right)\to g(-\infty )=\alpha$. Hence, we obtain approximate equality $R=-\left(1+\alpha \right)ln\left(\eta \right)$. Thus, we can find $\eta ,J\left(f\right),J\left(v\right):$

$$\eta =e^{-R/(1+\alpha )},J\left(f\right)={\displaystyle \frac{cR({\sigma }_3+k_4{m}_0)}{a_{1}b\alpha }},J\left(v\right)={\displaystyle \frac{cR}{a_1\alpha }}.$$

Note that $J\left(f\right)$ and $J\left(v\right)$ depend on the wave speed c. We will determine it in the next section.

3.2. Wave Speed

In order to determine the wave speed, we use the linearization method. We replace u by $u_0$ at $\infty .$ Therefore, we obtain the following linearized system of equations for f and $v:$

$$c{f}^{\prime }=a_1{u}_{0}v-\left({\displaystyle \frac{k_1{a}_{2}J\left(f\right)+k_{1}d{m}_0}{k_{2}J\left(f\right)+{\sigma }_4}}+k_3{m}_0+{\sigma }_1\right)f,$$

(3.15)

$$D{v}^{\prime \prime }+c{v}^{\prime }+b{f}_{\tau }-({\sigma }_3+k_4{m}_0)v=0.$$

(3.16)

The integral $J\left(f\right)$ is considered here as a given constant. Let us look for the solution of this system in the form $f\left(\xi \right)=p_1{e}^{-\lambda \xi },v\left(\xi \right)=p_2{e}^{-\lambda \xi }.$ Substituting these functions into equations (3.15), (3.16), we obtain:

$$-c\lambda p_1+a_1{u}_0{p}_2-\left({\displaystyle \frac{k_1{a}_{2}J\left(f\right)+k_{1}d{m}_0}{k_{2}J\left(f\right)+{\sigma }_4}}+k_3{m}_0+{\sigma }_1\right)p_1=0,$$

(3.17)

$$D{\lambda }^2{p}_2-c\lambda p_2+b{e}^{c\lambda \tau }{p}_1-({\sigma }_3+k_4{m}_0)p_2=0.$$

(3.18)

In order to find the minimal wave speed, we should determine the minimal value of c for which this system of equations has a positive solution $\lambda$ since the solution is decreasing at infinity. Introducing an independent parameter $\mu =\lambda c$ and excluding $p_1,p_2$ from the previous two equations, we obtain the following equation:

$$(D{(\mu /c)}^2-\mu -({\sigma }_3+k_4{m}_0))\left(\mu +H+k_3{m}_0+{\sigma }_1\right)+a_{1}b{u}_0{e}^{-\mu \tau }=0,$$

where

$$H={\displaystyle \frac{k_1{a}_{2}J\left(f\right)+k_{1}d{m}_0}{k_{2}J\left(f\right)+{\sigma }_4}}.$$

Hence,

$$c=\underset{\mu >{\mu }_0}{min}\mathsf{\Phi }\left(\mu \right),$$

(3.19)

where

$$\mathsf{\Phi }\left(\mu \right)={\displaystyle \frac{\sqrt{D}\mu \sqrt{\mu +{\sigma }_1+k_3{m}_0+H}}{\sqrt{(\mu +{\sigma }_3+k_4{m}_0)(\mu +{\sigma }_1+k_3{m}_0+H)-a_{1}b{u}_0{e}^{-\mu \tau }}}}=$$

$${\displaystyle \frac{\sqrt{D}\mu \sqrt{\mu /({\sigma }_1+k_3{m}_0)+P}}{\sqrt{{\sigma }_3+k_4{m}_0}\sqrt{(\mu /({\sigma }_3+k_4{m}_0)+1)(\mu /({\sigma }_1+k_3{m}_0)+1+P)-R{e}^{-\mu \tau }}}}.$$

Here

$$P={\displaystyle \frac{H}{{\sigma }_1+k_3{m}_0}}={\displaystyle \frac{k_1{a}_{2}J\left(f\right)+k_{1}d{m}_0}{({\sigma }_1+k_3{m}_0)(k_{2}J\left(f\right)+{\sigma }_4)}},$$

and ${\mu }_0>0$ is the value for which the denominator vanish. Such value exists if

$$R>1+P=1+{\displaystyle \frac{k_1{a}_{2}J\left(f\right)+k_{1}d{m}_0}{({\sigma }_1+k_3{m}_0)(k_{2}J\left(f\right)+{\sigma }_4)}}.$$

(3.20)

If condition (3.20) is not satisfied, then the function $\mathsf{\Phi }\left(\mu \right)$ is defined for all $\mu \ge 0$, and $\mathsf{\Phi }\left(0\right)=0$. The minimal wave speed in this case equals zero.

Note that this condition is not explicit. Indeed, the integral $J\left(f\right)$ depend on the unknown wave speed. Therefore, P and the function $F\left(\mu \right)$ also depend on it, and (3.19) is an equation with respect to c. However, P can be estimated by the respective values for $J\left(f\right)=0$ and $J\left(f\right)=\infty$:

$$P_0={\displaystyle \frac{k_{1}d{m}_0}{{\sigma }_4({\sigma }_1+k_3{m}_0)}},P_{\infty }={\displaystyle \frac{k_1{a}_2}{k_2({\sigma }_1+k_3{m}_0)}}.$$

If

$$k_{2}d{m}_0<a_2{\sigma }_4,$$

(3.21)

then $P_0<P_{\infty }$. If inequality (3.21) is opposite, then $P_0$ is larger than $P\infty$. We will use these values below to determine different regimes of infection progression.

3.3. Regimes of Infection Progression

We derived above two critical conditions that determine infection progression. Condition (3.10) provides the existence of solution of equation (3.9) in the interval $\eta \in (0,1)$. Condition (3.20) implies that the wave speed is positive. Both of them are necessary conditions of the existence of infection wave, but their sufficiency is not proved.

Note that $g\left(0\right)=P_0$ such that condition (3.10) coincides with condition (3.20) if $J\left(f\right)=0$. Hence, if

$$R<1+min(P_0,P_{\infty }),$$

(3.22)

then neither of conditions (3.10) and (3.20) are satisfied, and infection wave does not exist.

However, infection can spread in a transient regime with gradual extinction. In order to illustrate this regime, let us consider system (2.1)-(2.5) under the assumption that $h=0,a_2=0$. This means that the concentrations of inflammatory cytokines and macrophages remain zero (as for the initial condition). Hence, this model is the same as in the case of infection spreading in cell culture [4]. It propagates if $R_0=a_{1}b{u}_0/\left({\sigma }_1{\sigma }_3\right)<1$. Note that $R_0$ is obtained from the virus replication number R if $m_0=0$.

We now return to the full system with $h>0$ assuming that $R<1$. In the beginning of the simulations, when the concentration of macrophages is low, infection propagates as for the case without macrophages (Figure 3). However, since their concentration grows with time, viral load decreases and converges to zero. Thus, in this example, infection is progressively eliminated by the innate immune response.

On the other hand, if

$$R>1+\max (P_0,P_{\infty }),$$

(3.23)

then both conditions (3.10) and (3.20) are satisfied, and we observe propagation of infection wave. Note that the analytical results obtained above do not prove the existence and stability of such wave but numerical simulations confirm infection propagation as a reaction-diffusion wave (Figure 4). The analytical and numerical values of the viral load and wave speed are in good agreement. The analytic value of viral load is $J\left(v\right)=65454.545$ and its numerical value is $J\left(v\right)=65909.090.$ The analytic value of the wave speed ic $c=0.144$ and the numerical value is $c=0.145.$

Figure 5 shows the dependence of the wave speed and the viral load on time delay in virus replication. Both of them decay with increasing $\tau .$

In order to continue the analysis of the regimes of infection progression, suppose that condition (3.21) is satisfied. Then $P_0<P_{\infty }$. We will consider now the case where condition (3.10) is satisfied and condition (3.20) is not satisfied. This means that

$$R>1+g\left(0\right)=1+P_0=1+{\displaystyle \frac{k_{1}d{m}_0}{{\sigma }_4({\sigma }_1+k_3{m}_0)}},$$

and

$$R<1+P=1+{\displaystyle \frac{k_1{a}_{2}J\left(f\right)+k_{1}d{m}_0}{({\sigma }_1+k_3{m}_0)(k_{2}J\left(f\right)+{\sigma }_4)}}.$$

For any $J\left(f\right)>0$, $P>P_0$, so that these two conditions are compatible. The first condition provides infection progression, while the second condition means that the wave speed equals zero, that is, infection wave extincts. We have previously observed such behavior in [9] for the model with resident macrophages.

4. System with Degradation of Activated Macrophages $({\sigma }_2\ne 0)$

We consider the system (2.6)-(2.10). As in Section 3, we derive the equation for $\eta$:

$$R(\eta -1)={\displaystyle \frac{{\sigma }_3+k_{4}N}{{\sigma }_3+k_4{m}_0}}·{\displaystyle \frac{k_{1}S+{\sigma }_1+k_{3}N}{{\sigma }_1+k_3{m}_0}}ln\left(\eta \right),$$

(4.1)

where

$$J\left(v\right)=-{\displaystyle \frac{cln\left(\eta \right)}{a_1}},J\left(f\right)=-{\displaystyle \frac{c({\sigma }_3+k_{4}N)}{a_{1}b}}ln\left(\eta \right),S={\displaystyle \frac{a_{2}J\left(f\right)+dN}{k_{2}J\left(f\right)+{\sigma }_4}},N={\displaystyle \frac{h(J\left(v\right)+\beta S)m_0}{{\sigma }_2+h(J\left(v\right)+\beta S)}}.$$

Note that for ${\sigma }_2=0$, we have $N=m_0$, and we reduce this equation to the previous one.

Set

$$\mathsf{\Phi }\left(\eta \right)={\displaystyle \frac{{\sigma }_3+k_{4}N}{{\sigma }_3+k_4{m}_0}}·{\displaystyle \frac{k_{1}S+{\sigma }_1+k_{3}N}{{\sigma }_1+k_3{m}_0}}.$$

Let us behavior of this function as $\eta \to 0$ assuming that $c>0$. Note that $J\left(v\right)\left(\eta \right)\to \infty$ and $J\left(f\right)\left(\eta \right)\to \infty$, since $N\left(\eta \right)\ge 0$. Therefore, $J\left(v\right)\left(\eta \right)+S\left(\eta \right)\to infty$, since $S\left(\eta \right)\ge 0$. Dividing numerator and denominator of N by this expression, we conclude that $N\left(\eta \right)\to m_0$. Hence, $S\left(\eta \right)\to a_2/k_2$. We conclude that

$$\mathsf{\Phi }\left(0\right)=1+{\displaystyle \frac{a_2{k}_1}{k_2({\sigma }_1+k_3{m}_0)}}.$$

Thus, the right-hand side of equation (4.1) converges to $-\infty$ as $\eta \to 0$. Therefore, this equation has a solution $\eta \in (0,1)$ if

$${\displaystyle \frac{d}{d\eta }}{(\mathsf{\Phi }\left(\eta \right)ln\left(\eta \right))}_{\eta =1}<R.$$

(4.2)

By virtue of approximation $ln\left(\eta \right)\sim \eta -1$ in the vicinity of $\eta =1$, we conclude that ${\mathsf{\Phi }}^{\prime }\left(1\right)$ is bounded. Hence, condition (3.5) is equivalent to the condition

$$R>\mathsf{\Phi }\left(1\right).$$

(4.3)

In order to determine $\mathsf{\Phi }\left(1\right)$, we note that $J\left(v\right)\left(1\right)=0$ and $J\left(f\right)\left(1\right)=0$, since $N\le m_0$. Therefore,

$$S={\displaystyle \frac{dN}{{\sigma }_4}},N={\displaystyle \frac{h\beta m_{0}S}{{\sigma }_2+h\beta S}}.$$

We reduce this system to the single equation with respect to S:

$$S={\displaystyle \frac{dh\beta m_{0}S}{{\sigma }_4({\sigma }_2+h\beta S)}}.$$

(4.4)

It has a positive solution $S^{*}$ if

$${\displaystyle \frac{dh\beta m_0}{{\sigma }_2{\sigma }_4}}>1.$$

(4.5)

This solution $S^{*}$ can be found as a solution of the corresponding quadratic equation.

We proved the following proposition.

Proposition 4.1.  Suppose that $c>0$ and condition (4.5) is satisfied. If

$$R>{\displaystyle \frac{{\sigma }_3+k_4{N}^{*}}{{\sigma }_3+k_4{m}_0}}·{\displaystyle \frac{k_1{S}^{*}+{\sigma }_1+k_3{N}^{*}}{{\sigma }_1+k_3{m}_0}},$$

(4.6)

then equation (4.1) has a solution $\eta \in (0,1)$. Here $S^{*}$ is a positive solution of equation (4.4), and $N^{*}={\sigma }_4{S}^{*}/d$. In the limit of small ${\sigma }_2$ condition (4.6) coincides with condition (3.10).

Set

$$R_c\left({\sigma }_2\right)={\displaystyle \frac{{\sigma }_3+k_4{N}^{*}}{{\sigma }_3+k_4{m}_0}}·{\displaystyle \frac{k_1{S}^{*}+{\sigma }_1+k_3{N}^{*}}{{\sigma }_1+k_3{m}_0}}.$$

Since $S^{*}$ is a decreasing function of ${\sigma }_2$ as solution of equation (4.4), then this is also true for $N^{*}$ and for $R_c\left({\sigma }_2\right)$. Hence, the critical value of R decreases with increase of the macrophage death rate.

If equation (4.1) has a solution $\eta \in (0,1)$, then we can determine viral load and wave speed (cf. Section 3). However, the calculations becomes excessively complex and do not lead to explicit analytical formula. Hence, we restrict their investigation to numerical simulations (Figure 6). Both viral load and wave speed increase with the increase of ${\sigma }_2$ and decrease with increase of $k_3$. Hence, increasing the macrophage death rate promotes infection propagation while elimination of infected cells by macrophages down-regulates it.

5. Discussion

In this study, we developed and analyzed a spatial mathematical model that captures the early dynamics of viral infection in tissue, explicitly incorporating key elements of the innate immune response — namely, inflammatory cytokines and circulating macrophages.

The model consists of a hybrid system coupling three reaction-diffusion partial differential equations (for healthy cells, infected cells, and free virus particles) with two spatially homogeneous ordinary differential equations (for cytokines and activated macrophages). This structure allows us to explore both local and systemic effects, particularly travelling infection fronts.

The model integrates the most essential features of viral infections and immune response. In particuler, infection of healthy cells by viral particles and virus production by infected cells, clearance of infected cells via cytokines and macrophages, Cytokine production by infected cells and activated macrophages, activation of macrophages by both cytokines and virus, and optional degradation of activated macrophages (modeled via ${\sigma }_2$).

Together, these interactions form a nonlinear feedback system, where immune activation is dynamically modulated by infection burden. This feedback is biologically supported by experimental studies on cytokine signaling and macrophage polarization [6,13,14,16,28].

An important property of the model is its ability to represent the feedback loop formed by cytokines and activated macrophages. Infected cells trigger cytokine release, which promotes macrophage activation. In turn, these macrophages enhance cytokine production and contribute to the elimination of infected cells [16,28]. This forms a positive feedback loop that may facilitate rapid viral clearance.

However, if the degradation of activated macrophages is included ( ${\sigma }_2>0$), the feedback becomes self-limiting. This addition reflects biologically realistic processes such as macrophage exhaustion, apoptosis, or immune regulation [14,17]. The loss of effector cells dampens cytokine production, weakening immune control and allowing broader viral spread.

Characterization of infection progression.

In the idealized scenario ${\sigma }_2=0$, activated macrophages persist indefinitely, and the system simplifies. Considering travelling wave of infection spreading, we derived an explicit nonlinear equation for the fraction of surviving uninfected cells $\eta \in (0,1)$. This equation links $\eta$ to the total viral load $J\left(v\right)$ and wave speed c. It provides a necessary condition for the infection spreading in the tissue in the form of travelling wave.

Solvability of this equation provides an important characterization of the progression of viral infection. It is determined by the virus replication number R expressed through the parameters of the model. In the case without macrophages ($m_0=0$) or if they do not eliminate infected cells and free virions ($k_3=k_4=0$), the virus replication number introduced in this work coincides with the virus replication number $R_0$ considered in [4] for infection spreading in cell cultures without immune response. As it was shown in this previous work, viral infection spreads in cell culture if $R_0>1$ and extincts otherwise. Solvability condition of this equation in the case of immune response provided by circulating macrophages is given by the inequality $R>1+P_0$, where $P_0$ is a positive parameter determined by the coefficients of the model. Hence, the region in the parameter space corresponding to infection progression shrinks due to the immune response while the region of infection elimination increases.

It should be noted, however, that condition $R>1+P_0$ is a necessary condition for the solvability of this equation but not sufficient. It appears that in some parameter range, it can be solvable in a narrow interval $R^{*}<R<1+P_0$. Moreover, there are two solutions of this equation satisfying $\eta \in (0,1)$. We can expect that there are two waves exist in this case, one of them stable and another one unstable, and that initial conditions should be sufficiently large to provide convergence to the stable wave versus extinction of the infection and convergence to zero. This conjecture is not yet confirmed by direct numerical simulations since the parameter range of possible bistability is narrow, and the choice of initial conditions should be very precise.

Besides solvability of equation for $\eta$, we determine the wave speed by the linearization method. The minimal wave speed $c_0$ is positive if $R>1+P$, where P is the strength of the immune response. This parameter is expressed through the coefficients of the model and through the viral load. The latter satisfies a transcendental system of equation (together with the wave speed), and it does not admit a simple analytical expression, but it can be approximated by the method of successive approximation. Note that, depending on parameters, we can have $P>P_0$ or $P<P_0$.

These three parameters, $R,P_0$ and P determine regimes of infection progression. If $R>1+P_0$ and $R>1+P$, then infection spreads as a travelling wave, while if both inequalities are opposite, it extincts. Though analytical results give only necessary conditions of infection spreading, these conclusions are confirmed by direct numerical simulations.

An interesting question is whether it is possible to have only one of these conditions satisfied. In this case, we can expect that infection starts its progression but gradually extincts, as it was previously observed in the model with resident macrophages [22]. In numerical simulations of the actual model, we could not yet find this regime of infection spreading.

One more regime corresponds to transient infection spreading and its gradual extinction. Since initially there are no activated macrophages ($N\left(0\right)=0$), then infection propagates if $R_0>1$. However, since the concentration of macrophages grows with time, then the viral load decreases. In the case $R<1$, which is possible since $R>R_0$, infection is gradually eliminated by the innate immune response (Figure 3).

If ${\sigma }_2>0$, the model becomes more difficult to analyze. We determine the virus replication number and sufficient conditions of solvability of equation for $\eta$. We conclude that higher values of ${\sigma }_2$ facilitate infection progression. This conclusion is confirmed by direct numerical simulations which show that the infection front propagates more rapidly and reaches higher viral burdens.

Biological impliciations and model limitations.

The results of the work corroborate biological data about the role of macrophages in the development of viral infections. During respiratory viral infections, circulating monocytes are rapidly recruited into the lungs, where they differentiate into macrophages that participate in viral clearance. These macrophages lower viral load through multiple mechanisms: direct phagocytosis of viral particles and infected epithelial cells, production of antiviral mediators such as nitric oxide and type I interferons, and enhancement of adaptive immunity by presenting viral antigens to T lymphocytes. Experimental studies in influenza virus infection demonstrate that depletion of circulating or alveolar macrophages results in higher viral titers and worsened pathology, while their activation promotes faster viral clearance and improved survival [15,29].

These findings emphasize the dual role of macrophage regulation: while necessary to avoid overactivation and tissue damage, excessive degradation may impair viral containment, especially in early infection [17,18]. Our results underscore a fundamental trade-off in immune dynamics: persistent activation of macrophages is beneficial for rapid viral control, but requires tight regulation to avoid immunopathology. Conversely, premature degradation of effector cells facilitates viral spread. This suggests that therapeutic strategies enhancing macrophage stability—or delaying their exhaustion—could improve immune effectiveness in early infection phases. Excessive immune suppression might unintentionally promote viral persistence.

Immune response to viral infection is a very complex biological processes, and the models of immune response encounter necessarily essential simplifications. In this work we consider only circulating macrophages, as one of the most important actor of the innate immune response. Resident macrophages were considered on [22], and adaptive immune response was considered in [27]. Bringing together these different parts of the immune response will require further investigations.

Conclusions.

By integrating spatial dynamics and immune feedback, our model highlights the central role of macrophage-mediated responses in shaping infection spread. The analytical and numerical results illustrate how spatially distributed infections can be contained or amplified depending on macrophage turnover and cytokine signaling strength. Future extensions may include adaptive immunity, tissue heterogeneity, and stochastic effects to further refine the understanding of immune-virus interplay.