Space and Genotype-Dependent Virus Distribution during Infection Progression

: 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. frequent genotype, and it rapidly decays as the genotype moves away. This virus quasi-species progresses in the tissue (variable x ).


Introduction
Viral infections are often accompanied by rapid virus mutations, leading to the emergence of virus quasispecies [1][2][3][4][5] and their evolution [6][7][8]. Virus quasispecies are characterized by the presence of different but related genotypes [9][10][11][12], facilitating their adaptation to the environment, facilitating more efficient infection of host cells and tissues [13][14][15], escaping the immune system [16][17][18][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][23][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 Mathematics 2022, 10, 96 2 of 17 as a continuous variable, and of time t. Taking into account random virus motion, its mutations, multiplication and death, we obtain the following equation considered in the space interval 0 ≤ x ≤ L x and genotype interval |y| ≤ L y with the no-flux boundary conditions 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) = 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 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 σ(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 ≥ 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: for y ∈ 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 ± (y) = U(±∞, y) at infinity. Here, U ± (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. 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 [39] or genotype distribution [34]. Systems of two equations for the virus density and for the concentration of immune cells were considered in [40]. In this work we will consider the interaction of spatial virus distribution with its genotype characterization. Figure 1: A snapshot of solution u(x, y, t) of equation (1.1) without immune response (f (u) = 0) and with a piece-wise constant genotype-dependent mortality σ(y). The cross section in the y-direction (genetic variable) shows a stationary distribution describing virus quasispecies. 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).
Let us consider the stationary 1D equation in the y-direction: for y ∈ 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 solutions of this equation decaying at infinity corresponds to a virus quasi-species concentrated around some most frequent genotypes and decaying as the genotype goes away from this value [41].
3 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 σ(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).
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 = 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).

Wave Existence and Stability
We begin with a particular case b = 0, for which Equation (1) writes We consider it in the whole R 2 , and we look for its solution in the form: assuming that the initial condition u(x, y, 0) satisfies a similar representation: After the separation of variables where λ is a constant. Hence, Let λ 0 < 0 be the principle eigenvalue of the operator and φ 0 the corresponding eigenfunction. Then, φ 0 (y) > 0 for all y ∈ R. Without loss of generality, we suppose that I(φ) = 1. We set λ = −λ 0 > 0. Then, φ 0 (y) satisfies this equation. We can now write Equation (6) as follows: 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 σ(y) is continuous, and it has the limits σ ± = lim y→±∞ σ(y) at infinity, σ ± < 0. Let λ 0 be the principal eigenvalue of the operator L 0 φ = D y φ − σ(y)φ such that max(σ + , σ − ) < λ 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)φ 0 (y), where v(x, t) satisfies of Equation (8), where λ = −λ 0 , v(x, 0) = v 0 (x), φ 0 (y) is the eigenfunction of the operator L 0 corresponding to the principal eigenvalue, and φ 0 (y) > 0 for all y ∈ 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 λ 0 > max(σ + , σ − ) [38]. The assumption that the function σ(y) has limits at infinity allows for an explicit determination of the essential spectrum, S e = {λ ∈ R, λ ≤ max(σ + , σ − ). 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.
Such solutions exist for all c ≥ 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 reactiondiffusion equation in the monostable case [39]. According to the construction, travelling wave solutions of Equation (4) have the form U c (ξ, y) = w c (ξ)φ(y), where w c (ξ), c ≥ c 0 is a monotonically decreasing solution of the problem An example of an analytical solution is constructed in the Appendix A.
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.

1D Problem Depending on the Eigenvalue
We consider the one-dimensional equation where For b = 0, this equation is obtained from the 2D problem by separation of variables. For b = 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.

Properties of the Principal Eigenvalue
Consider the eigenvalue problem on the whole real axis. In order to get a more complete information about the dependence of the principal eigenvalue λ 0 on parameters, we consider a piece-wise constant function σ(y): where y 0 and σ 0 are some positive numbers. We will assume that σ 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 µ = −λ 0 /D y , σ 1 = σ 0 /D y . Then, µ satisfies the equation Consider the minimal positive solution of this equation located at the first growing branch of tangent for 0 < √ µ < π/(2y 0 ). Since the right-hand side of Equation (12) is a decreasing function of µ, such solution exists for any y 0 and σ 1 . The corresponding eigenfunction is given by the following expression: Assuming that I(φ 0 ) = 1, we find k 1 : Consider the limiting cases. If y 0 → 0, then it follows from Equation (12) that µ → σ 1 . Hence λ 0 → −σ 0 , k 1 → 0. If y 0 → ∞, then it follows from Equation (12) that µ → 0 and tan √ µy 0 → ∞. From the second equality for k 1 , we conclude that k 1 → 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 σ(y).

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 1 v + a 2 )e −a 3 v 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.
where λ = −λ 0 > 0, and λ 0 is the principal eigenvalue of problem (11). The function F λ (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).
Equation (10) has a solution v(x, t) = V(x − ct), that is, travelling wave with a speed c. For the function F λ (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(−∞) = v 1 and V(∞) = 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 λ (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 ≥ c 0 , where the minimal speed c 0 can be found by the formula [39]. These waves are stable in a properly chosen space [40].
Let us recall that λ is determined by the principal eigenvalue λ 0 = −λ of problem (11). It can be found from Equation (12) (µ = −λ 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 λ (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 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 λ (v) considered here, c 0 > c 1 . In this case, there are two consecutive waves propagating one after another. The monostable 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 (−∞) = v * , v(∞) = 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) ≡ 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.

Dependence on λ
We discussed above behavior of solutions of Equation (11) for a fixed value of λ for which the function F λ (v) has three positive zeros. Let us now analyze how this behavior depends on the value of λ. Let us recall that λ depends on y 0 and decreases from σ 0 to 0 as y 0 increases.
There are four critical values of λ that determine the behaviour of solutions ( Table 1). The first one, λ 1 = 1, is such that F λ (v) < 0 for all v > 0 if 1 < λ < σ 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 λ < λ 1 (= 1), then the function F λ (v) has a positive zero v 1 . There exist monostable [v 0 , v 1 ]-waves, as described above.
If λ becomes less than the next critical value λ 2 , then two more positive zeros v 2 and v 3 of the function F λ (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 λ < λ 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 λ 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 λ = λ 4 and disappear ( Figure 2, upper curve). We return to the monostable case with the [v 0 , v 3 ]-waves.

No.
λ Stationary Points Waves Speed Parameter a 1 characterizes the strength of immune response. Varying a 1 , we get the types of function F λ (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 λ (v) < 0 for all positive v and locally stable if F λ (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

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

Influence of the Genotype Distribution
Virus density distribution with respect to the genotype is determined by the function σ(y). We will vary the value of y 0 characterizing the width of this distribution, keeping σ 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 has a positive solution U 1 (y), and there is a wave propagating in the x-direction with the limits u(−∞, y) = U 1 (y), u(∞, 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).  (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:

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 1 c , then both solutions converge to the trivial solution. In both cases, we find y 1 c = 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 1 c . 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 2 c < y 0 < y 3 c , along with the monostable wave, there exists a bistable wave with a negative speed. The speed becomes positive for y > y 3 c (lines 3 and 4 in Table 1). For the chosen values of parameters, y 3 c ≈ 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).

Parameter Dependence
For the values of parameters considered above (a 1 = 15, a 2 = 0, a 3 = 6) and λ = 0, the function F λ (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 λ 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 ≈ 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 (λ ≈ 0.02).
Modelling and simulations presented above were carried out for a 2 = 0. Function F λ (v) becomes qualitatively different if a 2 > 1. The stationary solution v 0 = 0 of the equation dv/dt = F λ (v) becomes stable, and the function F λ (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 λ = 0.046 (a 1 = 13, a 2 = 1.1, a 3 = 6, σ 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.

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 λ 0 of problem (11). Let us note that proportional decrease of D y and y 0 does not change λ 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 ≤ x ≤ l x and |y| ≤ 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 √ 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 approximate 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 − λ(y 0 ), where λ(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 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 − λ(y 0 ) for the 1D problem. The values of parameters: D x = D y = 0.005, b = 1, k = 1, f (u) = (a 1 u + a 2 )e −a 3 u , a 1 = 15, and a 3 = 6.
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 D x (1 − λ(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 6. Dependence of virulence (wave speed) on immunity (a 2 ). The minimal wave speed is found by the formula c 0 = 2 D x k(1 − λ(y 0 ) − a 2 ) for the 1D problem. The upper curve corresponds to λ = 0 (limit of large y 0 ), the middle curve to λ = 0.17 (y 0 = 0.1), and the lower curve to λ = 0.41 (y 0 = 0.2). The values of parameters: 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 λ (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 (−∞) = v * , v 0 (∞) = 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 then some critical value, then there are no positive zeros of the function F λ (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 λ (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 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.

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.

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 : It is a discrete version of the diffusion equation. Diffusion coefficient D y = µ(∆y) 2 depends on the probability of mutations µ and on the characteristic space interval ∆y. The half-length y 0 of the admissible interval can be represented as y 0 = k∆y, where k is the number of mutations from the most frequent genotype to the end of viability interval.
The ratio D y /y 2 0 = µ/k 2 , independent of ∆y, determines the principle eigenvalue of the operator L 0 . Therefore, solutions are characterized by the value µ/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.

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.

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.

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 where c is the speed of the wave, and σ(y) = 0 , |y| ≤ 1 σ 0 , |y| > 1 .
This theorem affirms that for each solution λ 1 of Equation (A7) there is a set of solutions of Equation (A1). If σ 0 < π 2 /4, then λ 1 is unique. Increasing σ 0 leads to the increase of the number of solutions λ 1 .