Viral Infection Spreading and Mutation in Cell Culture

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


Introduction
Viruses constitute a separate specific kingdom of entities. They can have a plethora of geometrical forms, characterized by different shapes, symmetries and sizes, which cover the interval ranging from 20 nm (Porcine circovirus) up to 700 nm (Mimivirus). These properties depend mainly on the amount of nucleic acids they contain, their forms (single stranded positive-or negative-sense RNA, double stranded RNA, single stranded DNA or double stranded DNA), as well as on the physiology of the type of cells they use for their proliferation. Despite these essential differences, the survival and proliferation of all viruses (as quasi-species) is similar and based on the following strategy: enter the host cells, force the cells to multiplicate their genetic material (thus, in fact, to produce their multiple copies), leave the cells and invade new ones.
Experimental or clinical assessment of the progression of viral infection implies the evaluation of virus concentration in the infected tissue by means of conventional multiplicity of infection (MOI) assays (see, e.g., [1][2][3][4][5] and the references therein). After several consecutive dilutions, the virus-containing solution is poured in a cell culture leading to the formation of virus plaques. The number of such plaques determines the virus concentration measured in plaque forming units (PFU). The plaques consist of dead or modified cells, and they form due to virus replication inside the cells and its random motion (diffusion) between cells. The plaque growth rate characterizes the efficacy of virus penetration inside the cells, the intensity of its production and transmission between cells, thus, the virus virulence.
In spite of the importance and wide use of these experiments, there are relatively few modeling works devoted to the infection progression in cell culture taking into account spatial distributions of cells and virus particles. A reaction-diffusion model of plaque growth is considered in [6] in the case of reversible host infection. The plaque growth rate is 1.
Virus concentration decay during time delay in its replication in cells; 2.
Explosive growth of concentration, when infected cells begin to reproduce new viruses; 3.
Propagation of infection wave along the cell culture.
It was also found that the minimal speed of infection decreases with the delay time τ. This speed was effectively determined by finding the minimum of an explicitly constructed function containing the parameters of the model. However, the model proposed in [10] neglects other characteristic features of viruses survival strategy, namely their natural mutations and competition between their strains, which is crucial in the context of epidemiological analysis of virus infection spreading. An incorporation of these effects is the main objective of the present paper. We show that the model supplemented with the terms corresponding to mutation and competition phenomena, though significantly more complicated, can be studied by similar analytic methods and, together with the numerical calculations, can be a source of interesting information concerning the evolution of infection in the cell culture. The paper is organized as follows.
In Section 2, we analyze the model without mutation and competition terms and assuming that the mortality coefficients of the infected cells is nonzero, we show that the progression of infection is possible only if the virus replication number R v satisfies the inequality R v > 1. We also derive a lower bound for the minimal speed of the infection front. In Section 3, we propose a model taking into account mutation and competition phenomena. In Section 3.1, we show that, under some assumptions on the coefficients characterizing the considered strains of the virus, the new model can be reduced to one virus model (without mutation and competition). We also construct a function determining the lower bound of the infection propagation speed to the considered case. In Section 3.3, a system of reaction-diffusion equations describing multiple mutations and coinfections is proposed. In Section 4, we analyze different stages of infection progression in terms of time dependence of the total virus loads. The values and units of the assumed parameter values and their reference to culture experiments is summarized at the end of Section 4. Numerical simulations presented in this work are carried out using the finite element method with software freeFem++ [11]. Numerical accuracy of the simulations is controlled by decreasing the time and space discretization, and by the comparison with the analytical results. We begin with a conventional model of viral infection progression in cell culture without taking into account its spatial distribution:

Infection Progression with
Here, U(t) is the concentration of uninfected cells, I(t) of infected cells and V(t) of virus; a, b, β, σ are positive parameters. Parameter a characterizes infection rate of uninfected cells, b is the rate of virus replication in infected cells, β is the death rate of infected cells, and σ is the virus death rate. We suppose that Clearly, U = U 0 , I = 0, V = 0 is a stationary point of system (1)-(3). If β > 0, then it is a unique stationary point. In order to study its stability, we linearize the system about this stationary point and consider the corresponding eigenvalue problem: Denote the matrix in the left-hand side of this equality by A. Condition det A = 0 can be written as By analogy with epidemiological models, we will call it virus replication number (VRN). If R v < 1, then both eigenvalues of the matrix A are negative, and the stationary point is stable. If R v > 1, then this matrix has one positive eigenvalue, and the stationary point is unstable. In the first case, viral load V(t) decays, in the second case it grows. VRN represents a ratio of virus replication and elimination rates. Considering cell culture, we interpret parameters β and σ as death rates of infected cells and viruses. In the application to living tissue, which will be considered in the subsequent works, these parameters characterize the strength of the immune response with the elimination of infected cells by cytotoxic lymphocytes and virus neutralization by antibodies.

Final Concentrations of Cells and Total Viral Load
Next, assuming that R v > 1, we will determine the final value U f of uninfected cells as t → ∞. Taking a sum of Equations (1) and (2) and integrating from 0 to ∞, we obtain Integrating Equation (3), we have Finally, we divide Equation (1) by U and integrate: ln Excluding the integrals from Equations (5)-(7), we deduce the following equation where The derivation of Equation (8) is inspired by the derivation of the final size of epidemic in the model SIR, and the equation is the same, though the model is different. If we neglect the initial viral load and set V 0 = 0, then Equation (8) has a positive solution in the interval 0 < U < 1 if and only if R v > 1. This case corresponds to infection progression, and the final value of uninfected cells is less than its initial value. If V 0 > 0, then such solution exists for all values of R v .
Having found the solution of Equation (8), we can determine the total number of infected cells, I f = U 0 − U f and the total viral load

Model with Time Delay
If we take into account time delay in virus replication in the infected cells, then instead of Equation (3) we consider the equation where τ is a positive number. The initial condition for I(t) is now considered as I(t) = 0 for −τ ≤ t ≤ 0. In the stability analysis for system (1), (2), (9), we have det A = βσ − abU 0 e −λτ . The stability boundary, that is, the values of parameters for which λ = 0, remains the same as before, R v = 1. In the derivation of the final concentrations of cells, instead of Equation (6) we have now the equation then Equation (10) is reduced to Equation (6), and the formulas for the final call concentrations and total viral load remain valid for the model with time delay.

Spatial Infection Spreading
We continue with the model of viral infection progression in cell culture taking into account its spatial distribution [10]: The first term in the right-hand side of Equation (13) describes virus diffusion in the extracellular matrix, D is the diffusion coefficient.

Wave Propagation
Experimental data and previous modeling show [10] that infection spreads in cell cultures as a reaction-diffusion wave ( Figure 1). In order to study this solution, we consider a system of Equations (11)-(13) on the whole axis, and set U( where c is the wave speed. Then we obtain the following system of equations: where primes denote the derivatives with respect to the variable ξ = x − ct. We are interested in the solution of this system of equations with the following limits at infinity: In order to determine unknown value u f , we proceed as in the previous section. From Equation (14), we obtain taking a sum of (14), (15) and integrating: and from (16) Excluding the integrals from Equations (18)-(20), we obtain the equation with respect to ω = u f /u 0 . Equation (21) has a solution ω in the interval 0 < ω < 1 if and only if R v > 1. Thus, we obtain the following theorem.
Theorem 1. Inequality R v > 1 provides a necessary condition for the existence of a positive solution of problem (14)-(17). In this case, the final value u f can be found from Equation (21), and the total viral load Let us note that the notion of total viral load V X here is different from V T in the previous section. Here it signifies the total viral concentration in cell culture at any moment of time, while before it was the total viral load with respect to time (without space distribution). Furthermore, V X depends on the wave speed c, which is not yet found.

Model with Time Delay
In the model with time delay, Equation (13) is replaced by the equation and Equation (16) by the equation Clearly, integration of this equation gives (20), and all conclusions above remain applicable for this case.

Minimal Wave Speed
In order to determine the wave speed, we will use the linearization method widely used for monostable reaction-diffusion equations beginning from the KPP work [12]. The idea of the method consists to consider the system linearized at infinity and to find the minimal wave speed c * for which this system has a positive decreasing solution. Thus, instead of Equations (15) and (16) we consider the equations where u(ξ) is replaced by its value u 0 at infinity. We look for a solution of this system in the form with some positive numbers p, q and λ. Substituting these expressions into Equations (24) and (25) and excluding p and q, we obtain the following equation with respect to λ: We need to find a minimal positive value of c for which this equation has a positive solution. Set µ = λc. Then from the last equation we find Under the condition R v > 1, that is, σβ < abu 0 , function F(µ) becomes negative for small positive µ and tends to −∞ as µ µ 0 for a positive µ 0 for which the denominator vanish. Therefore, F(µ) > 0 for µ > µ 0 , F(µ) decreases in some interval µ 0 < µ < µ * and increases for µ > µ * . Hence, this function has a minimum for µ > µ 0 reached at µ = µ * . Set This equality provides the minimal wave speed found by the linearization method, and λ = µ * /c * . Let us note that, in general, this method gives an estimate of the minimal wave speed c 0 from below, c 0 ≥ c * . In some cases, it can be proved that c 0 = c * . This question is discussed in the next paragraph. Having found the wave speed, we can determine the total viral load V X in Theorem 1.

Wave Existence
If β = 0, then we can reduce system (11)-(13) to a system of two equations. Indeed, taking a sum of Equations (11) and (12), we conclude that U(x, t) + I(x, t) = U 0 , and the variable U can be excluded: This is a monotone system of equations for which it is possible to apply the maximum principle and various methods based on it. Existence of wave for this system is proved in [10] for all c ≥ c 0 for the values of τ limited from above. Thus, c 0 = c * for β = 0. If β > 0, this equality is verified in numerical simulations [10]. Knowing the value of speed, we can construct an approximate solution in the following way. Consider system (15), (16) with u = u 0 for ξ > 0 and u = u f for ξ < 0. Then we find the wave speed c * and solution w(ξ), v(ξ) for ξ > 0 as described in the previous paragraph. Let us note that this solution contains one unknown constant since there is a single relation between p and q. Next, we have a linear system for ξ < 0 with a given value c = c * . We can find its solution, and it depends on two unknown constants. These three constants can be determined from the continuity conditions at ξ = 0: . These calculations are quite cumbersome, and we do not present them here.

Viral Infection with Mutation
We will now consider the case where the first virus can mutate in another one. We consider the system of equations describing two viruses V 1 and V 2 in cell culture. Here U is the concentration of uninfected cells, I 1 of infected cells by virus V 1 and I 2 by virus V 2 , , is the mutation rate which shows that cells infected by virus V 1 can produce virus V 2 .
This theorem allows us to apply to system (29)-(33) the results of the previous section obtained for system (11)-(13) on wave existence and speed of propagation.

Choice between Three Solutions
System (29)-(33) can have one-virus solutions for which I 1 = 0, V 1 = 0 or I 2 = 0, V 2 = 0, and two-virus solution for which all components are positive. Conditions a i b i U 0 > β i σ i and the results of the previous section provide the existence of the first-virus solution for i = 1 and of the second-virus solution for i = 2. Conditions of the existence of two-virus solution are provided by Theorem 2. If these conditions are not satisfied, we can expect that the two-virus solution does not exist, though Theorem 2 does not state this non-existence result. It is confirmed by the results of numerical simulations. Analytical values of wave speed coincide with the wave speed in numerical simulations.

Method of Linearization
Reduction to the one-virus model is applied above under the assumptions D 1 = D 2 , β 1 = β 2 , σ 1 = σ 2 , τ 1 = τ 2 . If these conditions are not satisfied, we use the method of linearization to determine different regimes of infection spreading. We search solutions of system (29)-(33) in the form U(x, t) = u(x − ct), Then, After straightforward calculations, we obtain If q 1 = 0, then we deduce from (49) that Set µ = cλ. Then from the previous equality The wave speed c 0 is given by the equality In the particular case (50) and (51), (46): If q 1 = 0, then we have from (50) Then, Set µ = λc. The wave speed c 0 is given by the equality: It is similar to the previous one but with the parameters characterizing the second virus. From (53) we can conclude that a positive solution exists if a 1 b 1 > a 2 b 2 , that is, replication of the first virus is faster than of the second one. This condition appears also in Theorem 1. In the case of equality, the first virus is gradually replaced by the second one due to the mutation process ( Figure 2). If the inequality is opposite, the first virus vanishes faster since its replication is slower. We will return to this question in the next section.

Multiple Mutations and Coinfections
Generic model of virus mutation and coinfections can be written as follows: The coefficient b ij (positive or negative) show that viruses can promote or downregulate replication of other viruses in the case of coinfection. As before, under certain conditions on parameters, we can set I i = p i I, V i = q i V, and reduce this system of equations to the one-virus model. In the case of negative coefficients b ij the positiveness of solution in Equation (59) does not directly follow from the equation, and it should be controlled.

Single Virus
In the previous sections, we have studied reaction-diffusion waves of infection spreading in cell culture. The solution of the initial boundary value problem approaches the wave after the initial transient period. However, this initial period cannot be neglected in the investigation of viral infection since it determines further infection progression. It is particularly important from the point of view of interaction with the immune response, though we do not consider it in this work (see [13][14][15]).
In order to describe analytically the initial stage of infection development, instead of system (11)-(13), we consider the system ∂U ∂t where we replaced U by U 0 in Equation (61). This means that we neglect the depletion of uninfected cells in virus replication. This is a linear system which can be further simplified if we integrate the last two equations. In order to compare the results with the numerical simulations, we consider this system in a bounded interval 0 < x < L with the boundary conditions Using the notation J(t) = L 0 I(x, t)dx, W(t) = L 0 V(x, t)dx, we obtain the delay differential system of equations with the initial condition This system can be solved analytically. We restrict ourselves here to the two time intervals, t ∈ [0, τ): and t ∈ [τ, 2τ): where h(t) = (β − σ)(t − τ), used for the comparison with the numerical simulations. Figure 3 (left) shows the virus concentration distribution in consecutive moments of time with time difference between them equal 1 h. Let us note that since the length of the space interval is large and the diffusion coefficient is small, the descending part of the function V(x,t) looks steep in the graphical representation. It is smoother in Figure 4 and even more in Figure 1. Beyond the graphical representation, important question is about numerical accuracy. This descending part of solution contains about 130 discretization points, and the accuracy of the solution is verified.
The right panel of Figure 3 illustrates how the total viral load (the integral of virus concentration with respect to the whole space interval) changes in time. The total viral load is determined in direct numerical simulations of the space-dependent problem (11)- (13) and by two approximate analytical methods. The first one is based on solution of linear ODE (64), (65). It is applicable for a relatively small time (t < 4). The second analytical approximation uses Theorem 1 and the formula for the wave speed that allows us to determine the viral load in the propagating wave for t large enough.

Virus Mutation and Competition
In Section 3, we have established conditions of the coexistence of two viruses in cell culture. If these conditions are not satisfied, and the new variant replicates faster, it will eliminate the first one ( Figure 4).
Analysis of the total viral load ( Figure 5) shows that there are four stages of infection progression in this case instead of the three stages observed previously. After the infection decay and explosive growth, the distribution of the first virus reaches its bell-shape form and propagates along the culture. At the same time, the total viral load of the second virus variant exponentially grows. When it is reaches its maximum, the first viral load begins its exponential decay.
Denote by T the transition time defined as the moment of time at which V 2 (t) reaches 90% of its maximal value. Transition time depends on the mutation rate and on the replication rates determined by the parameters a i , b i , i = 1, 2 ( Figure 5, right). As it can be expected, T is larger for small , and it decreases if the replication rate of the second virus increases. The dependence of T on is well approximated by the function p/ q with some fitted values of p and q. Let us note that the for these values of q, the dependence of the transition time on is weak. Extrapolating the curves using the analytical approximation, we can see that decreasing the mutation rate by 10 6 increases the transition time only about twice.

Values of Parameters
The values of parameters for the single variant were determined in [10] by fitting the experimental data on the variation of viral load in time. Parameter estimation begins with the first stage of infection development during which virus is not yet reproduced by the infected cells, and its concentration exponentially decays. Duration of this stage and the decay rate allow us to uniquely determine parameters τ and σ by comparison with the total viral load in the experimental data. The next stage is characterized by explosive virus production with the maximum of the total viral load determined by the parameter b. At the same time, parameter a is chosen for the best fit of the growth curve. Finally, the third stage of infection development corresponds to its propagation in cell culture as a reaction-diffusion wave. The total viral load remains constant for β = 0, and it has a slow growth for a positive β, also determined by comparison with the total viral load in the experiment. Virus diffusion coefficient D is estimated from the literature. Let  The replication rate for the second virus (constants a 2 , b 2 versus a 1 , b 1 ) is taken in the same range or larger than that for the first virus. In order to evaluate the mutation rate , let us note that the rate of mutations of RNA viruses is from 10 −6 to 10 −4 per base per generation. Taking into account that only from 0.001 to 0.01 viruses represent plaque forming units, we estimate the probability of a given mutation during one virus replication of the order of 10 −9 to 10 −6 . The virus replication rate with respect to a single cell in the model is given by the constant b 1 = 8 × 10 4 . Multiplied by the probability determined above, we obtain the value of in the range from 10 −4 to 10 −1 . Further, we can reasonably assume that viruses of the same size have close diffusion coefficients. In particular, this is the case of two variants of the same virus. Therefore, we take the same values of diffusion coefficient for both viruses.

Discussion
Multiplicity of virus assays represent conventional tests for the assessment of virus virulence and other properties [3,5]. In a previous work [10], we have shown that infection spreading in cell culture can be described by a reaction-diffusion system with time delay.
Three stages of infection progression were identified: virus decay, explosive growth, wave propagation. In this work, we apply this modeling approach to study virus mutation during infection progression. This question has a particular importance in the context of the emergence of new virus variants in COVID-19 pandemic. In a simplified setting, without taking into account the influence of the immune response, we can model this process as infection progression in cell culture and analyze the conditions of the emergence of new virus strains.
We introduce in this work the notion of virus replication number R v (similar to the basic reproduction number in epidemiological models) and show that infection spreads in cell culture if this number calculated for the first virus is larger than 1. Furthermore, the total viral load can be expressed in terms of R v .
If the first virus variant spreads in cell culture and it can give a new virus variant due to mutations, then the two virus can co-exist or the second virus can eliminate the first one. The condition of their coexistence is given by the inequality a 1 b 1 > a 2 b 2 . Taking into account that the coefficients a i characterize the rate of cell infection and the coefficients b i the rate of virus replication in the infected cells, their product a i b i describes the total rate of virus production in cell culture. Hence, the two viruses coexist if the first one is produced faster than the second one. Otherwise, the second virus eliminates the first one.
The transition time from the first virus to the second virus depends on the mutation rate and on the replication rates. It is an important characterization of new variants since it shows whether they have enough time to develop in patients during the disease duration. The estimates of the transition time show that its dependence on the mutation rate is weak. Therefore, once new variant appears, it will rapidly replace the original one if it replicates faster, and the infected individual will become a career of the new variant.
This work is limited to the investigation of the emergence of new virus variants in the experimental setting of cell culture. The influence of the immune response on the infection mutation and spreading in the human organism will be considered in further works. However, some preliminary conclusions can be already done on the basis of the results of this work, in particular, about the action of vaccines on disease progression. In the case of SARS-CoV-2, the vaccine acts on virus spike preventing its penetration into the host cells. In the framework of the model, vaccine action can be described by decreasing the coefficient a. If it becomes small enough, such that the virus replication number R v is less than 1, then infection will not develop. If R v decreases but remains greater than 1, then infection will develop. However, infection spreading speed decreases. For the patients it signifies that disease symptoms are less severe.