# Space and Genotype-Dependent Virus Distribution during Infection Progression

Nicholas Bessonov, Gennady Bocharov, Vitaly Volpert

► To cite this version:

Nicholas Bessonov, Gennady Bocharov, Vitaly Volpert. Space and Genotype-Dependent Virus Distribution during Infection Progression. Mathematics, 2022, 10 (1), pp.96. 10.3390/math10010096. hal-04134233

HAL Id: hal-04134233

https://hal.science/hal-04134233v1

Submitted on 11 Mar 2025

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

# Space and Genotype-Dependent Virus Distribution during Infection Progression

Nicholas Bessonov 1,2, Gennady Bocharov 2,3,4* and Vitaly Volpert 2,5,6,7,*

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

check for
updates

**Citation:** Bessonov, N.; Bocharov, G.; Volpert, V. Space and Genotype-Dependent Virus Distribution during Infection Progression. Mathematics **2022**, *10*, 96. https://doi.org/10.3390/math10010096

Academic Editor: Andrea Scapellato

Received: 29 November 2021

Accepted: 23 December 2021

Published: 28 December 2021

**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:** © 2021 by the authors. License 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/).

**Abstract:** The paper is devoted to a nonlocal reaction-diffusion equation describing the development of viral infection in tissue, taking into account virus distribution in the space of genotypes, the antiviral immune response, and natural genotype-dependent virus death. It is shown that infection propagates as a reaction-diffusion wave. In some particular cases, the 2D problem can be reduced to a 1D problem by separation of variables, allowing for proof of wave existence and stability. In general, this reduction provides an approximation of the 2D problem by a 1D problem. The analysis of the reduced problem allows us to determine how viral load and virulence depend on genotype distribution, the strength of the immune response, and the level of immunity.

**Keywords:** virus density distribution; genotype; nonlocal interaction; wave propagation; infection progression

## 1. Introduction

Viral infections are often accompanied by rapid virus mutations, leading to the emergence of virus quasispecies [1–5] and their evolution [6–8]. Virus quasispecies are characterized by the presence of different but related genotypes [9–12], facilitating their adaptation to the environment, facilitating more efficient infection of host cells and tissues [13–15], escaping the immune system [16–19] and resisting antiviral drugs [20,21]. The genetic variation of viruses (e.g., HIV type 1) can generate mutants that escape from CTL recognition [22–24]. Mutations that are not detrimental and still lead to escape are limited in number and mathematically can be described by a confined domain in the genotype space. For example, the diversity of the HIV quasispecies observed during HIV infection increases at rate of about 1% per year (envelope gene) [25]. For a non-treated infection lasting on average 10 years, the overall diversity will be about 10%. Immune system can also adapt to virus evolution [26], resulting in the increase of target regions and in the shift in the immunodominance [27,28].

Mathematical models of virus evolution can be based on ODE systems [27,29] or stochastic models [30,31], and they take into account the interactions in viral populations and the disease pathogenesis [13,32]. Let us also mention some models that take into account spatial virus distribution with diffusion and chemotaxis [33].

In this work, we will consider virus density distribution $u(x, y, t)$ in some tissue of the organisms as a function of the space variable $x$, of the genotype characteristics $y$ considered

<footnote n="1">1 Institute of Problems of Mechanical Engineering, Russian Academy of Sciences, 199178 Saint Petersburg, Russia; bnm@ipme.ru</footnote>

<footnote n="2">2 Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences, 119333 Moscow, Russia; g.bocharov@inm.ras.ru</footnote>

<footnote n="3">3 Moscow Center of Fundamental and Applied Mathematics at INM RAS, 119333 Moscow, Russia</footnote>

<footnote n="4">4 Institute of Computer Science and Mathematical Modelling, Sechenov First Moscow State Medical University, 119991 Moscow, Russia</footnote>

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

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

<footnote n="7">7 Peoples' Friendship University of Russia (RUDN University), 6 Miklukho-Maklaya Str., 117198 Moscow, Russia</footnote>

<footnote>* Correspondence: volpert@math.univ-lyon1.fr</footnote>
as a continuous variable, and of time $t$. Taking into account random virus motion, its mutations, multiplication and death, we obtain the following equation

$$ \frac{\partial u}{\partial t} = D_x \frac{\partial^2 u}{\partial x^2} + D_y \frac{\partial^2 u}{\partial y^2} + ku(1 - bu) - f(I(u))u - \sigma(y)u \quad (1) $$

considered in the space interval $0 \le x \le L_x$ and genotype interval $|y| \le L_y$ with the no-flux boundary conditions

$$ x = 0, L_x: \frac{\partial u}{\partial x} = 0, \quad y = \pm L_y: \frac{\partial u}{\partial y} = 0. \quad (2) $$

Diffusion coefficients $D_x$ and $D_y$, and virus multiplication rate $k$, are positive constants; inverse carrying capacity $b$ is non-negative; and $I(u) = \int_{-L_y}^{L_y} u(x, y, t)dy$. In the mathematical analysis of this problem, we will also consider it in the whole two-dimensional space $\mathbb{R}^2$, assuming that both variables $x$ and $y$ can adopt all real values. This mathematical approximation is applicable if the solution of this equation evolves sufficiently far from the boundary of the domain.

The first term in the right-hand side of Equation (1) describes virus displacement in the tissue, and the second term characterizes its random mutations. Virus multiplication rate is determined by the logistic function with a normalized carrying capacity. Virus death occurs due to its elimination by the immune cells or due to its natural genotype-dependent mortality with the rate $\sigma(y)u$.

The concentration of immune cells $C$ depends on the virus density since adaptive immune response is stimulated by the antigen. Cytotoxic T-lymphocytes migrate towards the infected tissue attracted by the inflammatory cytokines, and their local concentration is proportional to the local concentration of the infected cells. Hence, in the simplest formulation we set $C = f(u)$, where $f(u)$ is a non-negative function determined for $u \ge 0$. Two cases can be considered: this function can grow with saturation or it can grow for $u$ sufficiently small and decrease for $u$ sufficiently large due to death of the immune cells taking place at high virus densities. Assuming that immune response is independent of virus genotype, we get its dependence on the total virus density, $f = f(I(u))$.

In our previous works, we have studied 1D problems with either spatial virus distribution [34] or genotype distribution [35]. Systems of two equations for the virus density and for the concentration of immune cells were considered in [36]. In this work, we will consider the interaction of spatial virus distribution with its genotype characterization.

Let us consider the stationary 1D equation in the $y$-direction:

$$ D_y \frac{d^2 u}{dy^2} + ku(1 - bu - f(I(u)) - \sigma(y)u = 0 \quad (3) $$

for $y \in \mathbb{R}$. This equation describes virus distribution with respect to its genotype, taking into account its reproduction and its elimination by the immune response and due to the genotype-dependent mortality. A positive solution to this equation decaying at infinity corresponds to a virus quasi-species concentrated around some of the most frequent genotypes and decaying as the genotype goes away from this value [37].

The goal of this work is to study infection spread in the tissue or organ, taking into account the virus distribution with respect to its genotype. Therefore, we introduce the second variable $x$ corresponding to the space coordinate in the tissue. Clearly, real biological tissues should be considered in 3D. The 1D approximation is justified if the solution does not depend on the other two variables. We study in this work how the virus quasi-species spreads in the tissue. From the mathematical point of view, we study propagation of travelling waves of Equation (1), that is, solutions of the form $u(x, y, t) = U(x - ct, y)$ with some limits $U_{\pm}(y) = U(\pm\infty, y)$ at infinity. Here, $U_{\pm}(y)$ are some solutions of Equation (3).
Figure 1 shows a typical structure of travelling wave propagating in the x-direction and characterizing the quasi-species distribution with respect to the genotype variable y.

<figure><img src="image_4.png" /><figcaption>Figure 1. A snapshot of solution $u(x, y, t)$ of Equation (1) without immune response ($f(u) = 0$) and with a piece-wise constant genotype-dependent mortality $\sigma(y)$. The cross section in the y-direction (genetic variable) shows a stationary distribution describing virus quasi-species. It has a maximum for some most frequent genotype, and it rapidly decays as the genotype moves away. This virus quasi-species progresses in the tissue (variable x).</figcaption></figure>

We begin with a particular case where $b = 0$ and the 2D problem can be reduced to some 1D problem depending on the x-variable (Section 2). It includes as a parameter the principal eigenvalue of some eigenvalue problem in the y-direction. This reduction allows us to carry out a detailed study of the dynamics of solutions. If $b \neq 0$, then the 1D problem depending on parameter provides an approximation of the 2D problem (Sections 3 and 4). We use it to investigate the main features of infection progression and to determine how virus virulence and viral load depend on the width of genotype distribution, the strength of immune response, and the level of immunity (Section 5).

## 2. Wave Existence and Stability

We begin with a particular case $b = 0$, for which Equation (1) writes

$$ \frac{\partial u}{\partial t} = D_x \frac{\partial^2 u}{\partial x^2} + D_y \frac{\partial^2 u}{\partial y^2} + ku(1 - f(I_y(u))) - \sigma(y)u. \quad (4) $$

We consider it in the whole $\mathbb{R}^2$, and we look for its solution in the form:

$$ u(x, y, t) = v(x, t)\phi(y), $$

assuming that the initial condition $u(x, y, 0)$ satisfies a similar representation:

$$ u(x, y, 0) = v_0(x)\phi(y). \quad (5) $$

After the separation of variables

$$ \left( -\frac{\partial v}{\partial t} + D_x \frac{\partial^2 v}{\partial x^2} + kv(1 - f(vI(\phi))) \right) \phi = -v(D_y\phi'' + \sigma(y)\phi) $$
we obtain the equality

$$ \frac{-\frac{\partial v}{\partial t} + D_x \frac{\partial^2 v}{\partial x^2} + kv(1 - f(vI(\phi)))}{v} = \frac{-D_y \phi'' + \sigma(y)\phi}{\phi} = \lambda, $$

where $\lambda$ is a constant. Hence,

$$ \frac{\partial v}{\partial t} = D_x \frac{\partial^2 v}{\partial x^2} + kv(1 - f(vI(\phi))) - \lambda v, \quad (6) $$

$$ D_y \phi'' - \sigma(y) \phi + \lambda \phi = 0. \quad (7) $$

Let $\lambda_0 < 0$ be the principle eigenvalue of the operator

$$ L_0 \phi = D \phi'' - \sigma(y) \phi, $$

and $\phi_0$ the corresponding eigenfunction. Then, $\phi_0(y) > 0$ for all $y \in \mathbb{R}$. Without loss of generality, we suppose that $I(\phi) = 1$. We set $\lambda = -\lambda_0 > 0$. Then, $\phi_0(y)$ satisfies this equation. We can now write Equation (6) as follows:

$$ \frac{\partial v}{\partial t} = D_x \frac{\partial^2 v}{\partial x^2} + kv(1 - f(v)) - \lambda v. \quad (8) $$

Thus, we have reduced 2D nonlocal Equation (4) to 1D local Equation (8). This reduction allows us to formulate the results on wave existence and stability.

**Theorem 1.** Suppose that the function $\sigma(y)$ is continuous, and it has the limits $\sigma_{\pm} = \lim_{y \to \pm\infty} \sigma(y)$ at infinity, $\sigma_{\pm} < 0$. Let $\lambda_0$ be the principal eigenvalue of the operator $L_0\phi = D_y\phi'' - \sigma(y)\phi$ such that $\max(\sigma_+, \sigma_-) < \lambda_0 < 0$. Then, solution $u(x, y, t)$ of the Cauchy problem (4), (5) can be represented in the form $u(x, y, t) = v(x, t)\phi_0(y)$, where $v(x, t)$ satisfies of Equation (8), where $\lambda = -\lambda_0$, $v(x, 0) = v_0(x)$, $\phi_0(y)$ is the eigenfunction of the operator $L_0$ corresponding to the principal eigenvalue, and $\phi_0(y) > 0$ for all $y \in \mathbb{R}$.

The proof of the theorem is a consequence of the consideration above. We take into account that the principal eigenvalue of the operator $L_0$ is real and simple, and that the corresponding eigenfunction is positive if $\lambda_0 > \max(\sigma_+, \sigma_-)$ [38]. The assumption that the function $\sigma(y)$ has limits at infinity allows for an explicit determination of the essential spectrum, $S_e = \{\lambda \in \mathbb{R}, \lambda \le \max(\sigma_+, \sigma_-)\}$. If the limits do not exist, the essential spectrum cannot be explicitly determined but it can be estimated. The theorem remains applicable if the principal eigenvalue of the operator $L_0$ is greater than its upper bound. Let us also note that solution of problem (4), (5) is unique.

Reduction of Equation (4) to the 1D equation provided by Theorem 1 leads to the following result on the wave existence.

**Theorem 2.** Under conditions of the previous theorem, let the function $F_{\lambda}(u) = ku(1 - f(u)) - \lambda u$ have a single positive zero $u_1$, such that $F'_{\lambda}(u_1) < 0$. Then, Equation (4) has travelling wave solutions $u(x, y, t) = U(x - ct, y)$, monotonically decreasing as a function of $x$, with the limits

$$ U(-\infty, y) = u_1\phi_0(y), \quad U(\infty, y) = 0. $$

Such solutions exist for all $c \ge c_0$, where $c_0 > 0$ is a minimal speed, and they do not exist for $c < c_0$.

The proof of the theorem follows from the corresponding results for the scalar reaction-diffusion equation in the monostable case [39]. According to the construction, travelling wave solutions of Equation (4) have the form $U_c(\xi, y) = w_c(\xi)\phi(y)$, where $w_c(\xi)$, $c \ge c_0$ is a monotonically decreasing solution of the problem
$$D_x w'' + cw' + F_\lambda(w) = 0, \quad w(-\infty) = u_1, \quad w(\infty) = 0. \qquad (9)$$

An example of an analytical solution is constructed in the Appendix A.

**Theorem 3.** Suppose that the function $F_\lambda(u)$ is such that $F_\lambda(0) = F_\lambda(u_1) = 0$, $F_\lambda(u) > 0$ for $0 < u < u_1$, and $F'_\lambda(0) > 0$. If $v_0(x) = 0$ for all $x$ sufficiently large and $\lim_{x \to -\infty} v_0(x) > 0$, then solution $u(x, y, t)$ of Equation (4) with initial condition (6) converges to the wave $U_{c_0}$ with the minimal speed $c_0$.

As before, this theorem follows from the corresponding results for the 1D Equation [39] (see [40] and the references therein). Convergence here occurs in form and in speed. More general results and convergence to the waves with other speeds can be obtained.

### 3. 1D Problem Depending on the Eigenvalue

We consider the one-dimensional equation

$$\frac{\partial v}{\partial t} = D_x \frac{\partial^2 v}{\partial x^2} + kv(K_\lambda - bv - f(v)), \qquad (10)$$

where $K_\lambda = 1 - \lambda/k$, and $-\lambda = \lambda_0$ is the principal eigenvalue of the operator $L_0v = D_yv'' - \sigma(y)v$. For $b = 0$, this equation is obtained from the 2D problem by separation of variables. For $b \neq 0$, this reduction does not work, and we will use this equation as an approximation of the 2D problem. We set here $b = 1$ and $k = 1$.

#### 3.1. Properties of the Principal Eigenvalue

Consider the eigenvalue problem

$$D_y \phi'' - \sigma(y) \phi = \lambda \phi, \quad \phi(\pm \infty) = 0 \qquad (11)$$

on the whole real axis. In order to get a more complete information about the dependence of the principal eigenvalue $\lambda_0$ on parameters, we consider a piece-wise constant function $\sigma(y)$:

$$\sigma(y) = \begin{cases} 0, & |y| \le y_0 \\ \sigma_0, & |y| > y_0 \end{cases},$$

where $y_0$ and $\sigma_0$ are some positive numbers. We will assume that $\sigma_0 > 1$ (for $k = 1$). Hence, the essential spectrum of problem (3) belongs to the left-half plane.

This eigenvalue problem can be solved analytically. Let us note that the eigenvalue with the maximal real part of this problem is negative. Set $\mu = -\lambda_0/D_y$, $\sigma_1 = \sigma_0/D_y$. Then, $\mu$ satisfies the equation

$$\tan(\sqrt{\mu} y_0) = \sqrt{\frac{\sigma_1}{\mu} - 1}, \quad \mu < \sigma_1. \qquad (12)$$

Consider the minimal positive solution of this equation located at the first growing branch of tangent for $0 < \sqrt{\mu} < \pi/(2y_0)$. Since the right-hand side of Equation (12) is a decreasing function of $\mu$, such solution exists for any $y_0$ and $\sigma_1$. The corresponding eigenfunction is given by the following expression:

$$\phi_0(y) = \begin{cases} k_1 \cos(\sqrt{\mu} y), & |y| \le y_0 \\ k_1 \cos(\sqrt{\mu} y_0) e^{-\sqrt{\sigma_1 - \mu}|y - y_0|}, & |y| > y_0 \end{cases}.$$

Assuming that $I(\phi_0) = 1$, we find $k_1$:

$$k_1 = \frac{\mu \sqrt{\sigma_1 - \mu}}{2\sigma_1 \cos(\sqrt{\mu} y_0)} = \frac{\sqrt{\mu}(\sigma_1 - \mu)}{\sigma_1 \sin(\sqrt{\mu} y_0)}.$$
Consider the limiting cases. If $y_0 \to 0$, then it follows from Equation (12) that $\mu \to \sigma_1$. Hence $\lambda_0 \to -\sigma_0$, $k_1 \to 0$. If $y_0 \to \infty$, then it follows from Equation (12) that $\mu \to 0$ and $\tan(\sqrt{\mu}y_0) \to \infty$. From the second equality for $k_1$, we conclude that $k_1 \to 0$. Hence, the maximum of the eigenfunction increases for small $y_0$ and decreases for large $y_0$. Similar properties hold for a smoothed function $\sigma(y)$.

## 3.2. Dynamics of 1D Waves

Dynamics of solutions of Equation (10) are studied for a wide class of functions $f(v)$ (see [40] and the references therein). We will consider for certainty the function $f(v) = (a_1v + a_2)e^{-a_3v}$ qualitatively describing the properties of immune response as a function of viral load. Here, $a_1, a_2, a_3$ are some non-negative numbers. Set $F_\lambda(v) = v(1 - \lambda - v - f(v))$, where $\lambda = -\lambda_0 > 0$, and $\lambda_0$ is the principal eigenvalue of problem (11). The function $F_\lambda(v)$ vanishes at $v_0 = 0$, and it has up to three positive zeros $v_1, v_2, v_3$ denoted in the growing order (Figure 2).

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Function $F_\lambda(v) = v(1 - \lambda - v - f(v))$ for $\lambda = 0$ (upper curve), $\lambda = 0.1$ (middle curve), and $\lambda = 0.2$ (lower curve), $f(v) = (a_1v + a_2)e^{-a_3v}$, $a_1 = 13, a_2 = 0, a_3 = 6$.</figcaption></figure>

Equation (10) has a solution $v(x, t) = V(x - ct)$, that is, travelling wave with a speed c. For the function $F_\lambda(v)$ shown in Figure 2 (middle curve), there are three positive zeros and two different types of waves. The $[v_0, v_1]$-waves are the waves with the limits $V(-\infty) = v_1$ and $V(\infty) = v_0(= 0)$ at infinity. This is so-called monostable case since the point $v_0$ is unstable with respect to the equation $dv/dt = F_\lambda(v)$, and the point $v_1$ is stable. Such waves can be monotone as functions of x and non-monotone. The latter are unstable, and we will not consider them here. The monotone waves exist for all $c \ge c_0$, where the minimal speed $c_0$ can be found by the formula $c_0 = 2\sqrt{D_xF'_\lambda(0)} = 2\sqrt{D_x(1 - \lambda)}$ since $F(v) \le F'(0)v$ [39]. These waves are stable in a properly chosen space [40].

Let us recall that $\lambda$ is determined by the principal eigenvalue $\lambda_0 = -\lambda$ of problem (11). It can be found from Equation (12) ($\mu = -\lambda_0/D_y$). Therefore, we can determine the minimal speed of the monostable wave. Its dependence on $y_0$ is shown in Figure 3 (dashed curve).

The second wave is the $[v_1, v_3]$-wave, which corresponds to the bistable case since both points $v_1$ and $v_3$ are stable as stationary points of equation $dv/dt = F_\lambda(v)$. Such wave exists for a unique value of speed $c_1$; it is monotone as a function of x, and it is globally asymptotically stable. Its speed cannot be found analytically except for the case $c_1 = 0$, which is reached if and only if the integral $J(v_1, v_3) = \int_{v_1}^{v_3} F_\lambda(v)dv$ equals 0. The case of zero speed separates the cases where infection spreads in the tissue (positive speed) and where it does not spread (negative speed). We will discuss this question below.

Behavior of solutions of Equation (11) depends on the relation between the speeds $c_0$ and $c_1$, and on the initial condition. For the function $F_\lambda(v)$ considered here, $c_0 > c_1$. In this case, there are two consecutive waves propagating one after another. The monostable $[v_0, v_1]$-wave is followed by the $[v_1, v_2]$-wave.
Denote the initial condition for Equation (11) by $v_0(x)$ (Cauchy problem). For simplicity of presentation, suppose that it is a monotonically decreasing function with the limits $v_0(-\infty) = v_*, v(\infty) = 0$ at infinity. Then, for $v_* < v_2$, the solution $v(x,t)$ converges to the monostable $[v_0, v_1]$-wave. We will restrict ourselves here to the convergence to the wave with the minimal speed that occurs, in particular, in the case where $v_0(x) \equiv 0$ for all $x$ that are sufficiently large. This case is most appropriate for applications and numerical simulations. If $v_* > v_2$, then the solution converges to the system of two waves, monostable and bistable, propagating one after another with different speeds.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Wave speed for the 1D Equation (10) (dashed line) and for the 2D Equation (1) (solid lines). Curves 1 and 3 correspond to the monostable and bistable waves, respectively, for the same values of parameters, and curve 2 shows the transition from the bistable wave to the monostable wave. The values of parameters are $D_x = D_y = 0.005$, $k = 1$, $\sigma_0 = 1.1$, $f(v) = a_1ve^{-a_3v}$, $a_3 = 6$, $a_1 = 15$ (curves 1, 3, and dashed curve), and $a_1 = 13$ (curve 2).</figcaption></figure>

### 3.3. Dependence on $\lambda$

We discussed above behavior of solutions of Equation (11) for a fixed value of $\lambda$ for which the function $F_\lambda(v)$ has three positive zeros. Let us now analyze how this behavior depends on the value of $\lambda$. Let us recall that $\lambda$ depends on $y_0$ and decreases from $\sigma_0$ to 0 as $y_0$ increases.

There are four critical values of $\lambda$ that determine the behaviour of solutions (Table 1). The first one, $\lambda_1 = 1$, is such that $F_\lambda(v) < 0$ for all $v > 0$ if $1 < \lambda < \sigma_0$. In this case, the waves do not exist. For any $v_0(x)$ (initial condition), the solution $v(x,t)$ converges to 0. If $\lambda < \lambda_1(= 1)$, then the function $F_\lambda(v)$ has a positive zero $v_1$. There exist monostable $[v_0, v_1]$-waves, as described above.

If $\lambda$ becomes less than the next critical value $\lambda_2$, then two more positive zeros $v_2$ and $v_3$ of the function $F_\lambda(v)$ appear. In this case, in addition to the monostable $[v_0, v_1]$-waves, there is a bistable $[v_1, v_3]$-wave with a negative speed. It becomes positive for $\lambda < \lambda_3$. This critical value is determined by the equality $J(v_1, v_3) = 0$. Let us note that the minimal speed increases with the decrease of $\lambda$ for the monostable wave, similar to the single speed of the bistable wave. If $0 < c_1 < c_0$, there is a system of two waves. As we discussed above, solution can converge either to the monostable wave or to the system of two waves depending on the initial condition. Zeros $v_1$ and $v_2$ merge for $\lambda = \lambda_4$ and disappear (Figure 2, upper curve). We return to the monostable case with the $[v_0, v_3]$-waves.
<table><caption>**Table 1.** Four critical values $\lambda_1, \lambda_2, \lambda_3, \lambda_4$, determine the intervals of $\lambda$ with different stationary points, waves, and their speeds.</caption><thead><tr><th>No.</th><th>λ</th><th>Stationary Points</th><th>Waves</th><th>Speed</th></tr></thead><tbody><tr><td>1</td><td>λ<sub>1</sub> &lt; λ &lt; σ<sub>0</sub></td><td>v<sub>0</sub> = 0</td><td>-</td><td>-</td></tr><tr><td>2</td><td>λ<sub>2</sub> &lt; λ &lt; λ<sub>1</sub></td><td>v<sub>0</sub>, v<sub>1</sub></td><td>monostable [v<sub>0</sub>, v<sub>1</sub>]-waves</td><td>c ≥ c<sub>0</sub> = 2√D<sub>x</sub>(1 - λ)</td></tr><tr><td>3</td><td>λ<sub>3</sub> &lt; λ &lt; λ<sub>2</sub></td><td>v<sub>0</sub>, v<sub>1</sub>, v<sub>2</sub>, v<sub>3</sub></td><td>monostable, bistable</td><td>c ≥ c<sub>0</sub>, c<sub>1</sub> &lt; 0</td></tr><tr><td>4</td><td>λ<sub>4</sub> &lt; λ &lt; λ<sub>3</sub></td><td>v<sub>0</sub>, v<sub>1</sub>, v<sub>2</sub>, v<sub>3</sub></td><td>monostable, bistable</td><td>c ≥ c<sub>0</sub>, 0 &lt; c<sub>1</sub> &lt; c<sub>0</sub></td></tr><tr><td>5</td><td>0 &lt; λ &lt; λ<sub>4</sub></td><td>v<sub>0</sub>, v<sub>3</sub></td><td>monostable [v<sub>0</sub>, v<sub>3</sub>]-waves</td><td>c ≥ c<sub>0</sub></td></tr></tbody></table>

### 3.4. Dependence on $f(v)$

Parameter $a_1$ characterizes the strength of immune response. Varying $a_1$, we get the types of function $F_\lambda(v)$ similar to those considered above. Parameter $a_2$ determines the presence of immunity from the previous infection or vaccine. If $a_2 > 1$, then $v_0$ is a stable stationary point. It is globally stable (for $v > 0$) if $F_\lambda(v) < 0$ for all positive $v$ and locally stable if $F_\lambda(v)$ changes sign. In the latter case, there is a bistable $[v_0, v_2]$-wave with a positive or negative depending on the sign of the integral $\int_{v_0}^{v_2} F_\lambda(v)dv$.

## 4. 2D Wave Dynamics

In this section, we study dynamics of solutions of Equation (1) with $b \neq 0$. We will show that it correlates with the dynamics of solutions of 1D equation considered in the previous section.

### 4.1. Influence of the Genotype Distribution

Virus density distribution with respect to the genotype is determined by the function $\sigma(y)$. We will vary the value of $y_0$ characterizing the width of this distribution, keeping $\sigma_0$ constant. For sufficiently small values of $y_0$, the trivial solution $u = 0$ of Equation (1) is stable and wave propagation is not observed. For larger values of $y_0$, the trivial solution loses its stability. Equation

$$D_y \frac{d^2 u}{dy^2} + ku(1 - u - f(I(u))) - \sigma(y)u = 0 \quad (13)$$

has a positive solution $U_1(y)$, and there is a wave propagating in the x-direction with the limits $u(-\infty, y) = U_1(y)$, $u(\infty, y) = U_0(= 0)$ at infinity (Figure 4, right). These limits are solutions of Equation (13). If we increase $y_0$, the speed of this wave also increases (curves 1 in Figure 3).

For $y_0 > 0.42$, another solution is observed for the same values of parameters. It consists of two consecutive waves (Figure 4, left). A small amplitude wave is the same as in Figure 4, right (the scaling in the figures is different). It is followed by a high amplitude wave with a smaller speed. The presence of this second wave indicates the existence of another positive stationary solution $U_3(y)$ of Equation (13), which is presumably stable. Therefore, we can expect that they are separated by an unstable solution $U_2(y)$ that is not observed in the simulations. Thus, there exists a monstable $[U_0, U_1]$-wave, and a bistable $[U_1, U_3]$-wave. If the initial condition is small enough, the single small-amplitude wave propagates. If the initial condition is sufficiently large, then we observe the system of two waves. The speed of the bistable wave increases with $y_0$ (curve 2 in Figure 3).
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Snapshots of solutions of Equation (1) with two different initial conditions and the same values of parameters. Two waves with different speeds propagate one after another for a sufficiently large initial condition (**left**). A faster monostable wave with a small amplitude is followed by a slower bistable wave with a large amplitude. If the initial condition is small enough, then only the fast monostable wave with a small amplitude is observed (**right**). The values of parameters: $D_x = D_y = 0.005, b = 1, k = 1, y_0 = 0.5, f(u) = (a_1u + a_2)e^{-a_3u}, a_1 = 15, a_2 = 0, a_3 = 6$.</figcaption></figure>

## 4.2. Comparison with 1D Problem

We will now compare the 2D problem with the 1D problem, depending on the eigenvalue for the same and fixed values of parameters and varying $y_0$. If $y_0$ is less than the first critical value, $0 < y_0 < y_c^1$, then both solutions converge to the trivial solution. In both cases, we find $y_c^1 = 0.02$. Let us recall that we determine the principal eigenvalue of problem (11) as a function of $y_0$. These values of $y_0$ correspond to the first line in Table 1.

The monostable wave is observed for both problems for $y_0 > y_c^1$. The speeds of such waves are shown in Figure 3 (curve 1 and dashed curve). The wave speeds have a similar dependence on $y_0$. We recall that this is the minimal speed of the 1D wave. In the case of 2D problem, we observe only one wave in numerical simulations. Therefore, strictly speaking, we cannot affirm that there are also waves with larger speeds. We can conjecture their existence by analogy with the 1D equation. This case corresponds to line 2 in Table 1.

In the next interval of values of $y_0$, $y_c^2 < y_0 < y_c^3$, along with the monostable wave, there exists a bistable wave with a negative speed. The speed becomes positive for $y > y_c^3$ (lines 3 and 4 in Table 1). For the chosen values of parameters, $y_c^3 \approx 0.42$ for the 2D problem and 0.38 for the 1D problem. The latter is found analytically from the condition $J(v_1, v_3) = 0$ (Section 3). The speed of the bistable wave increases as a function of $y_0$ (Figure 3, curve 3).

## 4.3. Parameter Dependence

For the values of parameters considered above ($a_1 = 15, a_2 = 0, a_3 = 6$) and $\lambda = 0$, the function $F_\lambda(v)$ has three positive zeros, $v_1, v_2$, and $v_3$. If we increase $a_1$, then zeros $v_2$ and $v_3$ approach each other, merge, and disappear. In this case, there is the only positive zero $v_1$ and the monostable $[v_0, v_1]$-waves. This case is similar to the case of large $\lambda$ considered above.

If we decrease $a_1$, then zeros $v_1$ and $v_2$ merge and disappear. In this case, there is the only positive zero $v_3$ and the monostable $[v_0, v_3]$-waves. This case was not considered above. Set $a_0 = 13$ and consider dynamics of 2D waves for different values of $y_0$. The difference in comparison with the previous case ($a_1 = 15$) manifests itself in the behavior of the bistable wave for $y_0 > 0.5$. Its speed rapidly grows and approaches the minimal speed of the monostable wave (Figure 3, curve 2). The bistable $[v_1, v_3]$-wave and the monostable $[v_0, v_1]$-wave disappear for $y_0 \approx 0.7$, resulting in the emergence of the monostable $[v_0, v_3]$-wave. In the 1D problem, this transition occurs for the same value of $y_0$ ($\lambda \approx 0.02$).
Modelling and simulations presented above were carried out for $a_2 = 0$. Function $F_\lambda(v)$ becomes qualitatively different if $a_2 > 1$. The stationary solution $v_0 = 0$ of the equation $dv/dt = F_\lambda(v)$ becomes stable, and the function $F_\lambda(v)$ may not have positive zeros or have two of them (the degenerate case of merged zeros is not considered). In the first case, the trivial solution of the 1D problem is globally asymptotically stable. In the second case, this problem has a bistable [$v_0, v_2$]-wave. The sign of its speed can be positive or negative depending on the other parameters. The transition between them is determined by the equality $J(v_0, v_2) = 0$. Numerical simulations of the 2D problem also reveal the existence of a large amplitude [$U_0, U_2$]-wave where $U_0 = 0, U_2(y) > 0$ is a solution of Equation (13). Its speed changes sign for $y_0 = 0.45$, for which $\lambda = 0.046$ ($a_1 = 13, a_2 = 1.1, a_3 = 6, \sigma_0 = 1.1, D_x = D_y = 0.005$, and $k = 1$), and the speed of the 1D is also close to zero.

Thus, dynamics of the 1D waves depending on the eigenvalue and of the 2D waves are similar to each other. They show the existence of the monostable and bistable waves, and of a monostable-bistable wave. Their parameter dependence is also analogous.

**5. The Properties of Infection Progression**

We will study in this section how infection spread in the tissue is influenced by the genotype distribution, the initial viral load, the strength of immune response, and the presence of immunity. We will characterize infection progression and immune response in terms of the parameters of the model.

* Genotype distribution. Virus genotype distribution depends on parameter $y_0$ characterizing the interval where the virus reproduction rate exceeds its natural mortality rate (without immune response), and on the mutation rate that determines the diffusion coefficient $D_y$. These two parameters determine the principal eigenvalue $\lambda_0$ of problem (11). Let us note that proportional decrease of $D_y$ and $y_0$ does not change $\lambda_0$, so that the results presented above are appropriate in a wide range of mutation rates.

* Initial viral load. The initial viral load corresponds to the initial condition for Equations (1) or (10). In 2D simulations, we set $u_0(x,y) = u_0^*$ for $0 \le x \le l_x$ and $|y| \le l_y$ with some positive numbers $l_x$ and $l_y$, and $u_0(x,y) = 0$ otherwise. In this context, the initial viral load $V_0$ is the integral of the initial condition, $V_0 = 2l_x l_y u_0^*$.

* Strength of immune response. Adaptive immune response proceeds by clonal expansion of lymphocytes due to the interaction with antigens-presenting cells (macrophages, dendritic cells). For small viral loads, increase of the level of pathogens in the organism intensifies the immune response. Some viral infections can affect the functioning of lymphocytes by downregulating their proliferation and increasing their death (e.g., HIV, LCMV [41] but not coronavirus). Therefore, the function $f(u)$ increases for small $u$ and decreases for large $u$. Stronger immune response corresponds to a larger function $f(u)$. In modelling, we characterize the strength of immune response by the value of parameter $a_1$.

* Immunity. Vaccination and previous infections can lead to the appearance of antibodies and memory cells responding to a new antigen. This response can be attenuated by the reduced affinity to the antigen of the immunity mediators (antibodies, T cell receptors). Immunity slows down infection progression and accelerates clonal expansion of immune cells. We model the presence of immunity by the coefficient $a_2$. If it is positive, then immune response starts from some positive value under the introduction of antigen, $f(0) > 0$. Infection-free equilibrium $v_0 = 0$ becomes stable for sufficiently narrow genotype distribution $y_0$, and infection is eliminated unless the initial viral load is sufficiently large to cause its persistent progression.

* Viral load. The level of infection in the tissue determines the severity of symptoms and the intensity of infection transmission to other individuals. Viral load is determined by all factors presented in the previous paragraphs. In modelling, it corresponds to the virus density distribution after the wave propagation.
* Virus virulence. Virus virulence is characterized by its spreading rate. In modelling, it is determined by the coefficient $k$ in Equation (1). In the absence of immune response, it is related to the minimal speed of the monostable wave $c_0 = 2\sqrt{D_x k}$. Multiplicity of infection tests in cell culture [42] characterize the virulence of infection by the speed of viral plaque growth, that is, by the wave speed. We will also characterize the virulence of infection by the wave speed in the presence of immune response.

We will determine how viral load and virulence depend on the strength of immune response, the level of immunity, and the width of the genotype distribution. Since the 2D problem can be approximated with the 1D problem depending on eigenvalue, we will use the latter for the investigation of infection progression. It will allow us to obtain a relatively simple analytical description of its main features.

If the disease-free equilibrium $v_0 = 0$ is stable, then infection does not develop for sufficiently small initial viral loads and solution uniformly converges to 0. If it is unstable, then a positive disease equilibrium appears and there is a wave of infection propagation. Figure 5 shows the stability region depending on the width of genotype distribution (parameter $y_0$) and on the level of immunity (coefficient $a_2 = f(0)$). If the latter is large enough, then the disease free equilibrium is stable. In the 1D problem, the stability boundary is determined by the function $a_2 = 1 - \lambda(y_0)$, where $\lambda(y_0)$ is the principal eigenvalue of problem (11) (up to sign). The stability boundary for the 1D problem is compared with the result of numerical simulation of the 2D problem. Let us note that the disease-free equilibrium can be globally stable or stable only to small perturbations. In the second case, there are two positive stationary points and a bistable wave describing disease progression for a sufficiently large initial viral load. We will return to this question below.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** The critical value of immunity (coefficient $a_2$) is shown as a function of the width of the genotype distribution ($y_0$). If the value of immunity exceeds the critical value, then infection is eliminated. Otherwise, it progresses. Solid line represents the results of 2D simulations, and dashed line is given by the formula $a_2 = 1 - \lambda(y_0)$ for the 1D problem. The values of parameters: $D_x = D_y = 0.005$, $b = 1$, $k = 1$, $f(u) = (a_1u + a_2)e^{-a_3u}$, $a_1 = 15$, and $a_3 = 6$.</figcaption></figure>

Next, we study how the virulence of infection depends on the level of immunity. We determine the virulence as the minimal speed of the monostable wave, $c_0 = 2\sqrt{D_x(1 - \lambda(y_0) - f(0))}$. For a fixed width of genotype distribution, it becomes a function of $a_2 = f(0)$, that is, a function of immune memory protection in the absence of pathogen. As it can be expected, the virulence of infection decreases with the increase of
immunity level (Figure 6). Infection does not develop if immunity level exceeds the critical value where the curve crosses the x-axis.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Dependence of virulence (wave speed) on immunity ($a_2$). The minimal wave speed is found by the formula $c_0 = 2\sqrt{D_x k(1 - \lambda(y_0) - a_2)}$ for the 1D problem. The upper curve corresponds to $\lambda = 0$ (limit of large $y_0$), the middle curve to $\lambda = 0.17$ ($y_0 = 0.1$), and the lower curve to $\lambda = 0.41$ ($y_0 = 0.2$). The values of parameters: $D_x = D_y = 0.005$, $b = 1$, $k = 1$, $f(v) = (a_1v + a_2)e^{-a_3v}$, $a_1 = 15$, and $a_3 = 6$.</figcaption></figure>

Viral load is an independent characteristic with respect to virus virulence. We determine it as the maximal level of virus concentration during infection progression in the tissue. Hence, in the 1D problem, it corresponds to positive zeros of the function $F_\lambda(v)$. There are two essentially different cases presented in Figure 7. Curve 1 corresponds to the bistable case with two positive zeros $v_1$ and $v_2$; $v_1 < v_2$; and a bistable [$v_0, v_2$]-wave. The value $v_1$ is on the lower branch of this curve, and $v_2$ on the upper branch. Let us recall that the initial condition in the 1D problem considered on the whole axis is a function $v_0(x)$ with the limits $v_0(-\infty) = v_*$, $v_0(\infty) = 0$. Here, $v_*$ is the initial viral load. If it exceeds $v_1$, then the solution converges to the bistable wave with the maximal value (viral load) $v_2$. If $v_* < v_1$, then the solution uniformly converges to 0 and infection is eliminated. If the strength of immune response (parameter $a_1$) is larger than some critical value, then there are no positive zeros of the function $F_\lambda(v)$. The disease-free equilibrium $v_0$ is globally stable, and infection is eliminated for any initial viral load.

The second curve corresponds to the case with up to three positive zeros of the function $F_\lambda(v)$ and the existence of either monostable or monostable-bistable waves. The lower branch of the curve shows the value of $v_1$, the middle branch $v_2$, and the upper branch $v_3$. For a weak immune response, there is only one positive solution at the upper branch and the corresponding large amplitude monostable [$v_0, v_3$]-wave. For a strong immune response, there is the monostable [$v_0, v_1$]-wave with small amplitude. For the intermediate values, there are three positive zeros and the monostable-bistable system of waves, where the monostable [$v_0, v_1$]-wave is followed by the bistable [$v_1, v_3$]-wave with a smaller speed. In this case, if the initial viral load is small enough, then the final viral load remains small. If the initial viral load is sufficiently large, then the viral load proceeds in two consecutive steps, first a small one, then a large one. In terms of disease progression, a mild form of the disease is followed by a severe form.
<figure><img src="image_10.png" /><figcaption>Figure 7. Dependence of the viral load (wave amplitude) on the strength of immune response ($a_1$) for the 1D problem. The wave amplitude is found as a solution of the equation $F_\lambda(v) = 0$. The values of parameters: $D_x = D_y = 0.005, b = 1, k = 1, f(v) = (a_1v + a_2)e^{-a_3v}$, Curve 1: $a_2 = 0.7, a_3 = 6, \lambda = 0.41$ ($y_0 = 0.1$). Curve 2: $a_2 = 0, a_3 = 6, \lambda = 0.17$ ($y_0 = 0.2$).</figcaption></figure>

# 6. Discussion

Viral infection development in the tissues of the host organism depends on many factors determined by the pathogen and by the immune response to the infection. In this work, we study the influence of virus genotype distribution on spatial infection spreading.

## 6.1. Model

The model suggested in this work is a generalization of the 1D model introduced in [34] for the infection spreading in the tissue. This 1D model describes spatiotemporal virus distribution taking into account its random motion, reproduction, and elimination due to the immune response.

In order to study virus density distribution with respect to its genotype, we need to introduce its genotypical characteristic as a continuous variable. The simplest way to do it is to consider a sequence of consecutive mutations and the corresponding genotypes denoted by $y_i$. Transition from the genotype $y_i$ to the genotype $y_{i+1}$ occurs due to one of the mutations of this sequence. Let $u_i$ be the virus density corresponding to the genotype $y_i$. Assuming that all mutations are reversible and have the same rate, we obtain the following equations for $u_i$:

$$ \frac{du_i}{dt} = \mu(u_{i+1} - 2u_i + u_{i-1}). $$

It is a discrete version of the diffusion equation. Diffusion coefficient $D_y = \mu(\Delta y)^2$ depends on the probability of mutations $\mu$ and on the characteristic space interval $\Delta y$. The half-length $y_0$ of the admissible interval can be represented as $y_0 = k\Delta y$, where $k$ is the number of mutations from the most frequent genotype to the end of viability interval.

The ratio $D_y/y_0^2 = \mu/k^2$, independent of $\Delta y$, determines the principle eigenvalue of the operator $L_0$. Therefore, solutions are characterized by the value $\mu/k^2$ and not by the individual parameters $D_y$ and $y_0$.

The combination of the space variable $x$ and the genotype variable $y$ allows us to study the influence of genotype distribution on infection spreading. The method developed in this work is also applicable for the equation with time delay.

## 6.2. Method and Results

The method developed in this work implies the reduction of the 2D problem to a 1D problem in the $x$-direction and to some eigenvalue problem in the $y$-direction. In a particular case, this reduction is exact. It allows us to study the existence and stability of
travelling waves. The case with time delay, though not considered in this work, can be studied by the same method. Dynamics of 1D solutions with time delay is investigated in [34,43].

In the general case, the exact reduction of the 2D problem to the 1D problem is not possible. We consider the 1D problem depending on the eigenvalue as an approximation of the 2D problem. In some limiting cases, this approximation might be possible to justify by more rigourous mathematical considerations but this question is beyond the scope of this work. Numerical simulations of the 2D problem are in a good agreement with the analytical results for the 1D problem.

We use the 1D problem depending on the eigenvalue in order to study the progression of viral infection. Virus virulence (wave speed) and viral load (wave amplitude) are determined by the width of the genotype distribution ($y_0$), strength of immune response ($a_1$), level of immunity ($f(0)$), and the initial viral load (initial condition). Detailed characterization of infection progression is presented in the previous section. Let us mention here that the increase of the width of the genotype distribution increases virus virulence and viral load.

### 6.3. Limitations and Perspectives

We consider a simplified model of infection development where immune response is reduced to a function $f(u)$ describing the immune cell concentration depending on virus density. This model follows from a more detailed model due to a quasi-stationary approximation [36]. Time delay characterizing the duration of the clonal expansion of immune cells is not taken into account in order to simplify the presentation and to concentrate on other aspects of the problem. The dependence of the virus density distribution on the genotype variable is considered in the simplest approximation of reversible consecutive mutations. Finally, the developed model provides the tool for comprehensive analysis and evaluation of the link between the viral evolution and immune escape from the immunodominant CD8+ T cell responses [44]. More complex mutation patterns can be considered.

**Author Contributions:** Investigation, V.V.; Methodology, G.B.; Software, N.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research was funded by the Russian Science Foundation (Grant no. 18-11-00171). G.B. was partly supported by Moscow Center of Fundamental and Applied Mathematics (agreement with the Ministry of Education and Science of the Russian Federation No. 075-15-2019-1624).

**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.

## Appendix A. Analytical Solution for Nonlocal Equation

We will construct here the analytical solution of the equation

$$D \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) + c \frac{\partial u}{\partial x} + ku(1 - I_y(u)) - \sigma(y)u = 0, \quad (A1)$$

where c is the speed of the wave, and

$$\sigma(y) = \begin{cases} 0, & |y| \le 1 \\ \sigma_0, & |y| > 1 \end{cases}. \quad (A2)$$
In what follows we set $D = 1$. We look for the solution of Equation (A1) in the form: $u(x,y) = v(x)\phi(y)$, where $\phi(y) = \phi_1(y)$ for $|y| \le 1$ and $\phi(y) = \phi_2(y)$ for $|y| > 1$. Substituting these expressions into the equation, we obtain for $|y| \le 1$:

$$v''\phi + v\phi'' + cv'\phi + kv\phi(1 - vI(\phi)) = 0,$$

where $v' = dv/dx$, $\phi' = d\phi/dy$, $I_y(u) = vI(\phi)$. Separating the variables in the last equation, we have

$$\frac{v'' + cv' + kv(1 - vI)}{v} = -\frac{\phi''}{\phi} = \lambda_1.$$

The left-hand side of this equality depends only on x, the next expression depends only on y. Therefore, both of them equal some constant $\lambda_1$. Hence,

$$\left\{ \begin{array}{ll} v'' + cv' + kv(1 - vI) - \lambda_1 v = 0 & , \quad |y| \le 1. \\ \phi'' + \lambda_1 \phi = 0 & , \quad |y| > 1. \end{array} \right. \qquad (A3)$$

Similarly,

$$\left\{ \begin{array}{ll} v'' + cv' + kv(1 - vI) - \sigma_0 v - \lambda_2 v = 0 & , \quad |y| > 1. \\ \phi'' + \lambda_2 \phi = 0 & , \quad |y| < 1. \end{array} \right. \qquad (A4)$$

We set

$$\lambda_1 = \lambda_2 + \sigma_0, \qquad (A5)$$

such that $v(x)$ satisfies the same equation in (A3) and in (A4). Next, we suppose that $\lambda_1 > 0$ and $\lambda_2 < 0$. Then

$$\phi_1(y) = a_1 \cos(\sqrt{\lambda_1}y), \quad \phi_2(y) = a_2 e^{-\sqrt{\mu}y},$$

where $\mu = -\lambda_2$, $a_1$ and $a_2$ are arbitrary constants. From the conditions $\phi_1(1) = \phi_2(1)$, $\phi'_1(1) = \phi'_2(1)$ (similar conditions at $y = -1$),

$$\left\{ \begin{array}{ll} a_1 \cos(\sqrt{\lambda_1}) = a_2 e^{-\sqrt{\mu}} & , \quad y \le -1 \\ -a_1 \sqrt{\lambda_1} \sin(\sqrt{\lambda_1}) = -\sqrt{\mu} a_2 e^{-\sqrt{\mu}} & , \quad -1 \le y \le 1 \\ a_1 \sqrt{\lambda_1} \cos(\sqrt{\lambda_1}) e^{-\sqrt{\mu}y} & , \quad y \ge 1 \end{array} \right. . \qquad (A6)$$

Dividing the second equation in (A6) by the first one and taking into account (A5), we obtain the equation with respect to $\lambda_1$:

$$\tan \sqrt{\lambda_1} = \sqrt{\frac{\sigma_0 - \lambda_1}{\lambda_1}}, \quad 0 < \lambda_1 < \sigma_0. \qquad (A7)$$

This equation has at least one solution for any $\sigma_0 > 0$. It can have more than one solution for $\sigma_0$ sufficiently large. Thus,

$$\phi(y) = \begin{cases} ae^{\sqrt{\mu}} \cos(\sqrt{\lambda_1})e^{\sqrt{\mu}y} & , \quad y \le -1 \\ a \cos(\sqrt{\lambda_1}y) & , \quad -1 \le y \le 1 \\ ae^{\sqrt{\mu}} \cos(\sqrt{\lambda_1})e^{-\sqrt{\mu}y} & , \quad y \ge 1 \end{cases}, \qquad (A8)$$

where $a$ is an arbitrary positive constant, $\lambda_1$ is a solution of Equation (A7), $\mu = -\lambda_2 = \sigma_0 - \lambda_1$,

$$I(\phi) = 2a \left( \frac{1}{\sqrt{\mu}} \cos(\sqrt{\lambda_1}) + \frac{1}{\sqrt{\lambda_1}} \sin(\sqrt{\lambda_1}) \right). \qquad (A9)$$

In order to determine solutions of the first equation in (A5), consider the problem

$$w'' + cw' + kw(1 - w) - \lambda_1 w = 0, \qquad (A10)$$

$$w(-\infty) = (k - \lambda_1)/k, \quad w(\infty) = 0. \qquad (A11)$$
Equation (A10) is similar to the KPP equation. It has a monotonically decreasing solution if $\lambda_1 < k$ for all $c \ge c_0$, where $c_0 = 2\sqrt{k - \lambda_1}$.

From the first equation in (A5) we obtain $v = w/I$. Hence, $u(x,y) = w(x)\phi(y)/I(\phi)$. This solution does not depend on the arbitrary constant $a$. The following theorem is proved.

**Theorem A1.** Suppose that $\lambda_1$ is a solution of Equation (A7), $0 < \lambda_1 < k$. Then, Equation (A1) has a solution $u(x,y) = w(x)\phi(y)/I(\phi)$ for all $c \ge c_0 = 2\sqrt{k - \lambda_1}$, where $w(x)$ is a solution of problem (A10), (A11), $\phi(y)$ and $I(\phi)$ are determined by expressions (A8), (A9). Its second derivative with respect to $x$ is continuous for all $x \in \mathbb{R}$, its second derivative with respect to $y$ is continuous for all $y \in \mathbb{R}$ except for $y = \pm 1$.

This theorem affirms that for each solution $\lambda_1$ of Equation (A7) there is a set of solutions of Equation (A1). If $\sigma_0 < \pi^2/4$, then $\lambda_1$ is unique. Increasing $\sigma_0$ leads to the increase of the number of solutions $\lambda_1$.

## Referentes

1. Biebricher, C.K.; Eigen, M. What is a quasispecies? *Curr. Top. Microbiol. Immunol.* **2006**, *299*, 1–31.

2. Domingo, E.; Perales, C. Viral quasispecies. *PLoS Genet.* **2019**, *15*, e1008271. [CrossRef]

3. Goodenow, M.; Huet, T.; Saurin, W.; Kwok, S.; Sninsky, J.; Wain-Hobson, S. HIV-1 isolates are rapidly evolving quasispecies: Evidence for viral mixtures and preferred nucleotide substitutions. *J. Acquir. Immune Defic. Syndr.* **1989**, *2*, 344–352.

4. Holland, J.J.; De La Torre, J.C.; Steinhauer, D.A. RNA virus populations as quasispecies. *Curr. Top. Microbiol. Immunol.* **1992**, *176*, 1–20. [PubMed]

5. Meyerhans, A.; Cheynier, R.; Albert, J.; Seth, M.; Kwok, S.; Sninsky, J.; Morfeldt-Manson, L.; Asjo, B.; Wain-Hobson, S. Temporal fluctuations in HIV quasispecies in vivo are not reflected by sequential HIV isolations. *Cell* **1989**, *58*, 901–910. [CrossRef]

6. Poirier, E.Z.; Vignuzzi, M. Virus population dynamics during infection. *Curr. Opin. Virol.* **2017**, *23*, 82–87. [CrossRef]

7. Kautz, T.F.; Forrester, N.L. RNA Virus Fidelity Mutants: A Useful Tool for Evolutionary Biology or a Complex Challenge? *Viruses* **2018**, *10*, 600. [CrossRef]

8. Karamitros, T.; Papadopoulou, G.; Bousali, M.; Mexias, A.; Tsiodras, S.; Mentis, A. SARS-CoV-2 exhibits intra-host genomic plasticity and low-frequency polymorphic quasispecies. *J. Clin. Virol.* **2020**, *131*, 104585. [CrossRef]

9. Keele, B.H.; Giorgi, E.E.; Salazar-Gonzalez, J.F.; Decker, J.M.; Pham, K.T.; Salazar, M.G.; Sun, C.; Grayson, T.; Wang, S.; Li, H.; et al. Identification and characterization of transmitted and early founder virus envelopes in primary HIV-1 infection. *Proc. Natl. Acad. Sci. USA* **2008**, *105*, 7552–7557. [CrossRef]

10. Plikat, U.; Nieselt-Struwe, K.; Meyerhans, A. Genetic drift can dominate short-term human immunodeficiency virus type 1 nef quasispecies evolution in vivo. *J. Virol.* **1997**, *71*, 4233–4240. [CrossRef]

11. Mattenberger, F.; Vila-Nistal, M.; Geller, R. Increased RNA virus population diversity improves adaptability. *Sci. Rep.* **2021**, *11*, 6824. [CrossRef] [PubMed]

12. Jary, A.; Leducq, V.; Malet, I.; Marot, S.; Klement-Frutos, E.; Teyssou, E.; Soulié, C.; Abdi, B.; Wirden, M.; Pourcher, V.; et al. Evolution of viral quasispecies during SARS-CoV-2 infection. *Clin. Microbiol. Infect.* **2020**, *26*, 1560.e1–1560.e4. [CrossRef] [PubMed]

13. Vignuzzi, M.; Stone, J.K.; Arnold, J.J.; Cameron, C.E.; Andino, R. Quasispecies diversity determines pathogenesis through cooperative interactions in a viral population. *Nature* **2006**, *439*, 344–348. [CrossRef]

14. Karakasiliotis, I.; Lagopati, N.; Evangelou, K.; Gorgoulis, V.G. Cellular senescence as a source of SARS-CoV-2 quasispecies. *FEBS J.* **2021**, *Epub ahead of print*. [CrossRef]

15. Sun, F.; Wang, X.; Tan, S.; Dan, Y.; Lu, Y.; Zhang, J.; Xu, J.; Tan, Z.; Xiang, X.; Zhou, Y.; et al. SARS-CoV-2 Quasispecies Provides an Advantage Mutation Pool for the Epidemic Variants. *Microbiol. Spectr.* **2021**, *9*, e0026121. [CrossRef]

16. Collier, D.A.; Monit, C.; Gupta, R.K. The Impact of HIV-1 drug escape on the global treatment landscape. *Cell Host Microbe* **2019**, *26*, 48–60. [CrossRef]

17. Esposito, I.; Trinks, J.; Soriano, V. Hepatitis C virus resistance to the new direct-acting antivirals. *Expert Opin. Drug Metab. Toxicol.* **2016**, *12*, 1197–1209. [CrossRef] [PubMed]

18. Kiepiela, P.; Leslie, A.J.; Honeyborne, I.; Ramduth, D.; Thobakgale, C.; Chetty, S.; Rathnavalu, P.; Moore, C.; Pfaerott, K.J.; Hilton, L. Dominant inuence of HLA-B in mediating the potential co-evolution of HIV and HLA. *Nature* **2004**, *432*, 769–775. [CrossRef]

19. Illingworth, C.J.R.; Raghwani, J.; Serwadda, D.; Sewankambo, N.K.; Robb, M.L.; Eller, M.A.; Redd, A.R.; Quinn, T.C.; Lythgoe, K.A. A de novo approach to inferring within-host fitness effects during untreated HIV-1 infection. *PLoS Pathog.* **2020**, *16*, e1008171. [CrossRef] [PubMed]

20. Phillips, R.E.; Rowland-Jones, S.; Nixon, D.F.; Gotch, F.M.; Edwards, J.P.; Ogunlesi, A.O.; Elvin, J.G.; Rothbard, J.A.; Bangham, C.R.; Rizza, C.R.; et al. Human immunodeficiency virus genetic variation that can escape cytotoxic T cell recognition. *Nature* **1991**, *354*, 453–459. [CrossRef]
21. Larder, B.A.; Kemp, S.D. Multiple mutations in HIV-1 reverse transcriptase confer high-level resistance to zidovudine (AZT). *Science* **1989**, *246*, 1155–1158. [CrossRef]

22. McMichael, A.J.; Borrow, P.; Tomaras, G.D.; Goonetilleke, N.; Haynes, B.F. The immune response during acute HIV-1 infection: Clues for vaccine development. *Nat. Rev. Immunol.* **2010**, *10*, 11–23. [CrossRef]

23. Yang, Y.; Ganusov, V.V. Kinetics of HIV-Specific CTL Responses Plays a Minimal Role in Determining HIV Escape Dynamics. *Front. Immunol.* **2018**, *9*, 140. [CrossRef] [PubMed]

24. Conway, J.M.; Ribeiro, R.M. Modeling the immune response to HIV infection. *Curr. Opin. Syst. Biol.* **2018**, *12*, 61–69. [CrossRef]

25. Bocharov, G.; Chereshnev, V.; Gainova, I.; Bazhan, S.; Bachmetyev, B.; Argilaguet, J.; Martinez J.; Meyerhans, A. Human Immunodeficiency Virus Infection: From Biological Observations to Mechanistic Mathematical Modelling. *Math. Model. Nat. Phenom.* **2012**, *7*, 78–104. [CrossRef]

26. Haas, G.; Plikat, U.; Debré, P.; Lucchiari, M.; Katlama, C.; Dudoit, Y.; Bonduelle, O.; Bauer, M.; Ihlenfeldt, H.G.; Jung, G.; et al. Dynamics of viral variants in HIV-1 Nef and specific cytotoxic T lymphocytes in vivo. *J. Immunol.* **1996**, *157*, 4212–4221.

27. Ganusov, V.V.; Goonetilleke, N.; Liu, M.K.; Ferrari, G.; Shaw, G.M.; McMichael, A.J.; Borrow, P.; Korber, B.T.; Perelson, A.S. Fitness costs and diversity of the cytotoxic T lymphocyte (CTL) response determine the rate of CTL escape during acute and chronic phases of HIV infection. *J. Virol.* **2011**, *85*, 10518–10528. [CrossRef]

28. Turnbull, E.L.; Wong, M.; Wang, S.; Wei, X.; Jones, N.A.; Conrod, K.E.; Aldam, D.; Turner, J.; Pellegrino, P.; Keele, B.F.; et al. Kinetics of expansion of epitope-specific T cell responses during primary HIV-1 infection. *J. Immunol.* **2009**, *182*, 7131–7145. [CrossRef]

29. Ganusov, V.V.; Neher, R.A.; Perelson, A.S. Mathematical modeling of escape of HIV from cytotoxic T lymphocyte responses. *J. Stat. Mech.* **2013**, P01010. [CrossRef] [PubMed]

30. Bocharov, G.A.; Ford, N.J.; Edwards, J.; Breinig, T.; Wain-Hobson, S.; Meyerhans S.A. A genetic algorithm approach to simulating human immunodeficiency virus evolution reveals the strong impact of multiply infected cells and recombination. *J. Gen. Virol.* **2005**, *86*, 3109–3118. [CrossRef] [PubMed]

31. Merleau, N.S.C.; Pénisson, S.; Gerrish, P.J.; Elena, S.F.; Smerlak, M. Why are viral genomes so fragile? The bottleneck hypothesis. *PLoS Comput. Biol.* **2021**, *17*, e1009128. [CrossRef]

32. Swanstrom, R.; Coffin, J. HIV-1 pathogenesis: The virus. *Cold Spring Harb. Perspect. Med.* **2012**, *2*, a007443. [CrossRef]

33. Bellomo, N.; Painter, K.J.; Tao, Y.; Winkler, M. Occurrence vs. Absence of Taxis-Driven Instabilities in a May-Nowak Model for Virus Infection. *SIAM J. Appl. Math.* **2019**, *79*, 1990–2010. [CrossRef]

34. Bocharov, G.; Meyerhans, A.; Bessonov, N.; Trofimchuk, S.; Volpert, V. Spatiotemporal dynamics of virus infection spreading in tissues. *PLoS ONE* **2016**, *11*, e0168576. [CrossRef] [PubMed]

35. Bessonov, N.; Bocharov, G.; Meyerhans, A.; Popov, V.; Volpert, V. Nonlocal reaction-diffusion model of viral evolution: Emergence of virus strains. *Mathematics* **2020**, *8*, 117. [CrossRef]

36. 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]

37. Bessonov, N.; Bocharov, G.A.; Leon, C.; Popov, V.; Volpert, V. Genotype-dependent virus distribution and competition of virus strains. *Math. Mech. Complex Syst.* **2020**, *8*, 101–126. [CrossRef]

38. Volpert, A.; Volpert, V. Spectrum of elliptic operators and stability of travelling waves. *Asymptot. Anal.* **2000**, *23*, 111–134.

39. 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.

40. Volpert, A.; Volpert, V.; Volpert, V. *Traveling Wave Solutions of Parabolic Systems*; Translation of Mathematical Monographs; American Mathematical Society: Providence, RI, USA, 1994; Volume 140.

41. Bocharov, G.; Volpert, V.; Ludewig, B.; Meyerhans, A. *Mathematical Immunology of Virus Infections*; Springer: Berlin/Heidelberg, Germany, 2018.

42. Holder, B.P.; Simon, P.; Liao, L.E.; Abed, Y.; Bouhy, X.; Beauchemin, C.A.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]

43. Bessonov, N.; Bocharov, G.; Touaoula, T.M.; Trofimchuk, S.; Volpert, V. Delay reaction-diffusion equation for infection dynamics. *Discret. Contin. Dyn. Syst. B* **2019**, *24*, 2073. [CrossRef]

44. Henn, M.R.; Boutwell, C.L.; Charlebois, P.; Lennon, N.J.; Power, K.A.; Macalalad, A.R.; Berlin, A.M.; Malboeuf, C.M.; Ryan, E.M.; Gnerre, S.; et al. Whole Genome Deep Sequencing of HIV-1 Reveals the Impact of Early Minor Variants Upon Immune Recognition During Acute Infection. *PLoS Pathog.* **2012**, *8*, e1002529. [CrossRef] [PubMed]