# Viral Infection Spreading and Mutation in Cell Culture

Latifa Ait Mahiout, Bogdan Kazmierczak, Vitaly Volpert

► To cite this version:

Latifa Ait Mahiout, Bogdan Kazmierczak, Vitaly Volpert. Viral Infection Spreading and Mutation in Cell Culture. Mathematics, 2022, 10, 10.3390/math10020256. hal-03840321

HAL Id: hal-03840321

https://hal.science/hal-03840321

Submitted on 5 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Viral Infection Spreading and Mutation in Cell Culture

Latifa Ait Mahiout<fnref n="1" />, Bogdan Kazmierczak<fnref n="2" /> and Vitaly Volpert<fnref n="3" />,<fnref n="4" />,<fnref n="5" />,*

**Abstract:** A new model of viral infection spreading in cell cultures is proposed taking into account virus mutation. This model represents a reaction-diffusion system of equations with time delay for the concentrations of uninfected cells, infected cells and viral load. Infection progression is characterized by the virus replication number $R_v$, which determines the total viral load. Analytical formulas for the speed of propagation and for the viral load are obtained and confirmed by numerical simulations. It is shown that virus mutation leads to the emergence of a new virus variant. Conditions of the coexistence of the two variants or competitive exclusion of one of them are found, and different stages of infection progression are identified.

**Keywords:** viral infection; mutation; cell culture; reaction-diffusion equations; time delay

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Ait Mahiout, L.; Kazmierczak, B.; Volpert, V. Viral Infection Spreading and Mutation in Cell Culture. Mathematics **2022**, *10*, 256. https://doi.org/10.3390/math10020256

Academic Editor: Wantong Li

Received: 8 December 2021

Accepted: 11 January 2022

Published: 14 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Viruses constitute a separate specific kingdom of entities. They can have a plethora of geometrical forms, characterized by different shapes, symmetries and sizes, which cover the interval ranging from 20 nm (Porcine circovirus) up to 700 nm (Mimivirus). These properties depend mainly on the amount of nucleic acids they contain, their forms (single stranded positive- or negative-sense RNA, double stranded RNA, single stranded DNA or double stranded DNA), as well as on the physiology of the type of cells they use for their proliferation. Despite these essential differences, the survival and proliferation of all viruses (as quasi-species) is similar and based on the following strategy: enter the host cells, force the cells to multiplicate their genetic material (thus, in fact, to produce their multiple copies), leave the cells and invade new ones.

Experimental or clinical assessment of the progression of viral infection implies the evaluation of virus concentration in the infected tissue by means of conventional multiplicity of infection (MOI) assays (see, e.g., [1–5] and the references therein). After several consecutive dilutions, the virus-containing solution is poured in a cell culture leading to the formation of virus plaques. The number of such plaques determines the virus concentration measured in plaque forming units (PFU). The plaques consist of dead or modified cells, and they form due to virus replication inside the cells and its random motion (diffusion) between cells. The plaque growth rate characterizes the efficacy of virus penetration inside the cells, the intensity of its production and transmission between cells, thus, the virus virulence.

In spite of the importance and wide use of these experiments, there are relatively few modeling works devoted to the infection progression in cell culture taking into account spatial distributions of cells and virus particles. A reaction-diffusion model of plaque growth is considered in [6] in the case of reversible host infection. The plaque growth rate is

<footnote n="1">¹ Laboratoire D'équations aux Dérivées Partielles non Linéaires et Histoire des Mathématiques, Ecole Normale Supérieure, B.P. 92, Vieux Kouba, Algiers 16050, Algeria; latifa.aitmahiout@g.ens-kouba.dz</footnote>

<footnote n="2">² Institute of Fundamental Technological Research Polish Academy of Sciences, 02-106 Warsaw, Poland; bkazmier@ippt.pan.pl or bkazmier@brama.ippt.pan.pl</footnote>

<footnote n="3">³ Institut Camille Jordan, UMR 5208 CNRS, University Lyon 1, 69622 Villeurbanne, France;</footnote>

<footnote n="4">⁴ INRIA Team Dracula, INRIA Lyon La Doua, 69603 Villeurbanne, France</footnote>

<footnote n="5">5 S.M. Nikolskii Mathematical Institute, Peoples Friendship University of Russia (RUDN University), 6 Miklukho-Maklaya St., 117198 Moscow, Russia</footnote>

<footnote>* Correspondence: volpert@math.univ-lyon1.fr</footnote>
determined by the method of linearization (cf. Section 3.2). Numerical simulations of viral plaque growth described by a reaction-diffusion model with time delay are presented in [7]. This model is different in comparison with the model considered below, and these different approaches are complementary. Individual based models of viral infection spreading in cell cultures are developed in [8,9].

By that means, the region of infected cells increases and, after some preliminary time, the boundary of the infected region may propagate as a traveling wave. This fact has been proven in [10] in the framework of a model, in which the process of cell infection and virus multiplication with possible delays $\tau$ between the moment of the cell infection and the moment when the multiplied viruses are released, were taken into account. In a more complete model, the phenomenon of anti-virus defense via the interferon production by the infected cells, was considered. By the theoretical and numerical analysis presented in [10] we showed that, within the proposed model, the following three phases of the virus infection could be distinguished (in agreement with the experimental results):

1. Virus concentration decay during time delay in its replication in cells;

2. Explosive growth of concentration, when infected cells begin to reproduce new viruses;

3. Propagation of infection wave along the cell culture.

It was also found that the minimal speed of infection decreases with the delay time $\tau$. This speed was effectively determined by finding the minimum of an explicitly constructed function containing the parameters of the model.

However, the model proposed in [10] neglects other characteristic features of viruses survival strategy, namely their natural mutations and competition between their strains, which is crucial in the context of epidemiological analysis of virus infection spreading. An incorporation of these effects is the main objective of the present paper. We show that the model supplemented with the terms corresponding to mutation and competition phenomena, though significantly more complicated, can be studied by similar analytic methods and, together with the numerical calculations, can be a source of interesting information concerning the evolution of infection in the cell culture. The paper is organized as follows.

In Section 2, we analyze the model without mutation and competition terms and assuming that the mortality coefficients of the infected cells is nonzero, we show that the progression of infection is possible only if the virus replication number $R_v$ satisfies the inequality $R_v > 1$. We also derive a lower bound for the minimal speed of the infection front. In Section 3, we propose a model taking into account mutation and competition phenomena. In Section 3.1, we show that, under some assumptions on the coefficients characterizing the considered strains of the virus, the new model can be reduced to one virus model (without mutation and competition). We also construct a function determining the lower bound of the infection propagation speed to the considered case. In Section 3.3, a system of reaction-diffusion equations describing multiple mutations and coinfections is proposed. In Section 4, we analyze different stages of infection progression in terms of time dependence of the total virus loads. The values and units of the assumed parameter values and their reference to culture experiments is summarized at the end of Section 4. Numerical simulations presented in this work are carried out using the finite element method with software freeFem++ [11]. Numerical accuracy of the simulations is controlled by decreasing the time and space discretization, and by the comparison with the analytical results.

## 2. Infection Progression with Single Virus

### 2.1. Virus Replication Number and Final Size of Infection

#### 2.1.1. Virus Replication Number

We begin with a conventional model of viral infection progression in cell culture without taking into account its spatial distribution:
$$ \frac{dU}{dt} = -aUV, \qquad (1) $$

$$ \frac{dI}{dt} = aUV - \beta I, \qquad (2) $$

$$ \frac{dV}{dt} = bI - \sigma V. \qquad (3) $$

Here, $U(t)$ is the concentration of uninfected cells, $I(t)$ of infected cells and $V(t)$ of virus; $a, b, \beta, \sigma$ are positive parameters. Parameter $a$ characterizes infection rate of uninfected cells, $b$ is the rate of virus replication in infected cells, $\beta$ is the death rate of infected cells, and $\sigma$ is the virus death rate. We suppose that

$$ U(0) = U_0, \quad I(0) = 0, \quad V(0) = V_0. \qquad (4) $$

Clearly, $U = U_0, I = 0, V = 0$ is a stationary point of system (1)–(3). If $\beta > 0$, then it is a unique stationary point. In order to study its stability, we linearize the system about this stationary point and consider the corresponding eigenvalue problem:

$$ \begin{pmatrix} -\beta & aU_0 \\ b & -\sigma \end{pmatrix} \begin{pmatrix} I \\ V \end{pmatrix} = \lambda \begin{pmatrix} I \\ V \end{pmatrix}. $$

Denote the matrix in the left-hand side of this equality by $A$. Condition det $A = 0$ can be written as $R_v = 1$, where

$$ R_v = \frac{abU_0}{\beta\sigma}. $$

By analogy with epidemiological models, we will call it virus replication number (VRN). If $R_v < 1$, then both eigenvalues of the matrix $A$ are negative, and the stationary point is stable. If $R_v > 1$, then this matrix has one positive eigenvalue, and the stationary point is unstable. In the first case, viral load $V(t)$ decays, in the second case it grows.

VRN represents a ratio of virus replication and elimination rates. Considering cell culture, we interpret parameters $\beta$ and $\sigma$ as death rates of infected cells and viruses. In the application to living tissue, which will be considered in the subsequent works, these parameters characterize the strength of the immune response with the elimination of infected cells by cytotoxic lymphocytes and virus neutralization by antibodies.

### 2.1.2. Final Concentrations of Cells and Total Viral Load

Next, assuming that $R_v > 1$, we will determine the final value $U_f$ of uninfected cells as $t \to \infty$. Taking a sum of Equations (1) and (2) and integrating from 0 to $\infty$, we obtain

$$ U_0 - U_f = \beta \int_0^\infty I(t)dt. \qquad (5) $$

Integrating Equation (3), we have

$$ V_0 = \sigma \int_0^\infty V(t)dt - b \int_0^\infty I(t)dt. \qquad (6) $$

Finally, we divide Equation (1) by $U$ and integrate:

$$ \ln \frac{U_f}{U_0} = -a \int_0^\infty V(t)dt. \qquad (7) $$

Excluding the integrals from Equations (5)–(7), we deduce the following equation

$$ R_v(U - 1 - V_0) = \ln U, \qquad (8) $$
where

$$U = \frac{U_f}{U_0}, \quad V_0 = \frac{\beta V_0}{bU_0}.$$

The derivation of Equation (8) is inspired by the derivation of the final size of epidemic in the model SIR, and the equation is the same, though the model is different.

If we neglect the initial viral load and set $V_0 = 0$, then Equation (8) has a positive solution in the interval $0 < U < 1$ if and only if $R_v > 1$. This case corresponds to infection progression, and the final value of uninfected cells is less than its initial value. If $V_0 > 0$, then such solution exists for all values of $R_v$.

Having found the solution of Equation (8), we can determine the total number of infected cells, $I_f = U_0 - U_f$ and the total viral load

$$V_T \equiv \int_{0}^{\infty} V(t)dt = -\frac{1}{a}\ln U.$$

### 2.1.3. Model with Time Delay

If we take into account time delay in virus replication in the infected cells, then instead of Equation (3) we consider the equation

$$\frac{dV(t)}{dt} = bI(t - \tau) - \sigma V(t), \quad (9)$$

where $\tau$ is a positive number. The initial condition for $I(t)$ is now considered as $I(t) = 0$ for $-\tau \le t \le 0$. In the stability analysis for system (1), (2), (9), we have $\det A = \beta\sigma - abU_0e^{-\lambda\tau}$. The stability boundary, that is, the values of parameters for which $\lambda = 0$, remains the same as before, $R_v = 1$. In the derivation of the final concentrations of cells, instead of Equation (6) we have now the equation

$$V_0 = \sigma \int_{0}^{\infty} V(t)dt - b \int_{0}^{\infty} I(t - \tau)dt. \quad (10)$$

Since

$$\int_{0}^{\infty} I(t - \tau)dt = \int_{-\tau}^{\infty} I(t)dt = \int_{0}^{\infty} I(t)dt,$$

then Equation (10) is reduced to Equation (6), and the formulas for the final call concentrations and total viral load remain valid for the model with time delay.

### 2.2. Spatial Infection Spreading

We continue with the model of viral infection progression in cell culture taking into account its spatial distribution [10]:

$$\frac{\partial U}{\partial t} = -aUV, \quad (11)$$

$$\frac{\partial I}{\partial t} = aUV - \beta I, \quad (12)$$

$$\frac{\partial V}{\partial t} = D \frac{\partial^2 V}{\partial x^2} + bI - \sigma V. \quad (13)$$

The first term in the right-hand side of Equation (13) describes virus diffusion in the extracellular matrix, $D$ is the diffusion coefficient.

#### 2.2.1. Wave Propagation

Experimental data and previous modeling show [10] that infection spreads in cell cultures as a reaction-diffusion wave (Figure 1). In order to study this solution, we consider a system of Equations (11)–(13) on the whole axis, and set $U(x, t) = u(x - ct)$, $I(x, t) =$
$w(x - ct), V(x, t) = v(x - ct)$, where c is the wave speed. Then we obtain the following system of equations:

$$ - cu' = -auv, \qquad (14) $$

$$ - cw' = auv - \beta w, \qquad (15) $$

$$ -cv' = Dv'' + bw - \sigma v, \qquad (16) $$

where primes denote the derivatives with respect to the variable $\zeta = x - ct$. We are interested in the solution of this system of equations with the following limits at infinity:

$$ u(-\infty) = u_f, \quad u(\infty) = u_0, \quad w(\pm\infty) = 0, \quad v(\pm\infty) = 0. \qquad (17) $$

In order to determine unknown value $u_f$, we proceed as in the previous section. From Equation (14), we obtain

$$ c \ln \frac{u_0}{u_f} = a \int_{-\infty}^{\infty} v(x) dx, \qquad (18) $$

taking a sum of (14), (15) and integrating:

$$ c(u_0 - u_f) = \beta \int_{-\infty}^{\infty} w(x) dx, \qquad (19) $$

and from (16):

$$ b \int_{-\infty}^{\infty} w(x) dx = \sigma \int_{-\infty}^{\infty} v(x) dx. \qquad (20) $$

Excluding the integrals from Equations (18)–(20), we obtain the equation

$$ R_v(\omega - 1) = \ln \omega, \qquad (21) $$

with respect to $\omega = u_f/u_0$. Equation (21) has a solution $\omega$ in the interval $0 < \omega < 1$ if and only if $R_v > 1$. Thus, we obtain the following theorem.

**Theorem 1.** Inequality $R_v > 1$ provides a necessary condition for the existence of a positive solution of problem (14)–(17). In this case, the final value $u_f$ can be found from Equation (21), and the total viral load

$$ V_X \equiv \int_{-\infty}^{\infty} v(x) dx = -\frac{c}{a} \ln \omega. $$

<figure><img src="image_4.png" /><figcaption>**Figure 1.** A snapshot of numerical solution of system (11)–(13). Infection spreads as a reaction-diffusion wave with a speed c, $U(x, t) = u(x - ct)$, $I(x, t) = w(x - ct)$, $V(x, t) = v(x - ct)$. Note that $u_f = 0$ here.</figcaption></figure>
Let us note that the notion of total viral load $V_X$ here is different from $V_T$ in the previous section. Here it signifies the total viral concentration in cell culture at any moment of time, while before it was the total viral load with respect to time (without space distribution). Furthermore, $V_X$ depends on the wave speed $c$, which is not yet found.

### 2.2.2. Model with Time Delay

In the model with time delay, Equation (13) is replaced by the equation

$$ \frac{\partial V(x,t)}{\partial t} = D \frac{\partial^2 V(x,t)}{\partial x^2} + bI(x,t - \tau) - \sigma V(x,t), \quad (22) $$

and Equation (16) by the equation

$$ -cv'(\xi) = Dv''(\xi) + bw(\xi + c\tau) - \sigma v(\xi). \quad (23) $$

Clearly, integration of this equation gives (20), and all conclusions above remain applicable for this case.

### 2.2.3. Minimal Wave Speed

In order to determine the wave speed, we will use the linearization method widely used for monostable reaction-diffusion equations beginning from the KPP work [12]. The idea of the method consists to consider the system linearized at infinity and to find the minimal wave speed $c_*$ for which this system has a positive decreasing solution. Thus, instead of Equations (15) and (16) we consider the equations

$$ cw' + au_0v - \beta w = 0, \quad (24) $$

$$ Dv'' + cv' + bw(\xi + c\tau) - \sigma v = 0, \quad (25) $$

where $u(\xi)$ is replaced by its value $u_0$ at infinity. We look for a solution of this system in the form

$$ w(\xi) = pe^{-\lambda\xi}, \quad v(\xi) = qe^{-\lambda\xi} $$

with some positive numbers $p, q$ and $\lambda$. Substituting these expressions into Equations (24) and (25) and excluding $p$ and $q$, we obtain the following equation with respect to $\lambda$:

$$ D\lambda^2 - c\lambda - \sigma + \frac{abu_0}{c\lambda + \beta} = 0. $$

We need to find a minimal positive value of $c$ for which this equation has a positive solution. Set $\mu = \lambda c$. Then from the last equation we find

$$ c^2 = F(\mu), \quad F(\mu) = \frac{D\mu^2(\mu + \beta)e^{\mu\tau}}{(\mu + \sigma)(\mu + \beta)e^{\mu\tau} - abu_0}. $$

Under the condition $R_v > 1$, that is, $\sigma\beta < abu_0$, function $F(\mu)$ becomes negative for small positive $\mu$ and tends to $-\infty$ as $\mu \nearrow \mu_0$ for a positive $\mu_0$ for which the denominator vanish. Therefore, $F(\mu) > 0$ for $\mu > \mu_0$, $F(\mu)$ decreases in some interval $\mu_0 < \mu < \mu_*$ and increases for $\mu > \mu_*$. Hence, this function has a minimum for $\mu > \mu_0$ reached at $\mu = \mu_*$. Set

$$ c_*^2 = \min_{\mu > \mu_0} F(\mu). \quad (26) $$

This equality provides the minimal wave speed found by the linearization method, and $\lambda = \mu_*/c_*$. Let us note that, in general, this method gives an estimate of the minimal wave speed $c_0$ from below, $c_0 \ge c_*$. In some cases, it can be proved that $c_0 = c_*$. This question is discussed in the next paragraph. Having found the wave speed, we can determine the total viral load $V_X$ in Theorem 1.
### 2.2.4. Wave Existence

If $\beta = 0$, then we can reduce system (11)–(13) to a system of two equations. Indeed, taking a sum of Equations (11) and (12), we conclude that $U(x, t) + I(x, t) = U_0$, and the variable $U$ can be excluded:

$$ \frac{\partial I}{\partial t} = a(U_0 - I)V - \beta I, \qquad (27) $$

$$ \frac{\partial V}{\partial t} = D \frac{\partial^2 V}{\partial x^2} + b I_{\tau} - \sigma V. \qquad (28) $$

This is a monotone system of equations for which it is possible to apply the maximum principle and various methods based on it. Existence of wave for this system is proved in [10] for all $c \ge c_0$ for the values of $\tau$ limited from above.

Thus, $c_0 = c_*$ for $\beta = 0$. If $\beta > 0$, this equality is verified in numerical simulations [10]. Knowing the value of speed, we can construct an approximate solution in the following way. Consider system (15), (16) with $u = u_0$ for $\zeta > 0$ and $u = u_f$ for $\zeta < 0$. Then we find the wave speed $c_*$ and solution $w(\zeta)$, $v(\zeta)$ for $\zeta > 0$ as described in the previous paragraph. Let us note that this solution contains one unknown constant since there is a single relation between $p$ and $q$. Next, we have a linear system for $\zeta < 0$ with a given value $c = c_*$. We can find its solution, and it depends on two unknown constants. These three constants can be determined from the continuity conditions at $\zeta = 0$: $w(-0) = w(+0)$, $v(-0) = v(+0)$, $v'(-0) = v'(+0)$. These calculations are quite cumbersome, and we do not present them here.

## 3. Viral Infection with Mutation

We will now consider the case where the first virus can mutate in another one. We consider the system of equations

$$ \frac{\partial U}{\partial t} = -a_1 UV_1 - a_2 UV_2, \qquad (29) $$

$$ \frac{\partial I_1}{\partial t} = a_1 UV_1 - \beta_1 I_1, \qquad (30) $$

$$ \frac{\partial I_2}{\partial t} = a_2 UV_2 - \beta_2 I_2, \qquad (31) $$

$$ \frac{\partial V_1}{\partial t} = D_1 \frac{\partial^2 V_1}{\partial x^2} + b_1 I_1 \tau_1 - \sigma_1 V_1, \qquad (32) $$

$$ \frac{\partial V_2}{\partial t} = D_2 \frac{\partial^2 V_2}{\partial x^2} + \epsilon I_1 \tau_1 + b_2 I_2 \tau_2 - \sigma_2 V_2 \qquad (33) $$

describing two viruses $V_1$ and $V_2$ in cell culture. Here $U$ is the concentration of uninfected cells, $I_1$ of infected cells by virus $V_1$ and $I_2$ by virus $V_2$, $I_{\tau_1}(x) = I(x, t - \tau_1)$, $I_{\tau_2}(x) = I(x, t - \tau_2)$, $\epsilon$ is the mutation rate which shows that cells infected by virus $V_1$ can produce virus $V_2$.

### 3.1. Reduction to the One-Virus Model

System (37)–(39) can be reduced to the one-virus model by setting $I_1 = 0, V_1 = 0$ or $I_2 = 0, V_2 = 0$. We are interested in the case where all components of the solution are positive. Under the assumptions $\beta_1 = \beta_2, \sigma_1 = \sigma_2, D_1 = D_2, \tau_1 = \tau_2$, we look for the solution of system (29)–(33) in the form $I_2(x, t) = k_1 I_1(x, t), V_2(x, t) = k_2 V_1(x, t)$. Then we rewrite Equations (31) and (33) as follows:

$$ \frac{\partial I_1}{\partial t} = \frac{a_2 k_2}{k_1} UV_1 - \beta_1 I_1, \qquad (34) $$

$$ \frac{\partial V_1}{\partial t} = D_1 \frac{\partial^2 V_1}{\partial x^2} + \frac{\epsilon}{k_2} I_1 \tau_1 + \frac{b_2 k_1}{k_2} I_1 \tau_1 - \sigma_1 V_1. \qquad (35) $$
Comparing these equations with Equations (30) and (32), we conclude that

$$ \frac{a_2 k_2}{k_1} = a_1, \quad \frac{\epsilon + b_2 k_1}{k_2} = b_1. $$

Hence,

$$ k_1 = \frac{\epsilon a_2}{a_1 b_1 - a_2 b_2}, \quad k_2 = \frac{\epsilon a_1}{a_1 b_1 - a_2 b_2}. \qquad (36) $$

Since we look for positive solutions, that is $k_1, k_2 > 0$, then we should impose the condition $a_1 b_1 > a_2 b_2$.

Set $I = I_1 + I_2 = (1 + k_1)I_1, V = V_1 + V_2 = (1 + k_2)V_1$. Then we can write Equation (29), sum of Equations (30) and (31), and of Equations (32) and (33) as follows:

$$ \frac{\partial U}{\partial t} = -aUV, \qquad (37) $$

$$ \frac{\partial I}{\partial t} = aUV - \beta I, \qquad (38) $$

$$ \frac{\partial V}{\partial t} = D \frac{\partial^2 V}{\partial x^2} + bI\tau_1 - \sigma V, \qquad (39) $$

where $\beta = \beta_1, \sigma = \sigma_1, D = D_1$,

$$ a = \frac{a_1 + a_2 k_2}{1 + k_2}, \quad b = \frac{\epsilon + b_1 + b_2 k_1}{1 + k_1}, $$

where $k_1$ and $k_2$ are given by equalities (36). Hence, we obtain system (11)–(13) (or (22) in the case of time delay). We can formulate the following result.

**Theorem 2.** Suppose that system (37)–(39) has a traveling wave solution $U(x, t) = u(x - ct)$, $I(x, t) = w(x - ct)$, $V(x, t) = v(x - ct)$, that is, solution of problem (14)–(17). If $abU_0 > \beta\sigma$ and $a_1 b_1 > a_2 b_2$, then system (29)–(33) has a traveling wave solution with positive components $U(x, t) = u(x - ct)$, $I_i(x, t) = w_i(x - ct)$, $V_i(x, t) = v_i(x - ct)$, $i = 1, 2$ such that

$$ w_1(\xi) = \frac{w(\xi)}{1 + k_1}, \quad w_2(\xi) = \frac{k_1 w(\xi)}{1 + k_1}, \quad v_1(\xi) = \frac{v(\xi)}{1 + k_2}, \quad v_2(\xi) = \frac{k_2 v(\xi)}{1 + k_2}. \qquad (40) $$

This theorem allows us to apply to system (29)–(33) the results of the previous section obtained for system (11)–(13) on wave existence and speed of propagation.

### Choice between Three Solutions

System (29)–(33) can have one-virus solutions for which $I_1 = 0, V_1 = 0$ or $I_2 = 0, V_2 = 0$, and two-virus solution for which all components are positive. Conditions $a_i b_i U_0 > \beta_i \sigma_i$ and the results of the previous section provide the existence of the first-virus solution for $i = 1$ and of the second-virus solution for $i = 2$. Conditions of the existence of two-virus solution are provided by Theorem 2. If these conditions are not satisfied, we can expect that the two-virus solution does not exist, though Theorem 2 does not state this non-existence result. It is confirmed by the results of numerical simulations. Analytical values of wave speed coincide with the wave speed in numerical simulations.

## 3.2. Method of Linearization

Reduction to the one-virus model is applied above under the assumptions $D_1 = D_2, \beta_1 = \beta_2, \sigma_1 = \sigma_2, \tau_1 = \tau_2$. If these conditions are not satisfied, we use the method of linearization to determine different regimes of infection spreading. We search solutions
of system (29)–(33) in the form $U(x,t) = u(x-ct), I_1(x,t) = w_1(x-ct), I_2(x,t) = w_2(x-ct), V_1(x,t) = v_1(x-ct), V_2(x,t) = v_2(x-ct)$. Then,

$$ cu' - a_1uv_1 - a_2uv_2 = 0, \quad (41) $$

$$ cw'_{1} + a_{1}uv_{1} - \beta_{1}w_{1} = 0, \quad (42) $$

$$ cw'_{2} + a_{2}uv_{2} - \beta_{2}w_{2} = 0, \quad (43) $$

$$ D_1 v_1'' + c v_1' + b_1 w_1 (\xi + c \tau_1) - \sigma_1 v_1 = 0, \quad (44) $$

$$ D_2 v_2'' + c v_2' + \epsilon w_1 (\xi + c \tau_1) + b_2 w_2 (\xi + c \tau_2) - \sigma v_2 = 0. \quad (45) $$

Set $u = u_0$,

$$ v_1(x) = p_1 e^{-\lambda x}, v_2(x) = p_2 e^{-\lambda x}, w_1(x) = q_1 e^{-\lambda x}, w_2(x) = q_2 e^{-\lambda x}. $$

Then,

$$ p_1 = \frac{\lambda c + \beta_1}{a_1 u_0} q_1, \quad p_2 = \frac{\lambda c + \beta_2}{a_2 u_0} q_2, \quad (46) $$

$$ (D_1 \lambda^2 - c\lambda - \sigma) p_1 + b_1 q_1 e^{-\lambda c \tau_1} = 0, \quad (47) $$

$$ (D_2 \lambda^2 - c\lambda - \sigma) p_2 + \epsilon q_1 e^{-\lambda c \tau_1} + b_2 q_2 e^{-\lambda c \tau_2} = 0. \quad (48) $$

After straightforward calculations, we obtain

$$ (D_1 \lambda^2 - c\lambda - \sigma)(c\lambda + \beta_1)e^{\lambda c \tau_1} q_1 = -b_1 a_1 u_0 q_1, \quad (49) $$

$$ ((D_2 \lambda^2 - c\lambda - \sigma)(\lambda c + \beta_2) + b_2 a_2 u_0 e^{-\lambda c \tau_2}) q_2 + \epsilon a_2 u_0 e^{-\lambda c \tau_1} q_1 = 0. \quad (50) $$

If $q_1 \neq 0$, then we deduce from (49) that

$$ (D_1 \lambda^2 - c\lambda - \sigma)(\lambda c + \beta_1)e^{\lambda c \tau_1} = -a_1 u_0 b_1. \quad (51) $$

Set $\mu = c\lambda$. Then from the previous equality

$$ \left(D_1 \frac{\mu^2}{c^2} - \mu - \sigma\right) (\mu + \beta_1) e^{\mu \tau_1} = -a_1 u_0 b_1. $$

The wave speed $c_0$ is given by the equality

$$ c_0^2 = \min_{\mu>0} F(\mu), \quad F(\mu) = \frac{D_1 \mu^2 (\mu + \beta_1) e^{\mu \tau_1}}{(\mu + \sigma)(\mu + \beta_1) e^{\mu \tau_1} - b_1 a_1 u_0}. \quad (52) $$

In the particular case $D_1 = D_2(= D)$, $\tau_1 = \tau_2(= \tau)$, $\beta_1 = \beta_2(= \beta)$, we obtained from (50) and (51), (46):

$$ \frac{q_2}{q_1} = \frac{\epsilon a_2}{b_1 a_1 - b_2 a_2}, \quad \frac{p_2}{p_1} = \frac{a_1 q_2}{a_2 q_1}. \quad (53) $$

If $q_1 = 0$, then we have from (50)

$$ ((D_2 \lambda^2 - c\lambda - \sigma)(\lambda c + \beta_2) + b_2 a_2 u_0 e^{-\lambda c \tau_2}) q_2 = 0. \quad (54) $$

Then,

$$ (D_2 \lambda^2 - c\lambda - \sigma)(\lambda c + \beta_2)e^{\lambda c \tau_2} = -b_2 a_2 u_0 \quad (55) $$

Set $\mu = \lambda c$. The wave speed $c_0$ is given by the equality:

$$ c_0^2 = \min_{\mu>0} \frac{D_2 \mu^2 (\mu + \beta_2) e^{\mu \tau_2}}{(\mu + \sigma)(\mu + \beta_2) e^{\mu \tau_2} - b_2 a_2 u_0}. \quad (56) $$
It is similar to the previous one but with the parameters characterizing the second virus.

From (53) we can conclude that a positive solution exists if $a_1b_1 > a_2b_2$, that is, replication of the first virus is faster than of the second one. This condition appears also in Theorem 1. In the case of equality, the first virus is gradually replaced by the second one due to the mutation process (Figure 2). If the inequality is opposite, the first virus vanishes faster since its replication is slower. We will return to this question in the next section.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Numerical simulations of system (29)–(33) in the case where $a_1 = a_2 = b_1 = b_2$. Both virus spread in the form of waves, but the concentration of the first virus slowly decreases while of the second virus increases. The values of parameters: $a_1 = 0.01, a_2 = 0.01, b_1 = b_2 = 80,000; D_1 = D_2 = 10^{-3}, \beta_1 = \beta_2 = 0.1, \sigma_1 = \sigma_2 = 1, \tau = 2, u_0 = 1, i_{10} = i_{20} = 0, v_{1,0} = 4.5, v_{2,0} = 0, x_0 = 0.1$.</figcaption></figure>

### 3.3. Multiple Mutations and Coinfections

Generic model of virus mutation and coinfections can be written as follows:

$$ \frac{\partial U}{\partial t} = -U \sum_{i=1}^{n} a_i V_i, \qquad (57) $$

$$ \frac{\partial I_i}{\partial t} = a_i U V_i - \beta_i I_i, \quad i = 1, \dots, n, \qquad (58) $$

$$ \frac{\partial V_i}{\partial t} = D_i \frac{\partial^2 V_i}{\partial x^2} + \sum_{j=1}^{n} b_{ij} I_j - \sigma_i V_i, \quad i = 1, \dots, n. \qquad (59) $$

The coefficient $b_{ij}$ (positive or negative) show that viruses can promote or downregulate replication of other viruses in the case of coinfection. As before, under certain conditions on parameters, we can set $I_i = p_i I$, $V_i = q_i V$, and reduce this system of equations to the one-virus model. In the case of negative coefficients $b_{ij}$ the positiveness of solution in Equation (59) does not directly follow from the equation, and it should be controlled.

## 4. Different Stages of Infection Progression

### 4.1. Single Virus

In the previous sections, we have studied reaction-diffusion waves of infection spreading in cell culture. The solution of the initial boundary value problem approaches the wave after the initial transient period. However, this initial period cannot be neglected in the investigation of viral infection since it determines further infection progression. It is particularly important from the point of view of interaction with the immune response, though we do not consider it in this work (see [13–15]).

In order to describe analytically the initial stage of infection development, instead of system (11)–(13), we consider the system

$$ \frac{\partial U}{\partial t} = -aU V, \qquad (60) $$
$$ \frac{\partial I}{\partial t} = aU_0V - \beta I, \qquad (61) $$

$$ \frac{\partial V}{\partial t} = D \frac{\partial^2 V}{\partial x^2} + bI_{\tau} - \sigma V, \qquad (62) $$

where we replaced $U$ by $U_0$ in Equation (61). This means that we neglect the depletion of uninfected cells in virus replication. This is a linear system which can be further simplified if we integrate the last two equations. In order to compare the results with the numerical simulations, we consider this system in a bounded interval $0 < x < L$ with the boundary conditions

$$ x = 0, L: \frac{\partial V}{\partial x} = 0. \qquad (63) $$

Using the notation $J(t) = \int_0^L I(x,t)dx$, $W(t) = \int_0^L V(x,t)dx$, we obtain the delay differential system of equations

$$ \frac{dJ}{dt} = aU_0W - \beta J, \qquad (64) $$

$$ \frac{dW}{dt} = bJ_{\tau} - \sigma W \qquad (65) $$

with the initial condition

$$ J(t) = 0, \quad -\tau \le t \le 0, \quad W(0) = W_0 = \int_0^L V(x,0)dx. \qquad (66) $$

This system can be solved analytically. We restrict ourselves here to the two time intervals, $t \in [0, \tau)$:

$$ J(t) = \frac{aU_0W_0}{-\beta + \sigma}(e^{-\beta t} - e^{-\sigma t}), \quad W(t) = W_0e^{-\sigma t}, \qquad (67) $$

and $t \in [\tau, 2\tau)$:

$$ W(t) = W_0e^{-\sigma t} + \frac{abU_0W_0}{(\beta - \sigma)^2} [e^{-\beta(t-\tau)}(1 + e^{h(t)}(-1 + h(t)))], $$

$$ J(t) = \frac{aU_0W_0}{(\beta - \sigma)}(e^{-\sigma t} - e^{-\beta t}) + \frac{a^2bU_0^2W_0}{(\beta - \sigma)^3}(e^{-\sigma(t-\tau)}(-2 + h(t) + e^{-h(t)}(2 + h(t)))) $$

where $h(t) = (\beta - \sigma)(t - \tau)$, used for the comparison with the numerical simulations.

Figure 3 (left) shows the virus concentration distribution in consecutive moments of time with time difference between them equal 1 h. Let us note that since the length of the space interval is large and the diffusion coefficient is small, the descending part of the function $V(x,t)$ looks steep in the graphical representation. It is smoother in Figure 4 and even more in Figure 1. Beyond the graphical representation, important question is about numerical accuracy. This descending part of solution contains about 130 discretization points, and the accuracy of the solution is verified.

The right panel of Figure 3 illustrates how the total viral load (the integral of virus concentration with respect to the whole space interval) changes in time. The total viral load is determined in direct numerical simulations of the space-dependent problem (11)–(13) and by two approximate analytical methods. The first one is based on solution of linear ODE (64), (65). It is applicable for a relatively small time ($t < 4$). The second analytical approximation uses Theorem 1 and the formula for the wave speed that allows us to determine the viral load in the propagating wave for $t$ large enough.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** **Left:** virus distribution in consecutive moments of time (every 1 h). **Right:** total viral load $W(t) = \int_0^L V(x,t)dx$ (in the logarithmic scale) found as solution of system (64), (65) (green), of system (11)–(13) (blue) and by the formula in Theorem 1 (red). The values of parameters: $a = 0.01, b = 8 \times 10^4, D = 10^{-3}, \beta = 0.1, \sigma = 1, \tau = 2$.</figcaption></figure>

## 4.2. Virus Mutation and Competition

In Section 3, we have established conditions of the coexistence of two viruses in cell culture. If these conditions are not satisfied, and the new variant replicates faster, it will eliminate the first one (Figure 4).

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Distributions of $V_1(x,t)$ (**left**) and $V_2(x,t)$ (**right**) in numerical simulations of system (29)–(33). In the beginning of simulations, there is only the first virus variant. After some time, the second variant completely eliminate the first one. The values of parameters: $a_1 = 0.01, a_2 = 0.1, b_1 = b_2 = 8 \times 10^4, D_1 = D_2 = 10^{-3}, \beta_1 = \beta_2 = 0.1, \sigma_1 = \sigma_2 = 1, \tau = 2, u_0 = 1, i_{10} = i_{20} = 0, v_{1,0} = 4.5, v_{2,0} = 0, x_0 = 0.1$.</figcaption></figure>

Analysis of the total viral load (Figure 5) shows that there are four stages of infection progression in this case instead of the three stages observed previously. After the infection decay and explosive growth, the distribution of the first virus reaches its bell-shape form and propagates along the culture. At the same time, the total viral load of the second virus variant exponentially grows. When it reaches its maximum, the first viral load begins its exponential decay.

Denote by $T$ the transition time defined as the moment of time at which $V_2(t)$ reaches 90% of its maximal value. Transition time depends on the mutation rate and on the replication rates determined by the parameters $a_i, b_i, i = 1, 2$ (Figure 5, right). As it can be expected, $T$ is larger for small $\epsilon$, and it decreases if the replication rate of the second virus increases. The dependence of $T$ on $\epsilon$ is well approximated by the function $p/\epsilon^q$ with some fitted values of $p$ and $q$. Let us note that for these values of $q$, the dependence of the transition time on $\epsilon$ is weak. Extrapolating the curves using the analytical approximation,
we can see that decreasing the mutation rate by 10⁶ increases the transition time only about twice.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Left: $V_1(t)$ (red) and $V_2(t)$ (green) in the logarithmic scale for $b_2 = 8 \times 10^4$; $V_2(t)$ (blue) for $b_2 = 8 \times 10^5$. Right: transition time $T$ as a function of mutation rate $\epsilon$ for $b_2 = 8 \times 10^4$ (upper curve, crosses) and $b_2 = 8 \times 10^5$ (lower curve, crosses). Smooth lines represent analytical approximation with the functions $p/\epsilon^\theta$, $p = 25$, $q = 0.057$ (upper curve), $p = 17$, $q = 0.04$ (lower curve). The values of parameters: $a_1 = 0.01, a_2 = 0.1, b_1 = 8 \times 10^4$; $D_1 = D_2 = 10^{-3}$, $\beta_1 = \beta_2 = 0.1, \sigma_1 = \sigma_2 = 1, \tau = 2, u_0 = 1, i_{10} = i_{20} = 0, v_{1,0} = 4.5, v_{2,0} = 0, x_0 = 0.1$.</figcaption></figure>

### 4.3. Values of Parameters

The values of parameters for the single variant were determined in [10] by fitting the experimental data on the variation of viral load in time. Parameter estimation begins with the first stage of infection development during which virus is not yet reproduced by the infected cells, and its concentration exponentially decays. Duration of this stage and the decay rate allow us to uniquely determine parameters $\tau$ and $\sigma$ by comparison with the total viral load in the experimental data. The next stage is characterized by explosive virus production with the maximum of the total viral load determined by the parameter $b$. At the same time, parameter $a$ is chosen for the best fit of the growth curve. Finally, the third stage of infection development corresponds to its propagation in cell culture as a reaction-diffusion wave. The total viral load remains constant for $\beta = 0$, and it has a slow growth for a positive $\beta$, also determined by comparison with the total viral load in the experiment. Virus diffusion coefficient $D$ is estimated from the literature. Let us recall the dimensions of the main parameters: $[L] \sim \text{cm}$, $[time] \sim \text{hours}$, $[\tau] \sim \text{hours}$, $[\sigma] \sim 1/\text{hour}$, $[D] \sim \text{cm}^2/\text{hour}$.

The replication rate for the second virus (constants $a_2, b_2$ versus $a_1, b_1$) is taken in the same range or larger than that for the first virus. In order to evaluate the mutation rate $\epsilon$, let us note that the rate of mutations of RNA viruses is from $10^{-6}$ to $10^{-4}$ per base per generation. Taking into account that only from 0.001 to 0.01 viruses represent plaque forming units, we estimate the probability of a given mutation during one virus replication of the order of $10^{-9}$ to $10^{-6}$. The virus replication rate with respect to a single cell in the model is given by the constant $b_1 = 8 \times 10^4$. Multiplied by the probability determined above, we obtain the value of $\epsilon$ in the range from $10^{-4}$ to $10^{-1}$. Further, we can reasonably assume that viruses of the same size have close diffusion coefficients. In particular, this is the case of two variants of the same virus. Therefore, we take the same values of diffusion coefficient for both viruses.

## 5. Discussion

Multiplicity of virus assays represent conventional tests for the assessment of virus virulence and other properties [3,5]. In a previous work [10], we have shown that infection spreading in cell culture can be described by a reaction-diffusion system with time delay.
Three stages of infection progression were identified: virus decay, explosive growth, wave propagation. In this work, we apply this modeling approach to study virus mutation during infection progression. This question has a particular importance in the context of the emergence of new virus variants in COVID-19 pandemic. In a simplified setting, without taking into account the influence of the immune response, we can model this process as infection progression in cell culture and analyze the conditions of the emergence of new virus strains.

We introduce in this work the notion of virus replication number $R_v$ (similar to the basic reproduction number in epidemiological models) and show that infection spreads in cell culture if this number calculated for the first virus is larger than 1. Furthermore, the total viral load can be expressed in terms of $R_v$.

If the first virus variant spreads in cell culture and it can give a new virus variant due to mutations, then the two virus can co-exist or the second virus can eliminate the first one. The condition of their coexistence is given by the inequality $a_1b_1 > a_2b_2$. Taking into account that the coefficients $a_i$ characterize the rate of cell infection and the coefficients $b_i$ the rate of virus replication in the infected cells, their product $a_i b_i$ describes the total rate of virus production in cell culture. Hence, the two viruses coexist if the first one is produced faster than the second one. Otherwise, the second virus eliminates the first one.

The transition time from the first virus to the second virus depends on the mutation rate and on the replication rates. It is an important characterization of new variants since it shows whether they have enough time to develop in patients during the disease duration. The estimates of the transition time show that its dependence on the mutation rate is weak. Therefore, once new variant appears, it will rapidly replace the original one if it replicates faster, and the infected individual will become a career of the new variant.

This work is limited to the investigation of the emergence of new virus variants in the experimental setting of cell culture. The influence of the immune response on the infection mutation and spreading in the human organism will be considered in further works. However, some preliminary conclusions can be already done on the basis of the results of this work, in particular, about the action of vaccines on disease progression. In the case of SARS-CoV-2, the vaccine acts on virus spike preventing its penetration into the host cells. In the framework of the model, vaccine action can be described by decreasing the coefficient $a$. If it becomes small enough, such that the virus replication number $R_v$ is less than 1, then infection will not develop. If $R_v$ decreases but remains greater than 1, then infection will develop. However, infection spreading speed decreases. For the patients it signifies that disease symptoms are less severe.

**Author Contributions:** Conceptualization, V.V.; software, L.A.M.; formal analysis, B.K.; investigation, L.A.M., B.K. and V.V.; writing V.V.; All authors have read and agreed to the published version of the manuscript.

**Funding:** The second author was supported by the National Science Centre grant 2016/21/B/ST1/03071. The last author has been supported by the RUDN University Strategic Academic Leadership Program.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest

## References

1. Czerkies, M.; Korwek, Z.; Prus, W.; Kochanczyk, M.; Jaruszewicz-Blonska, J.; Tudelska, K.; Blonski, S.; Kimmel, M.; Brasier, A.R.; Lipniacki, T. Cell fate in antiviral response arises in the crosstalk of IRF, NF-kB and JAK/STAT pathways. *Nat. Commun.* **2018**, *9*, 493. [CrossRef] [PubMed]

2. Dykes, C.; Wang, J.; Jin, X.; Planelles, V.; An, D.S.; Tallo, A.; Huang, Y.; Wu, H.; Demeter, L.M. Evaluation of a multiple-cycle, recombinant virus, growth competition assay that uses flow cytometry to measure replication efficiency of human immunodeficiency virus type 1 in cell culture. *J. Clin. Microbiol.* **2006**, *44*, 1930–1943. [CrossRef] [PubMed]
3. Jegouic, S.; Joffret, M.L.; Blanchard, C.; Riquet, F.B.; Perret, C.; Pelletier, I.; Colbere-Garapin, F.; Rakoto-Andrianarivelo, M.; Delpeyroux, F. Recombination between polioviruses and co-circulating coxsackie A viruses: Role in the emergence of pathogenic vaccine-derived polioviruses. *PLoS Pathog.* **2009**, *5*, e1000412. [CrossRef] [PubMed]

4. Lindsay, S.M.; Timm, A.; Yin, J. A quantitative comet infection assay for influenza virus. *J. Virol. Methods* **2012**, *179*, 351–358. [CrossRef] [PubMed]

5. Sims, A.C.; Baric, R.S.; Yount, B.; Burkett, S.E.; Collins, P.L.; Pickles, R.J. Severe acute respiratory syndrome coronavirus infection of human ciliated airway epithelia: Role of ciliated cells in viral spread in the conducting airways of the lungs. *J. Virol.* **2005**, *79*, 15511–15524. [CrossRef] [PubMed]

6. Yin, J.; McCaskill, J.S. Replication of viruses in a growing plaque: A reaction-diffusion model. *Biophys. J.* **1992**, *61*, 1540–1549. [CrossRef]

7. Holder, B.P.; Simon, P.; Liao, L.E.; Abed, Y.; Bouhy, X.; Beauchemin, C.A.; Boivin, G. Assessing the in vitro fitness of an Oseltamivir-resistant seasonal A/H1N1 influenza strain using a mathematical model. *PLoS ONE* **2011**, *6*, e14767. [CrossRef] [PubMed]

8. Akpinar, F.; Inankur, B.; Yin, J. Spatial-temporal patterns of viral amplification and interference initiated by a single infected cell. *J. Virol.* **2016**, *90*, 7552–7566. [CrossRef] [PubMed]

9. Rodriguez-Brenes, I.A.; Hofacre, A.; Fan, H.; Wodarz, D. Complex dynamics of virus spread from low infection multiplicities: Implications for the spread of oncolytic viruses. *PLoS Comput. Biol.* **2017**, *13*, e1005241. [CrossRef] [PubMed]

10. Ait Mahiout, L.; Bessonov, N.; Kazmierczak, B.; Sadaka, G.; Volpert, V. Infection spreading in cell culture as a reaction-diffusion wave. *EDP Sci.* **2022**, submitted.

11. Hecht, F. New development in FreeFEM++. *J. Numer. Math.* **2012**, *20*, 251–256. [CrossRef]

12. Kolmogorov, A.N.; Petrovskii, I.G.; Piskunov, N.S. Etude de l'équation de la chaleur avec croissance de la quantité de matière et son application à un problème biologique. *Bull. Moskov. Gos. Univ. Mat. Mekh.* **1937**, *1*, 1–25.

13. Bessonov, N.; Bocharov, G.; Meyerhans, A.; Popov, V.; Volpert, V. Existence and Dynamics of Strains in a Nonlocal Reaction-Diffusion Model of Viral Evolution. *SIAM J. Appl. Math.* **2021**, *81*, 107–128. [CrossRef]

14. Bocharov, G.; Meyerhans, A.; Bessonov, N.; Trofimchuk, S.; Volpert, V. Modelling the dynamics of virus infection and immune response in space and time. *Int. J. Parallel Emergent Distrib. Syst.* **2019**, *34*, 341–355. [CrossRef]

15. Bocharov, G.; Meyerhans, A.; Bessonov, N.; Trofimchuk, S.; Volpert, V. Interplay between reaction and diffusion processes in governing the dynamics of virus infections. *J. Theor. Biol.* **2018**, *457*, 221–236. [CrossRef] [PubMed]