1 INTRODUCTION

Besides public health, epidemics of viral infections can have important social, economic, and political consequences [5–7]. Ongoing COVID-19 pandemic urged more profound and comprehensive investigations of viral infections. Human organism reacts to viral infections with the immune response started with the innate immune response and continued with the adaptive immune response. The innate immune response develops immediately when a pathogen enters the body. It provides the body’s primary reaction to the pathogen and triggers the adaptive response mechanism. The adaptive immune response consists in the production of antibodies specific for a given antigen and of cytotoxic lymphocytes (CTL) eliminating infected cells. Both innate and adaptive immune responses have cellular and humoral components. The humoral component of the innate immune response includes interferon, which is a common name for a number of proteins secreted by the cells of the body and down-regulating virus replication. In the adaptive immune response, the humoral component includes antibodies, which are molecules that serve to neutralize the pathogen, and the cellular component includes B-cells (plasma cells) producing antibodies, CTL which kill the infected cells, T-helpers and some other cells.

Mathematical modeling is widely used for the investigation of viral infections [8, 11, 16–18]. It allows the identification of parameters and conditions of the viral infection spreading, which can be difficult to identify in the experiments, and evaluation of the influence of the immune response and of medical treatment on the infection progression. Conventional models of viral infections are based on ordinary differential equations (ODE) for the concentrations of virus, uninfected, and infected cell concentrations in time [8]. There are different developments of such models targeting various features of virus replication and cell infection for different viruses, and accounting for the antiviral therapy (see, e.g., [8] and references therein). In such models, it is implied that virus and susceptible cells are homogeneously distributed in space, but in solid tissues this assumption may not be satisfied and it can lead to systemic errors in the estimates of parameters underlying acute infection dynamic [13].

Reaction–diffusion models taking into account spatial distribution of virus concentration allow more accurate prediction of infection dynamics. In [14], the existence of solutions of an initial-boundary value problem for a reaction-diffusion model of viral infection is proved. The reaction–diffusion model of the virus plaque growth in the case of reversible host infection is considered in [15], and its wave propagation speed is determined by the linearization method. Numerical simulations of the reaction–diffusion system describing virus propagation with time delay are carried out in [16].

Innate and adaptive immune responses are often modeled independently of each other. As such, in [14, 17], the models take into account the dynamics of interferon. In [14], the existence of a solution is proved, but the effect of interferon on the development of viral infection is not investigated. In [17], the model includes interferon and it is used to investigate the early viral infection dynamics, i.e., the dynamics which is not observed in physiological experiments. This model takes into account spatial distribution of virus and interferon not with diffusion but in the frame of metapopulation approach: the model consists of \(n\) patches and they are coupled to each other by populations of cells migrating between them. In individual patches, the local dynamics is governed by an ODE system. Adaptive immune response (CTL) is taken into account in [13] with an ODE model, while spatial dynamics is considered through the migration between adjacent patches. The spatial dynamics of viral load and cells of immune response (T-helpers and CTL) is investigated in [18] with two-dimensional spatial patterns of viruses and T-cells. These models are mainly based on ODEs, they take into account only a part of the immune response (either innate or adaptive) and do not take into account time delay in viral replication. In addition, these models only consider the local influence of the immune response in the infected tissue.

In our previous works, it has been shown that infection propagates in cell culture as a reaction-diffusion wave [4]. The main characteristics of this wave are the spreading speed and the total viral load, i.e., the total quantity of virus in the culture at every moment of time. The spreading speed determines the part of the infected tissue and, as a consequence, the severity of symptoms of virus-induced disease. The total viral load in the upper respiratory tract for respiratory infections is proportional to virus infectivity, i.e., the rate of infection transmission between individuals. These characteristics are different, and smaller spreading speed can be associated both with lower or higher total viral load depending on virus type and cell culture [3, 9].

In this work, we further develop the approach introduced earlier [1, 4] and investigate how the innate and adaptive immune responses change the spreading speed and the total viral load. In the interpretation of the results, we refer to the respiratory viruses, primarily, to SARS-Cov-2 variants, but the conclusions remain valid for other infections for which this model is applicable.

Viral infection progression in   vivo is first confronted by the innate immune response including interferon produced by infected cells and down-regulating virus replication. In Section 2, we consider the influence of locally produced interferon and its diffusion in the surrounding tissue on infection spreading. We show that in this case the viral propagation speed is not influenced by the interferon concentration, but the total viral load decreases in comparison with case without interferon. Next, we consider globally circulating interferon in the organism (Section 3) and show that it decreases both infection spreading speed and viral load. In this case, depending on the values of parameters one or two stable waves can exist. We determine the conditions for these regimes, and study the influence of interferon concentration on the viral propagation speed and total viral load. We note that circulating interferon decreases the viral propagation speed in contrast to its local production. Finally, we include in the model the adaptive immune response mediated by antibodies and CTL (Section 4). In this case, both the severity and infectivity of the disease inversely depend on the strength of the immune response. We also observe that the concentration of globally circulating interferon depends on the humoral part of the adaptive immune response (antibodies) and on its cellular part (CTL).

The analytical results presented above are completed by numerical simulations. Numerical algorithm is based on the explicit first-order Euler scheme with \(P_{1}\) finite element spatial approximation with the mesh step \(L/10^{4}\) where \(L\) is the length of the modeled interval. We use free software FreeFEM [10] that offers a large variety of triangular finite elements (linear and quadratic Lagrangian elements, discontinuous \(P_{1}\), Raviart–Thomas elements, etc.) to solve partial differential equations. The values of parameters for the model of SARS-CoV-2 infection are taken from the previous works [1, 2, 4].

2 INTERFERON DIFFUSION IN CELL CULTURE

During viral infection, infected cells produce interferon which prevents virus replication. Some viruses, including SARS-CoV-2, can down-regulate interferon production. In order to describe infection progression in cell culture, we consider the system of time–delay differential equations

$$\frac{\partial U}{\partial t}=-aUV,$$
(1)
$$\frac{\partial I}{\partial t}=aUV-\beta I,$$
(2)
$$\frac{\partial V}{\partial t}=D_{1}\frac{\partial^{2}V}{\partial x^{2}}+\frac{b_{1}}{1+k_{1}C}\;I_{\tau_{1}}-\sigma_{1}V,$$
(3)
$$\frac{\partial C}{\partial t}=D_{2}\frac{\partial^{2}C}{\partial x^{2}}+\frac{b_{2}}{1+k_{2}V}\;I_{\tau_{2}}-\sigma_{2}C$$
(4)

for the concentrations of uninfected cells \(U\), infected cells \(I\), virus \(V\) and interferon \(C\). Equations for \(V\) and \(C\) contain production rates which depend on the concentration of infected cells with time delay, \(I_{\tau_{i}}(x,t)=I(x,t-\tau_{i}),\;i=1,2\). This system describes mutual inhibition of virus and interferon: the virus replication rate in equation (3) is inversely proportional to the interferon concentration \(C\), and virus can slow down the interferon production. The cell concentrations \(U\) and \(I\) are measured in cell/mL, the virus concentration \(V\) is measured in copy/mL, and the interferon concentration \(C\) is measured in unit/mL. Time \(t\) and space \(x\) are measured in hours and cm, respectively, and parameters have the corresponding units, which are specified below.

2.1 Interferon Does Not Change the Wave Speed

Let us recall that traveling wave solution of system (1)–(4) is a solution of the form \(\bar{u}=\bar{u}(x-ct)\) considered on the whole real axis. Here, \(\bar{u}\) is the vector of unknowns and c is the wave speed. In order to determine the wave speed, we introduce the moving reference frame \(\xi=x-ct\), where \(c\) is the wave speed, and linearize the system about its value at infinity, \(U=U_{0}\), \(I=V=C=0\), where \(U_{0}\) is the initial concentrations of uninfected cells. Then, we obtain the following system of equations with respect to the functions \(v(\xi)\), \(w(\xi)\), \(z(\xi)\):

$$cw^{\prime}+au_{0}v-\beta w=0,$$
(5)
$$D_{1}v^{\prime\prime}+cv^{\prime}+g_{1}(0)w(\xi+c\tau_{1})-\sigma_{1}v=0,$$
(6)
$$D_{2}z^{\prime\prime}+cz^{\prime}+g_{2}(0)w(\xi+c\tau_{2})-\sigma_{2}z=0,$$
(7)

where prime denotes derivation with respect to \(\xi\), \(g_{1}(z)=\frac{b_{1}}{1+k_{1}z}\), and \(g_{2}(v)=\frac{b_{2}}{1+k_{2}v}\). Consider its solution in the form \(v=pe^{-\lambda\xi}\), \(w=qe^{-\lambda\xi}\) and substitute it into the equations (5), (6). The linear homogeneous system of equations has a nonzero solution if and only if the corresponding determinant equals zero:

$$\left|\begin{matrix}-c\lambda-\beta&au_{0}\\ g_{1}(0)e^{-\lambda c\tau_{1}}&-c\lambda+D_{1}\lambda^{2}-\sigma_{1}\end{matrix}\right|=0.$$

Solving this equation we get

$$(c\lambda-D_{1}\lambda^{2}+\sigma_{1})(c\lambda+\beta)-au_{0}g_{1}(0)e^{-\lambda c\tau_{1}}=0.$$

Setting \(\mu=c\lambda\) we get the expression for the wave speed, for which there exists a positive solution for system (5)–(6), decaying at infinity

$$c^{2}=\min_{\mu>\mu_{0}}\frac{D_{1}\mu^{2}(\mu+\beta)}{(\mu+\sigma_{1})(\mu+\beta)-au_{0}g_{1}(0)e^{-\mu\tau_{1}}},$$
(8)

where \(\mu_{0}>0\) is the value for which the denominator vanishes. Since \(g_{1}(0)=b_{1}\), formula (8) is the same as for the model without interferon [3]. Therefore, we can conclude that interferon does not change the wave speed in the linearized model. This conclusion was confirmed by numerical simulations of system (1)–(4).

We obtain the following result.

Proposition 1. If there exists a travelling wave solution of system (1)–(4), then the minimal speed \(c_{0}\) satisfies the inequality \(c_{0}\geq c\), where \(c\) is given by equality (8).

The wave existence is studied in [4] in a particular case without interferon and with \(\beta=0\).

2.2 Interferon Decreases Viral Load

If the coefficients \(\beta\) and \(\sigma_{1}\) in system (1)–(4) are positive, that is, the death rates of infected cells and virus are nonzero, then the total viral load \(\int_{0}^{L}V(x,t)dx\) tends to a stationary value (Fig. 1, lower panel). We use numerical simulation to determine the value of the total viral load and its dependence on the interferon-related parameters. An example of such simulation is shown in Fig. 1, and the results are presented in Fig. 2. Increase of interferon production rate leads to the decrease of total viral load (circles), while increasing rate of interferon degradation leads to increasing viral load (pluses). Parameter \(k_{2}\) in (4) determines how intensively virus slows down the interferon production. Increasing \(k_{2}\) leads to increasing viral load (asterisks). Parameter \(k_{1}\) determines down-regulation of virus replication by interferon. Increasing \(k_{1}\) results in decreasing the total viral load (crosses) with the dependence similar to interferon production (circles). Thus, the total viral load inversely depends on the interferon amount. The viral propagation speed equals 0.029 cm/hour in all these simulations. Let us note that interferon decreases the total viral load, (cf. Fig. 2, \(b_{2}=0\) corresponds to the case without interferon), while the viral propagation speed remains the same.

Fig. 1
figure 1

Example of numerical simulations of system (1)–(4). (Upper left panel) Virus concentration distribution in space for different moments of time. (Upper right panel) Interferon concentration distribution in space for different moments of time. (Lower panel) The total viral load as a function of time. The values of parameters are as follows: \(a=0.01\) mL/(h copy), \(b_{1}=80\,000\) copy/(h cell), \(b_{2}=13.5\) unit/(h cell), \(\beta=0.1\) h\({}^{-1}\), \(\sigma_{1}=0.1\) h\({}^{-1}\), \(\sigma_{2}=3.5\) h\({}^{-1}\), \(D_{1}=0.001\) cm\({}^{2}\)/h, \(D_{2}=0.16\) cm\({}^{2}\)/h, \(k_{1}=1\) mL/unit, \(k_{2}=10^{-5}\) mL/copy, \(\tau_{1}=10\) h, \(\tau_{2}=5\) h, \(L=10\) cm and initial values \(U_{0}=1\) cell/mL, \(I_{0}=0\) cell/mL, \(V_{0}=5000\) copy/mL, \(C_{0}=0\) unit/mL. Final time \(T=200\) h, the space step \(\delta x=0.1\) cm.

Full size image
Fig. 2
figure 2

Dependence of the total viral load on interferon production rate \(b_{2}\) (circles), interferon degradation rate \(\sigma_{2}\) (pluses), intensity of slowing down the interferon production by virus \(k_{2}\) (asterisks), intensity of virus replication down-regulation by interferon \(k_{1}\) (crosses). Interferon parameters are normalized (i.e., divided) by the corresponding values from Fig. 1. Other parameters are the same as in Fig. 1.

Full size image

3 CIRCULATING INTERFERON IN THE ORGANISM

3.1 Model

In the case of infection progression in human organism, interferon is produced not only locally in the infected tissue but also elsewhere. It is redistributed in the body by global blood circulation. In this case we consider global concentration of circulating interferon without its space distribution in the infected tissue

$$\frac{\partial U}{\partial t}=-aUV,$$
(9)
$$\frac{\partial I}{\partial t}=aUV-\beta I,$$
(10)
$$\frac{\partial V}{\partial t}=D_{1}\frac{\partial^{2}V}{\partial x^{2}}+\frac{b_{1}}{1+k_{1}Z}\;I_{\tau_{1}}-\sigma_{1}V,$$
(11)
$$\frac{dZ}{dt}=b_{2}J(I)e^{-k_{2}J(V)}-\sigma_{2}Z.$$
(12)

Here

$$J(I)=\int\limits_{-\infty}^{\infty}I(x,t)dx,\quad J(V)=\int\limits_{-\infty}^{\infty}V(x,t)dx.$$

Interferon concentration \(Z\) depends on time but not on the space variable. Its production is determined by the total quantity of infected cells \(J(I)\). It is down-regulated by virus depending on its total concentration. Let us note that interferon can be produced by macrophages and other cells outside the infected tissue, and virus can also be presented in other tissues and organs. We suppose that both of them are proportional to the corresponding quantities in the tissue. In this model, the units for the state variables \(U\), \(I\), \(V\), and \(Z\) remain the same as before.

3.2 Wave Speed and Viral Load

For the steady propagating wave we obtain the following equations:

$$cu^{\prime}-auv=0,$$
$$cw^{\prime}+auv-\beta w=0,$$
$$Dv^{\prime\prime}+cv^{\prime}+\frac{b_{1}}{1+k_{1}z}\;w(\xi+c\tau)-\sigma_{1}v=0,$$
$$b_{2}J(w)e^{-k_{2}J(v)}-\sigma_{2}z=0$$
(13)

considered on the whole real axis. Here \(\xi=x-ct\),

$$J(w)=\int\limits_{-\infty}^{\infty}w(\xi)d\xi,\quad J(v)=\int\limits_{-\infty}^{\infty}v(\xi)d\xi.$$

Then,

$$cu^{\prime}-auv=0,$$
(14)
$$cw^{\prime}+auv-\beta w=0,$$
(15)
$$Dv^{\prime\prime}+cv^{\prime}+b(J)w(x+c\tau)-\sigma_{1}v=0,$$
(16)

where

$$b(J)=\frac{b_{1}}{1+k_{1}b_{2}J(w)e^{-k_{2}J(v)}/\sigma_{2}}.$$

Integrating equation (16) and taking into account zero conditions for \(v\) at infinity, we obtain \(b(J)=\sigma_{1}J(v)/J(w).\) Substituting this expression into the previous equality, we obtain

$$\sigma_{1}J(v)=\frac{b_{1}J(w)}{1+k_{1}b_{2}J(w)e^{-k_{2}J(v)}/\sigma_{2}}.$$

Hence

$$J(w)=\frac{\sigma_{1}J(v)}{b_{1}-\alpha J(v)e^{-k_{2}J(v)}},\quad\alpha=\frac{k_{1}b_{2}\sigma_{1}}{\sigma_{2}},$$

and

$$b(J)=\frac{b_{1}(b_{1}-\alpha J(v)e^{-k_{2}J(v)})}{b_{1}-\alpha J(v)e^{-k_{2}J(v)}+k_{1}b_{2}\sigma_{1}J(v)e^{-k_{2}J(v)}/\sigma_{2}}\;=b_{1}-\alpha J(v)e^{-k_{2}J(v)}.$$
(17)

On the other hand, applying previously obtained results to system (14)–(16), we get for large virus replication number \(R_{v}=abu_{0}/\beta\sigma\) [3]:

$$J(v)=\frac{cu_{0}}{\beta\sigma_{1}}\;b(J),$$

where \(c\) is the wave speed. Together with (17), this gives the equation with respect to \(J(v)\):

$$J(v)=\frac{cu_{0}}{\beta\sigma_{1}}\;\left(b_{1}-\alpha J(v)e^{-k_{2}J(v)}\right).$$
(18)

Let us recall that \(c\) depends on \(b(J)\),

$$c^{2}=\min_{\mu>\mu_{0}}\frac{D\mu^{2}(\mu+\beta)}{(\mu+\beta)(\mu+\sigma_{1})-ab(J)u_{0}e^{-\tau\mu}},$$
(19)

where \(b(J)\) is determined in (17). Substituting (19) into (18), we obtain an equation with respect to the viral load \(J(v)\). Next, having found \(J(v)\), we can determine the wave speed \(c\). However, this dependence is weak. In the first approximation, we can consider equation (18) with a constant \(c\). Depending on parameters, this equation can have from one to three solutions. Therefore, problem (14)–(16) also has from one to three solutions.

We obtain the following result.

Proposition 2. If there exists a travelling wave solution of system (9)–(12), then the wave speed and total viral load satisfy equality (18). The minimal wave speed \(c_{0}\) satisfies the inequality \(c_{0}\geq c\), where \(c\) is given by equality (19).

Equation (18) can have from one to three solutions, thus, there can be from one to three travelling wave solutions of system (9)–(12) depending on the values of parameters.

3.3 Simulation and Examples

An example of numerical simulations of system (9)–(12) is shown in Fig. 3. As previously shown, infection spreads as a reaction-diffusion wave. Concentration of circulating interferon \(Z(t)\) converges to some constant value after an initial outbreak.

Fig. 3
figure 3

Numerical simulations of system (9)–(12) with the function \(V(x,t)\) shown every 10 time units (left) and the function \(Z(t)\) (right). The values of parameters: \(a=0.01\), \(b_{1}=80\,000\), \(b_{2}=13.5\) unit/(h cm cell), \(k_{1}=1\), \(k_{2}=10^{-5}\) cm\({}^{2}\)/copy, \(D=10^{-3}\), \(\beta=0.1\), \(\sigma_{1}=0.1\), \(\sigma_{2}=3.5\), \(\tau=10\), \(U_{0}=1\), \(V_{0}=5000\), \(x_{0}=0.1\). Other units are the same as at Fig. 1.

Full size image

In order to compare numerical and analytical results, set \(y=k_{2}J(v)\) and write equation (18) in the form

$$py=q-\alpha ye^{-ry},\quad p=\frac{\beta\sigma_{1}}{c},\quad q=b_{1}k_{2},\quad r=1.$$
(20)

Taking the value of the wave speed \(c=0.03\) from numerical simulations, we determine the values of parameters in equation (20) and find \(y\) as a solution of this equation (Fig. 4, left). It gives \(y=2.09\) (intersection of red and green curves). Therefore, \(J(v)=2.09\times 10^{5}\). In numerical simulations, we find a close value \(J(v)=2.07\times 10^{5}\). We can now determine the analytical value of \(z\) from equation (13), \(z=0.15\). It coincides with the numerical value in Fig. 3 (right) for large time when the wave is established.

Fig. 4
figure 4

Graphical solution of Eq. (20) (left) for the values of parameters \(p=1/3\), \(q=0.8\), \(\alpha=0.4\), \(r=1\), left-hand side of this equation is shown with red line and right-hand side with green curve. Blue curve is the right-hand side for \(\alpha=10\), \(r=4\). The minimal wave speed is determined as the square root of the minimum of the function in the right-hand side of Eq. (19) (right).

Full size image

Next, we determine \(b(J)\) from (17) and the value of the wave speed from (19) (Fig. 4, right). As in the numerical simulations we obtain \(c=0.03\).

Let us now consider equation (20) for \(\alpha=10\) and \(r=4\) (blue line in Fig. 4, left). This equation has three solutions. Numerical simulations of system (9)–(12) confirm the existence of two stable waves corresponding to the minimal and the maximal solutions.

Interferon production decreases the viral propagation speed and the total viral load, while the degradation of interferon and increased down-regulation of it by virus result in the increasing of the viral propagation speed and the total viral load (Fig. 5).

Fig. 5
figure 5

The influence of interferon-related parameters on the viral propagation speed (upper left), total viral load (upper right), and interferon (lower left) in the model (9)–(12). The lower right figure shows the dynamics of global interferon depending on production of interferon. The values of parameters are \(a=0.01\), \(b_{1}=80\,000\), \(\beta=0.1\), \(\sigma=0.1\), \(D_{1}=0.001\), \(\tau_{1}=10\), \(u_{0}=1\). In the horizontal axes, there are parameters from the legend normalized by \(b_{2}=13.5\), \(\sigma_{2}=3.5\), \(k_{2}=1.0\times 10^{-5}\), \(k_{1}=1\). The units are the same as in Fig. 3. Dots show the results of numerical simulation and lines are analytical values obtained with formulas (18), (19).

Full size image

The interferon amount in the organism increases with increasing of its production rate, and decreases with increased degradation and down-regulation of interferon production by virus (Fig. 5, lower panel).

Parameter \(k_{1}\) which describes the rate of virus production by infected cells has almost no effect on the viral propagation speed, total viral load, and interferon amount in the chosen set of parameters (Fig. 5). However, as the virus production rate inversely depends on it, increasing of this parameter decreases the viral propagation speed and total viral load, and increases the interferon amount (Fig. 5). All these results are in agreement with the numerical simulations (Fig. 5, dots).

4 ADAPTIVE IMMUNE RESPONSE

4.1 Model

Some time after infection onset, human organism initiates the adaptive immune response based on antigen specific T lymphocytes and antibodies produced by B lymphocytes. Since antibodies and CTL are produced in lymphoid organs and not in the infected tissue, we take into account global circulation of these elements of the adaptive immune response expressed by the integral terms

$$\frac{\partial U}{\partial t}=-aUV,$$
(21)
$$\frac{\partial I}{\partial t}=aUV-\beta(T)I,$$
(22)
$$\frac{\partial V}{\partial t}=D_{1}\frac{\partial^{2}V}{\partial x^{2}}+\frac{b_{1}}{1+k_{1}Z}\;I_{\tau_{1}}-\sigma_{1}(A)V,$$
(23)
$$\frac{dZ}{dt}=b_{2}J(I)e^{-k_{2}J(V)}-\sigma_{2}Z,$$
(24)
$$\frac{dT}{dt}=b_{3}J(V)-\sigma_{3}TJ(I),$$
(25)
$$\frac{dA}{dt}=b_{4}J(V)-\sigma_{4}AJ(V).$$
(26)

Here \(T\) is the concentration of cytotoxic T-lymphocytes (CTL), and \(A\) is the concentration of antibodies. In equation (22), there is an additional term characterizing elimination of infected cells by CTL, and in equation (23) the term describing neutralization of virus by antibodies. In equation (25) for the concentration of CTL, the first term in the right-hand side describes the stimulation of their production by the total virus quantity, while the second term corresponds to the decrease of CTL concentration due to the loss of their cytotoxic potential depleted in the elimination of infected cells. The last equation describes the concentration of antibodies with the first term in the right-hand side signifying their production stimulated by the antigen through antigen-specific B-lymphocytes, and the second term characterizes depletion of antibodies proportional to the total virus quantity. In this model, the concentration of CTL is measured in cell/mL, and the concentration of antibodies is measured in unit/mL.

4.2 Wave Propagation

Writing system (21)–(26) in the moving coordinate frame, \(\xi=x-ct\), we obtain the following system of equations

$$cu^{\prime}-auv=0,$$
(27)
$$cw^{\prime}+auv-\beta(\theta)w=0,$$
(28)
$$D_{1}v^{\prime\prime}+cv^{\prime}+\frac{b_{1}}{1+k_{1}z}w(\xi+c\tau_{1})-\sigma_{1}(y)v=0,$$
(29)
$$b_{2}J(w)e^{-k_{2}J(v)}=\sigma_{2}z,$$
(30)
$$b_{3}J(v)=\sigma_{3}\theta J(w),$$
(31)
$$b_{4}J(v)=\sigma_{4}yJ(v),$$
(32)

where \(U(x,t)=u(x-ct)=u(\xi)\), \(V(x,t)=v(x-ct)=v(\xi)\), \(I(x,t)=w(x-ct)=w(\xi)\), \(Z(t)=z(t)\), \(T(t)=\theta(t)\), \(A(t)=y(t)\), \(J(w)=\int_{-\infty}^{+\infty}w(\xi)\mathrm{d}\xi\), \(J(v)=\int_{-\infty}^{+\infty}v(\xi)\mathrm{d}\xi\). From the last equation

$$y=\frac{b_{4}}{\sigma_{4}}.$$

Hence,

$$\sigma_{1}(y)=\sigma_{1}\left(\frac{b_{4}}{\sigma_{4}}\right)=\sigma_{1}(=const).$$
(33)

Following the method from the previous section, we express \(z\) from the equation (30) and substitute it in (29):

$$D_{1}v^{\prime\prime}+cv^{\prime}+b(J)w(\xi+c\tau_{1})-\sigma_{1}v=0,$$
(34)

where

$$b(J)=\frac{b_{1}}{1+k_{1}b_{2}J(w)e^{-k_{2}J(v)}/\sigma_{2}}.$$
(35)

Integrate equation (34)

$$b(J)J(w)-\sigma_{1}(y)J(v)=0,$$

and get

$$b(J)=\frac{\sigma_{1}(y)J(v)}{J(w)},$$

where \(\sigma_{1}\) is defined in (33). We substitute this equality into in (35) and we get the following expression for \(J(w)\):

$$J(w)=\frac{\sigma_{1}(y)J(v)}{b_{1}-\alpha(y)J(v)e^{-k_{2}J(v)}},\quad\alpha(y)=\frac{k_{1}b_{2}\sigma_{1}(y)}{\sigma_{2}},$$
(36)

and for \(b(J)\):

$$b(J)=\frac{b_{1}\left(b_{1}-\alpha(y)J(v)e^{-k_{2}J(v)}\right)}{b_{1}-\alpha(y)J(v)e^{-k_{2}J(v)}+k_{1}b_{2}\sigma_{1}(y)J(v)e^{-k_{2}J(v)}/\sigma_{2}}=b_{1}-\alpha(y)J(v)e^{-k_{2}J(v)}.$$
(37)

Applying previously obtained results [1] for the system (27)–(32), we get for the total viral load

$$J(v)=\frac{cu_{0}}{\beta(\theta)\sigma_{1}(y)}b(J).$$
(38)

From equation (31) with respect to \(J(w)\) we can express \(\theta\) through \(J(v)\):

$$\theta=\frac{b_{3}J(v)}{\sigma_{3}J(w)}=\frac{b_{3}\left(b_{1}-\alpha(y)J(v)e^{-k_{2}J(v)}\right)}{\sigma_{3}\sigma_{1}(y)}.$$
(39)

Substituting \(b(J)\) into (38), we get

$$J(v)=\frac{cu_{0}}{\beta(\theta)\sigma_{1}(y)}b(J)=\frac{cu_{0}}{\beta(\theta)\sigma_{1}(y)}\left(b_{1}-\alpha(y)J(v)e^{-k_{2}J(v)}\right),$$
(40)

where \(\beta(\theta)=\beta_{0}+\beta_{1}\theta\), \(\sigma(y)=\sigma_{10}+\sigma_{11}y\), \(y=\frac{b_{4}}{\sigma_{4}}\), \(c\) depends on \(b(J)\) as defined in (19), and \(\theta\) is determined in (39).

Denote \(q=k_{2}J(v)\) and rewrite equation (40) as follows

$$p_{1}q=b_{1}+\frac{\alpha(y)}{k_{2}}qe^{-q}\left(p_{2}q-1\right),$$
(41)

where

$$p_{1}=\frac{\sigma_{1}(y)\beta_{0}+\beta_{1}b_{1}b_{3}/\sigma_{3}}{k_{2}cu_{0}},\quad p_{2}=\frac{\beta_{1}b_{3}}{\sigma_{3}k_{2}cu_{0}}.$$

This is equation with respect to \(q=k_{2}J(v)\) with a linear function in the left-hand side and a decreasing exponential multiplied by the second-order polynomial in the right-hand side. Denote the latter by \(F(q)\),

$$F(q)=b_{1}+\frac{\alpha(y)}{k_{2}}qe^{-q}\left(p_{2}q-1\right).$$

We have \(F(q)\to+\infty\) for \(q\to-\infty\), \(F(q)\to 0\) for \(q\to+\infty\), and \(F(q)\) intersects the vertical axis at \(b_{1}\). The function has minimum at \(q=q_{1}\) and maximum at \(q_{2}\), the inflection points are \(q_{3}\) and \(q_{4}\), where

$$q_{1}=\frac{1}{2}\left(2+1/p_{2}-\sqrt{4+1/p_{2}^{2}}\right),\quad q_{2}=\frac{1}{2}\left(2+1/p_{2}+\sqrt{4+1/p_{2}^{2}}\right),$$
$$q_{3}=\frac{1}{2}\left(4+1/p_{2}-\sqrt{14+6/p_{2}+1/p_{2}^{2}}\right),\quad q_{4}=\frac{1}{2}\left(4+1/p_{2}+\sqrt{14+6/p_{2}+1/p_{2}^{2}}\right),$$

and \(q_{1}<q_{3}<q_{2}<q_{4}\). Using these results, we can conclude that equation (41) has from one to three solutions (Fig. 6, left).

Fig. 6
figure 6

Left: graphical solution of equation (41). The blue curve represents the right-hand side of this equation, and other curves its left-hand side for different values of parameters. Depending on parameters, this equation can have from one to three solutions. Right: the plot of the right-hand side of equation (42).

Full size image

The wave speed is given by the formula

$$c^{2}=\min_{\mu>\mu_{0}}\frac{D\mu^{2}(\mu+\beta(\theta))}{(\mu+\beta(\theta))(\mu+\sigma_{1}(y))-ab(J)u_{0}e^{-\tau\mu}},$$
(42)

where

$$\beta(\theta)=\beta_{0}+\beta_{1}\theta,\quad\theta=\frac{b_{3}J(v)}{\sigma_{3}J(w)},\quad\sigma_{1}(y)=\sigma_{10}+\sigma_{11}y,\quad y=\frac{b_{4}}{\sigma_{4}}.$$

A typical form of right-hand side of this equation as a function of \(\mu\) is presented in Fig. 6 (right). It has one minimum for \(\mu>\mu_{0}>0\), where \(\mu=\mu_{0}\) is the zero of the denominator.

We obtain the following result.

Proposition 3. If there exists a travelling wave solution of system (21)–(26), then the wave speed and total viral load satisfy equality (40). The minimal wave speed \(c_{0}\) satisfies the inequality \(c_{0}\geq c\), where \(c\) is given by equality (42).

Equation (40) can have from one to three solutions, thus, there can be from one to three traveling wave solutions of system (21)–(26) depending on the values of parameters.

We use these results to calculate the total viral load, wave speed and interferon level as follows. First, we determine the total viral load \(J(v)\) and the coefficient \(b(J)\) with formulas (41) and (37), respectively. Then, we find the analytical value of the minimal wave speed with formula (42). In (41), the value of \(c\) is already required, so we use an initial approximation for it, and then correct it if this approximation does not match the wave speed determined with (42).

The results of numerical simulation are presented in Fig. 7. In the calculations, the interferon concentration converges to the value \(0.091\), and it corresponds to the analytical value \(0.09\) given by formulas (24), (36), (40). The total viral load in numerical simulation is \(133.5\), and it well corresponds to \(132.7\) obtained in the analytical solution with formula (40). The wave speed equals \(0.0181\) both in numerical simulation and in the analytical solution.

Fig. 7
figure 7

Numerical solution of system (21)–(26). Left: virus concentration distribution in space in different moments of time. Right: interferon concentration as a function of time. The values of parameters are as follows: \(a=0.01\), \(\beta_{0}=0.1\) h\({}^{-1}\), \(\beta_{1}=0.1\times 10^{-5}\) mL/(h cell), \(D=10^{-3}\), \(b_{1}=80\,000\), \(b_{2}=350\), \(b_{3}=13.5\) cell/(h cm copy), \(b_{4}=350\) unit/(h cm copy), \(\sigma_{10}=0.1\) h\({}^{-1}\), \(\sigma_{11}=0.004\) mL/(h unit), \(\sigma_{2}=3.5\), \(\sigma_{3}=0.1\) cm\({}^{2}\)/(h cell), \(\sigma_{4}=3.5\) cm\({}^{2}\)/(h copy), \(k_{1}=1\), \(k_{2}=10^{-5}\), \(\tau=10\), \(U_{0}=1\). Other units are the same as in Fig. 1.

Full size image

The dependence of the viral propagation speed, total viral load, and interferon concentration on the parameters of the immune response is presented in Fig. 8, and the results of its analysis are given in Table 1. Virus neutralization by antibodies (dotted yellow line) has the same effect on all controlled values (viral propagation speed, total viral load, interferon) as the production of antibodies (solid blue line), and elimination of infected cells by CTL (dotted purple line) has the same effect as the production of CTL (solid red line). It can be explained by the linear dependence of degradation (death) terms on antibodies and CTL concentrations in equations (23) and (22), respectively.

Fig. 8
figure 8

Dependence of the wave speed (upper left), total viral load (upper right), and interferon concentration (lower panel) on the parameters of the immune response. Blue lines correspond to the production (solid) and death (dashed) of antibodies, red lines correspond to the production (solid) and death (dashed) of CTL. Dotted lines represent the elimination of the infected cells by CTL (purple) and virus neutralization by antibodies (yellow). Markers show the corresponding numerical results: circles for the production of antibodies (blue) and CTL (red), stars for the death of antibodies (blue) and CTL (red), pluses for the neutralization of virus by antibodies (blue) and elimination of infected cells by CTL (red). In the horizontal axes, there are parameters from legend normalized on corresponding values, \(b_{4}=350\), \(\sigma_{4}=3.5\), \(b_{3}=13.5\), \(\sigma_{3}=0.1\), \(\beta_{1}=0.1\times 10^{-5}\), \(\sigma_{11}=0.004\). Other parameters are listed in Fig. 7.

Full size image
Table 1. Influence of the immune response parameters on the viral propagation speed, total viral load, and interferon concentration. The results are presented for the increasing of listed parameters. The corresponding graphs are presented in Fig. 8
Full size table

5 CONCLUSIONS

Infection progression in cell culture and tissue is characterized by the spreading speed and viral load. The former determines degree of tissue damage and correlates with the severity of symptoms while the latter determines infection transmission rate between different individuals. In this work we study how these characteristics are influenced by the immune response.

The model considered in Section 2 takes into account the reduction of virus replication as a result of local interferon production. The modeling results show that this process does not change the spreading speed, but the total viral load decreases. In the case of globally circulating interferon considered in Section 3, the results are different. In this case, both the spreading speed and the total viral load decrease with increase of global interferon level. Figure 5 (lower right panel) shows essential increase of the global interferon level. It has positive effect since the viral propagation speed and total viral load can be further decreased. However, interferon has its negative side effects as inhibiting cell proliferation, so its concentration should be controlled. This limitation of the model can be addressed in future works.

Simple model of the adaptive immune response in Section 4 allows investigating the system not only numerically, but also analytically. The influence of the intensity of the adaptive immune response on viral propagation speed and total viral load is evaluated. It is shown that more intensive immune response results in reduced viral propagation speed and total viral load. From these results we can conclude that intensive immune response decreases both infectivity of virus and severity of symptoms, i.e. the total viral load and the viral propagation speed respectively.

Modeling results show that different parts of the adaptive immune response effect the final concentration of interferon in different ways. As such, an increase in the amount of antibodies correlates with an increase in total interferon concentration, and an increase in the amount of CTL correlates with its decrease. Taking into account the previous remark regarding the possible side effects of interferon, we can conclude that the model allows the choice of the preferred treatment mechanism to reduce the infectiousness and severity of symptoms of the disease without increasing the global level of interferon. In this model, such a treatment strategy will be provided by an increase in CTL level.