Viruses competition in the genotype space

This paper is devoted to the study of persistence and evolution of two viruses taking into account virus mutation, reproduction, and genotype dependent mortality, either natural or determined by an antiviral treatment. The model describes the virus density distribution u(x; t) for the first virus and v(y; t) for the second one as functions of genotypes x and y considered as continuous variables and of time t. The model consists of a system of reaction-diffusion equations with integral terms characterizing virus competition for host cells. The analysis of the model shows the conditions of the existence of virus strains.


Introduction
Viruses are in a constant evolution owing to variation of their genetic structure as the result of the interaction of the replication, recombination or mutation [1]. Further, these changes in the virus genetic structure contribute to the appearance of new variations thereof. The variability of the viruses is considered to be one of the key factors in the pathogonesis of the respective infectious disease, besides the high genetic variety and quickness of mutation of some viruses makes the immune response ineffective and it can lead to the emergence of the resistance to the antiviral therapy [2]. A better understanding of infection development would allow us the improvement techniques of viral treatments and their effectiveness.
When a virus enters the host organism, it begins to contaminate uninfected cells starting its replication. The concentration of virus in the host organism can be affected by its elimination, either by the immune response, virus natural death or some antiviral treatment. In [3], a model that describes the evolution of virus density depending on the genotype is introduced and the conditions for the existence of virus strains are determined.
In this study we propose a mathematical model which considers the aforementioned processes that affects the concentration of virus in the host organism without taking into account the immune response, in the case of the existence of two viruses u, v in the host. We consider the system of equations: Theses equations describe the evolution of virus densities depending on the genotypes x and y respectively, considered as continuous variables and on time. Here D 1 , D 2 are the coefficients of diffusion and parameters α 1,2 , β 1,2 , γ 1,2 are positive constants. The first terms in the righthand side of these equations characterize virus mutation, the second terms its reproduction, and the last terms the virus death. We now detail each of these terms of the first equation.
They are similar in for the second equation.
• Assuming that there is a sequence of reversible mutations with consecutive genotypes x i , we can write the equation for the density u i of virus with genotype x i : where µ is the frequency of mutations. This equation represents a discretization of the diffusion equation with the diffusion coefficient proportional to µ.
• Virus multiplication term is proportional to the virus density u and to the quantity of uninfected host cells (1 − β 1 I(u) − γ 1 I(v)) which is valid for the acute stage of infection (virus multiplication phase). Here 1 is a dimensionless total number of cells, β 1 I(u) and γ 1 I(v) are the number of infected cells by virus u and v, respectively, which are proportional to the total viruses quantity: • The last term in the right-hand side of equation (1) describes virus natural death or its elimination by some antiviral treatment. Let us note that the death rate can depend on virus genotype x or y.
We consider the virus strain as density distribution concentrated around some genotype value. It is a non-negative solution of system (1) that decays at infinity. We will determine the conditions of the existence of such stationary solutions. The detailed analysis of the reactiondiffusion type equations as well as the principles of formulating a mathematical model of virus infections and immunology are presented in [4] and [5].

Existence of stationary solutions
Stationary solutions of system (1) considered on the whole axis satisfy the following system of equations: where In order to find an analytical solution of this system, consider piece-wise constant functions σ i (x): when |y| ≥ y 1 0, when |y| < y 1 , σ 1,2 > 1. We look for a non-negative bounded solution of system of equations (2). We suppose for simplicity that α 1 = α 2 , β 1 = β 2 , γ 1 = γ 2 , and we omit the subscript in what follows. Clearly, non-negative bounded solutions of equations (2) can exist only if Then we can rewrite system (2) as follows: We look for its solutions in the form: where c 1 , c 2 , c 3 and c 4 are some non-negative constants and From the continuity of the solution and its first derivative at x = ±x 1 and y = ±y 1 , we obtain the following equalities: Dividing the equations at right by the equations at left, we get equations with respect to k: where . From equality (3) we can see that α = α(βI(u) + γI(v)) + k 2 , where I(u), I(v) > 0 and α, β, γ are positive constants. Therefore k 2 < α and we look for a solution k < √ α of this equation.
We can identify three cases:

Let u(x)
0, v(y) = 0. Then the value of k can be found from the first equality of system (5).
2. Let u(x) = 0, v(y) 0. Then the value of k can be found from the second equality of system (5).
Let us recall that we look for a solution k < √ α. In the case where only one component of the solution is different from 0, such solution exists if x * 1 > ξ 1 (in the first case) or y * 1 > ξ 2 (in the second case) is greater than the critical value and it does exist if x * 1 ≤ ξ 1 or y * 1 < ξ 2 , respectively in each case.

Let u(x)
0, v(y) 0. In this case the value k(< √ α) can be found from the first equality in (5). This value of k inserted in the second equality in (5) allows us to determine the relation between y * 1 and σ 2 for which there exists a solution (Figure 1, left). In the case where σ 1 and σ 2 are given, we can determine the relation between y * 1 and x * 1 from equations of system (5). If σ 1 = σ 2 , then y * 1 = x * 1 (black line in Figure 1, right). If σ 1 < σ 2 , then y * 1 > x * 1 and such relation exists only if σ 1 > k 2 . If σ 1 > σ 2 , then y * 1 < x * 1 and such relation exists only if σ 2 > k 2 . We can now determine the integrals I(u) and I(v): and From (6) and (7), taking into account the first relations in (4), we have: The coeffcients c 1 and c 3 can be determined from equation (3): where We can solve linear system of equations (8) with respect to c 1 and c 3 . Similarly, we can find the coefficients c 2 and c 4 . An example of the analytical solution is shown in Figure 2. We see that u(x) is positive and v(y) is negative. We have not found the values of parameters for which both components of the soltuion are positives.

Numerical simulations of system (1)
For numerical simulations, system (1) is considered in the interval [0; L] with Neumann boundary conditions. The system is discretized by using an explicit finite differences scheme. The functions σ 1 (x) and σ 2 (y) are piece-wise constant. The integral is approximated by the trapeze method. Figure 3 shows the behavior of the solution of system (1) for the values of parameters L = 12, α = β = γ = D 1 = D 2 = 1; σ 1 (x) = 0 for 3 < x < 9 and = 2 otherwise, initial condition = 1 for 5 < x < 7; σ 2 (y) = 0 for 2.637 < x < 9.363 and = 8 otherwise, initial condition = 1 for 5.5 < x < 6.5. The behavior of this solution is characterized by a slow convergence to a stationary solution for which u(x) > 0 and v(y) = 0. The analytical solution considered in the previous section shows the existence of a stationary solution with a positive component u(x) and a negative component v(y). However, since the initial condition is positive for both components, the solution of system (1) remains also positive and converges to another stationary solution having one positive component and one zero component.

Discussion
In this work, we study the competition of two viruses in the host organism. The model takes into account virus mutation, reproduction, and mortality. The nonlocal terms describe virus multiplication and competition for uninfected host cells. Virus mortality depends on the genotype and it can be due to antiviral treatment or to its natural death.
We consider virus strain as a density distribution concentrated around some genotype value. This is a non-negative solution decaying at infinity. The analysis of the model shows the existence of such solutions under some conditions on parameters. Namely, the admissible interval where virus multiplication rate is larger than its mortality rate should be sufficiently long, and the mutation rate should be small enough.
We identify three types of nonzero solutions: in two of them, one component of the solution is positive, and another one is zero. There is one more solution with one positive and one negative components. This solution does not have biological meaning, and it is not considered in the applications. Thus, the competition between two viruses results in the disappearance of one of them and persistence of another one. We have not found a solution with two positive components. This can be due to the assumptions of the equality of the coefficients α 1 = α 2 , β 1 = β 2 , γ 1 = γ 2 . More general case where these coefficients are different will be considered in the future works.
Let us note that variation of the death function in time allows the investigation of the influence of an antiviral treatment and virus evasion. In the case of one equation model this question was studied in [3]. For the two equation model considered in this work this question remains open.