INFECTION SPREADING IN CELL CULTURE AS A REACTION-DIFFUSION WAVE

. Infection spreading in cell culture occurs due to virus replication in infected cells and its random motion in the extracellular space. Multiplicity of infection experiments in cell cultures are conventionally used for the characterization of viral infection by the number of viral plaques and the rate of their growth. We describe this process with a delay reaction-diffusion system of equations for the concentrations of uninfected cells, infected cells, virus, and interferon. Time delay corresponds to the duration of viral replication inside infected cells. We show that infection propagates in cell culture as a reaction-diffusion wave, we determine the wave speed and prove its existence. Next, we carry out numerical simulations and identify three stages of infection progression: infection decay during time delay due to virus replication, explosive growth of viral load when infected cells begin to reproduce it, and finally, wave-like infection progression in cell culture characterized by a constant or slowly growing total viral load. The modelling results are in agreement with the experimental data for the coronavirus infection in a culture of epithelial cells and for some other experiments. The presence of interferon produced by infected cells decreases the viral load but does not change the speed of infection progression in cell culture. In the 2D modelling, the total viral load grows faster than in the 1D case due to the increase of plaque perimeter.


Infection progression in cell culture
Virus replication occurs inside host cells due to a complex multi-stage process resulting in production of new virions expelled from the infected cell into the extracellular space [3].Each cell can produce hundreds or even thousands of new virus particles often dying after that, which depends on cell type and on specific properties of viral infection.Let us note that some viruses, including the coronavirus SARS-CoV-2, develop special mechanisms influencing cell cycle and provide a more efficient replication [17].
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., [6,7,11,14,16] 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 it forms 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 modelling 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 [22] in the case of reversible host infection.The plaque growth rate is determined by the method of linearization (cf.Sect.3.1).Numerical simulations of viral plaque growth described by a reaction-diffusion model with time delay are presented in [10].This model is different in comparison with the model (1.1)-(1.3)considered below, and these different approaches are complementary (see Discussion for more detail).Existence of solutions of an initial-boundary value problem for a reaction-diffusion model of viral infection is proved in [12].Individual based models of viral infection spreading in cell cultures are developed in [1,15].
In this work we continue to study viral plaque growth with a delay reaction-diffusion system of equation, in line with the previous investigations but with a different model, which seems to us more appropriate, and more complete results.The model is presented in the next section and the main results in Section 1.3.

Model formulation
We consider the system of equations describing virus spread in cell culture.Here  is the concentration of uninfected cells,  of infected cells,  of virus,   () = (,  −  ).Below in this work we will also consider a more complete model taking into account interferon production by infected cells, and a 2D formulation of this model.If we neglect death of infected cells ( = 0), and we will see below that this assumption is appropriate for some experiments, then from equations (1.1), (1.2), we get  +  =  0 , where  0 is some constant.In this case we can reduce system (1.1)-(1.3) to the following system of two equations: =   2   2 +   − . (1.5) In the study of travelling wave solution, we substitute (, ) = ( − ),  (, ) = ( − ) and obtain the following system of equations: ′′ +  ′ + ( +  ) −  = 0 (1.7) with respect to functions () and () having limits at infinity: In the analytical part of the work we will study the existence of solutions of this problem, and the speed of wave propagation.Numerical simulations will be carried out for system (1.1)-(1.3) in a bounded interval, and for a more complete system containing also the interferon concentration.

The contents and main results
Viral infection can spread in cell culture as a reaction-diffusion wave.We will study in this work the properties of infection waves.One of their main characteristics is the speed of propagation.In order to determine it, we will use the linearization method applicable in the so-called monostable case where the waves exist for all speeds greater than or equal to the minimal speed.In general, this method allows the estimation of the minimal wave speed from below, that is,  * ≤  0 , where  * is the linearization speed and  0 is the minimal speed (see [20] and the references therein).However, it is not known a priori whether it gives the exact value of the minimal wave speed.For example, for a single reaction-diffusion equation with the logistic nonlinearity,  0 =  * , but for some other nonlinearities, it is possible to have the exact inequality  0 >  * .We will prove that for system (1.6), (1.7) the waves exist for all speeds greater than or equal to the linearization speed  * .Therefore, it gives indeed the minimal speed.In Section 2 we obtain this result for the model with cell diffusion, and in Section 3 without cell diffusion, passing to the limit as the diffusion coefficient tends to 0.
The wave existence is proved for the model without cell death ( = 0) where it can be reduced to a monotone system of two equations.Wave existence is known for monotone systems without time delay (see [20] and the references therein).Time delay in equation (1.5) enters in such a way that the system still satisfies the maximum principle and the comparison theorems allowing the proof of wave existence.Namely, we will use the method of upper and lower functions.However, the presence of time delay changes some spectral properties of the corresponding differential operators, and the wave existence is proved under the additional condition that  is limited from above.
The analytical results on the wave existence and the speed of propagation are completed by numerical simulations of system (1.1)-(1.3) in Section 4. They allow us to identify three stages of infection progression.In the beginning, during time delay in virus replication, the total viral load decreases, while the viral particles introduced initially diffuse in the culture infecting the surrounding cells.At the second stage, infected cells begin to produce new virus particles, and there is an explosive growth of their concentration.The third stage of infection progression is characterized by wave propagation along the cell culture.If infected cells die, then the total viral load remains constant.In the case without cell death, the total viral load grows linearly in time proportionally to the wave speed.We compared these numerical results with the experimental data from the literature for coronavirus and poliovirus [11,16].In the 2D simulations, the total viral load increases faster than in the 1D case because the rate of cell infection increases proportionally to the plaque perimeter.We complete the numerical part of this work with a model including the interferon production by the infected cells.Interferon slows down virus replication and decreases the total viral load.However, the speed of infection wave in the cell culture remains the same.

Existence of waves
Consider the problem with respect to functions () and () having limits at infinity: 3) It differs from the previous one by the presence of the diffusion term in the first equation.
We will study the existence of solutions of problem (2.1)-( 2.3) and the corresponding values of speed .As usual in the monostable case, there is a minimal wave speed  0 which will be determined from the linearized problem under the condition that its solution monotonically decreases at infinity.If the wave speed is less than  0 , the solutions of the linear problem oscillate and change sign.Therefore, they do not have biological significance.The determination of the minimal wave speed  0 is quite straightforward.However, the proof of wave existence for all  ≥  0 is much more involved.We will use the method of upper and lower functions, and we will begin with an auxiliary problem (2.4)-(2.6)presented below, where the term  with an unknown speed  in equation (2.2) is replaced by a given constant ℎ.First, we will prove the wave existence for ℎ = 0, then for ℎ > 0, and will conclude with ℎ =  .In Section 3, we will consider the limiting case  1 → 0 and obtain wave existence for problem (1.6)- (1.8).
Along with problem (2.1)-(2.3)consider the auxiliary problem ) with respect to functions  1 () and  2 () having limits at infinity: Here ℎ is a real constant.We will use the method of upper and lower functions constructed as solutions of the corresponding linear problem: ) with respect to functions  1 () and  2 () having limits at infinity: We look for its solution in the form assuming that , ,  > 0. Let us begin with the case  1 =  2 = 1 in order to get an explicit expression for the minimal wave speed.
Multiplying the first equation by , the second by  and subtracting, we get Hence, Since  and  are positive, we take the sign plus in this equality.We obtain from (2.11): In what follows we will denote the minimal speed by  0 .Notation  01 used here indicates that this value is obtained for  1 =  2 = 1.
Proof.We begin with the linearized system (2.7), (2.8).The function provides its solution converging to 0 as  → ∞ for any  1 and  2 .Let us note that  1 <  2 since we suppose that  >  0 .For any  1 and  2 such that Then (2.17) and (2.18) hold for any .If we decrease , then solution of the nonlinear system (2.4), (2.5) approaches the solution of the linear system (2.7), (2.8).Hence, there are solutions of system (2.4), (2.5) satisfying the assertion of the theorem.
In order to have a more explicit proof of the existence of such solutions, we will use the implicit function theorem.Let us write system (2.4), (2.5) as follows: ) (ℎ = 0).Then for  = 0 we obtain the linear system (2.7), (2.8), and for  = 1 the system (2.4), (2.5).We consider system (2.19), (2.20) on the half-axis  ≥ 0 with the boundary condition Then Thus, for all sufficiently small , we obtain a solution of system (2.4), (2.5) on the half-axis  > 0 satisfying the boundary condition (2.24).
It remains to verify that we can apply the implicit function theorem to problem (2.19)- (2.21).In what follows, we will consider the corresponding operator [19].We recall that the norm in these spaces are defined by the equalities First of all, we show that  = 0 does not belong to the essential spectrum of the operator , that is, this operator satisfies the Fredholm property.The essential spectrum is determined as the set of all complex  for which the problem considered on the whole axis, has a nonzero bounded solution.Assuming that  = 0 and applying the Fourier transform, we obtain the equalities where ṽ1 and ṽ2 are the Fourier transforms of the corresponding functions.Equating the determinant of the matrix of this system to 0, we get From the second equation we conclude that  = 0. Then the first equation leads to a contradiction since , ,  0 > 0. Thus, we conclude that the essential spectrum does not contain the origin.Next, consider the eigenvalue problem for equations (2.25), (2.26) on the half-axis  ≥ 0 with the boundary conditions  1 (0) =  2 (0) = 0.This problem has a unique (up to a constant factor) nonzero solution for  = 0. Indeed, it is sufficient to set  1 = − 2 in order to satisfy the boundary conditions.Moreover, a similar solution can be constructed for all sufficiently small real .
Let us recall that the index  of a Fredholm operator is defined as  =  − , where  is the dimension of the kernel of the operator and  is the codimension of the image, that is, the number of solvability conditions of the non-homogeneous problem.Both numbers  and  are constant inside each connected component of the complex plane where the operator satisfies the Fredholm property, except possibly for some isolated points (eigenvalues).We will show in Lemma 2.5 below that  = 0, that is, the non-homogeneous problem is solvable for any right-hand side.In this case, we can apply the implicit function theorem ( [20], Thm.1.20, Chap.4, p. 233).The lemma is proved.
Remark.In the construction of the lower function, it is sufficient that it is positive on some interval.However, it can be verified that it is positive on the half-axis  ≥ 0. Indeed, its positiveness on bounded intervals follows from the fact that it is close to the solution of the linear problem.Its behavior at infinity is determined by the exponential with the minimal exponent  1 .For a more detailed analysis, we can consider weighted spaces with the exponential weight exp( 1 ) and apply the implicit function theorem in the weighted space.

The case ℎ ̸ = 0
Substitute solution (2.10) into equations (2.7), (2.8): Equating the determinant of this system to 0, we obtain the equation with respect to : (2.32) Equation ( 2.32) with respect to  has from two to four solutions (Fig. 1, right).We are interested in its positive solutions in the interval 0 <  < / 1 .The positivity of  provides decay of functions (2.10) at infinity.From equation (2.30) we deduce that / = ( −  1 )/( 0 ).Therefore,  and  have the same sign if  belongs to the indicated interval.
We can now proceed to the construction of upper and lower functions.Lemmas 2.2 and 2.3 remain applicable.Therefore, we can affirm the existence of an upper function.Some changes should be done in the proof of Lemma 2.5.
Lemma 2.10.There exists ℎ 0 () > 0 such that for all ℎ ∈ (0, ℎ 0 ()] the problem ) ) Proof.Two properties should be verified: (1) the essential spectrum does not cross the origin if ℎ decreases from the given value to 0. Then the index of the operator  does no change, and  = 1 as for ℎ = 0.Moreover, its  and  characteristics can change only in isolated points, (2)  = 1.Then  = 0.The essential spectrum is determined as the set of all complex  for which the problem ) considered on the whole axis, has a nonzero bounded solution.Assuming that  = 0 and applying the Fourier transform, we obtain the equalities where ṽ1 and ṽ2 are the Fourier transforms of the corresponding functions.Equating the determinant of the matrix of this system to 0, we get There is such ℎ 0 > 0 that this system does not have solutions for any 0 ≤ ℎ < ℎ 0 .Indeed, if  0 ℎ < , then the second equation has a unique solution  = 0. Then the first equation leads to a contradiction.Hence, ℎ 0 ≥ /( 0 ).It remains to verify that the homogeneous problem (2.33)-( 2.35) has a unique linearly independent solution.Since equation (2.32) has three positive solutions  1 <  2 <  3 , general solution of system (2.33), (2.34) has the following form: where   /  = ( 2  −  )/( 0 ),  = 1, 2, 3. Since the right-hand sides of these equalities are different for different   , then there are two linearly independent vectors (  ,   ).From the boundary condition (0) = 0, we express  1 and  2 through  3 and obtain a unique linearly independent solution of the problem (2.33)-(2.35).Since  1 <  2 ,  3 , then the sign of the solution at infinity is determined by the sign of  1 (the same as  1 ).The lemma is proved.
Remark.Let us note that for ℎ sufficiently large, system of equations (2.38) can have a solution (Fig. 1, left) signifying that the essential spectrum passes through the origin.Moreover, equation (2.32) can have two positive solutions in the interval 0 <  <  (Fig. 1, right).From equations (2.38) we get the equality It has a single solution  2 depending on .Having found this solution, we can determine a sequence of values of ℎ for which equations (2.38) have solution.The minimal positive ℎ = ℎ 0 from this sequence satisfies the condition of the lemma.

The case ℎ = 𝑐𝜏
In this case we get the equation instead of equation (2.32).Since the function  (, ) increases with respect to  in the interval 0 <  < , and the exponential in the right-hand side decreases, then the conclusions above about the existence of the minimal speed  0 remain valid.Lemma 2.12.There exists a value  0 > 0 such that for  >  0 the system of equations (2.39) with respect to  has exactly three real positive solutions  1 ,  2 ,  3 ,  1 <  2 <  3 , and  1,2 ∈ (0, / 1 ).
We now obtain the following theorem as a corollary of Theorem 2.11.

Analytical estimate of the wave speed
Linearizing system (1.6), (1.7) at ∞, we obtain the system We look for the solution in the form We express  from equation (3. 3) and substitute in equation (3.4): and, hence, There exits  0 > 0 such that this equation has two solutions  for any  >  0 .The minimal wave speed  0 is given by the equality Example.We find the minima of these functions for different values of  and determine the wave speed: The corresponding numerical values are obtained in numerical simulations of system (1.1)-(1.3)presented in the next section.
We note from formula (3.6) that the wave speed depends on √ , as usually for reaction-diffusion problems.It depends on the product , and this dependence is weak.

Existence of waves for 𝐷 1 = 0
We begin with some estimates of solutions of problem (2.1)-(2.3) in order to consider the limiting case  1 → 0. We need to estimate the second derivative  ′′ independently of  1 .Consider this question in a more general setting.Proposition 3.1.Suppose that solution () of the system of equations on the whole axis, where  = ( holds for any fixed  ≥  > 0 with some constant  2 independent of   and depending on . Proof.Consider, for certainty, the first equation of system (3.7).Let  0 be the point of maximum of the function  ′ 1 ().Then  ′′ 1 ( 0 ) = 0, and | ′ 1 ( 0 )| ≤ /.Let us note that the solution is supposed to be bounded.Therefore, its derivative converges to 0 at infinity, at least along some sequence of .The same estimate can be obtained at each maximum of the derivative providing the estimate for all  ∈ R. Similar estimates can be obtained for other components of the vector  ′ ().Differentiating the first equation with respect to  and applying the same approach, we obtain the estimate of the second derivative.Theorem 3.2.Suppose that  ≥  0 , where  0 is defined by equality (3.6), and  ≤  0 (), where  0 (), is defined in Theorem 2.13.Then problem (1.6)-(1.8)has a monotonically decreasing (component-wise) solution.
Proof.Let us note that the minimal speed  0 =  0 ( 1 ) determined in Lemma 2.1 converges to the minimal speed  0 =  0 (0) in equality (3.6) as  1 → 0. Therefore, for any  >  0 (0) and  1 sufficiently small,  >  0 ( 1 ).The existence of waves for such  1 > 0 follows from Theorem 2.13.We will consider below the limit  1 → 0 for the fixed value of  and, consequently, prove the wave existence for  1 = 0. Having proved the wave existence for any  >  0 (0), we can conclude about the existence of waves for  =  0 (0) passing to the limit.In order to prove the convergence of solutions as  1 → 0, we will begin with the convergence on compact subintervals of the form [−, ], and pass to the limit  → ∞.
For  ∈ (0, 1], let  1 =  and (  ,   ) denote the family of solutions from Theorem 2.13.As the system does not depend explicitly on , then we can for definiteness fix its position in the -space by imposing the condition   (0) = 1 2  0 .Consider a sequence {  } =∞ =1 such that   > 0 and lim →∞   = 0 and the corresponding sequence of solutions (  ,   ) := (   ,    ),  = 1, 2, . . .From Proposition 3.1, it follows that these solutions are uniformly bounded in R together with their third derivatives independently of  1 .
Using the Arzela-Ascoli lemma we conclude that on every interval   = [−, ],  > 0, we can find a subsequence {  } =∞ =1 and the corresponding subsequence {   } =∞ =1 such that the functions (   ,    ) converge to some limiting function ( 0 ,  0 ) together with their second derivatives Since all functions in the sequence (   (),    ()) tend to the same limits as  → ±∞.This, together with the fact that their derivatives (up to the third order) are uniformly bounded implies the convergence to the limiting function in  2 (R).This limiting function satisfies problem (1.6)-(1.8).The theorem is proved.These initial conditions mean that cell culture initially contains only uninfected cells, and some initial viral load is introduced at the left part of the computational interval.Numerical implementation of the problem is presented in the appendix.Figure 2 shows an example of numerical simulations of system (1.1)-(1.3)with the spatial distributions of virus and infected cells in consecutive moments of time.The distance between the curves is not the same in the beginning of the simulation because of decaying oscillations due to time delay.Analysing the evolution of virus concentration distribution in space and in time, we can identify three stages of infection progression.At the first stage, virus does not yet reproduce due to the time delay in its replication.It diffuses in space while the total viral load decays because of its death.This first stage can be seen in Figure 4 as a linear decay of the function log   (), where   () = ∫︀  0  (, )d is the total viral load.This decay corresponds to exponential decrease of the function   ().At the next stage, initially infected cells begin to reproduce new virus particles leading to explosive growth of its total concentration.The last stage of infection progression corresponds to the propagation of infection front with a constant speed and linear growth of the total virus load.It can be seen as a logarithmic growth in Figure 4 (right).

Three stages of infection progression
If infected cells die after some time ( > 0), then their distribution has a spike-like shape, or a pulse moving along the cell culture (Fig. 3, upper left).Virus distribution has a similar form.In this case, the total virus concentration remains constant at the third stage of infection progression (Fig. 4, left).

Comparison with the experiments
Conventional multiplicity of infection (MOI) test is based on the experimental evaluation of virus spread in cell culture.We describe the experimental results on the spread of coronavirus in the culture of cilia cells [16] and of poliovirus [11].Virus spread in the cilia cells occurs preferentially through the apical surface.The dependence of virus titer on time   () has three stages: it decreases till time 10 h, then rapidly increases from 10 to 20 h, and finally its growth rate is slower after 20 h.The first stage corresponds to time delay of virus replication inside cells.New viruses are not yet produced, and the initial viral load decays due to virus death.Two other stages correspond to virus replication.During the second stage, the initial virus load (decreased in the first stage) converges to the wave.The last stage corresponds to the wave propagation.
The choice of the model and parameters.Let us note that there are two different wave propagation regimes from the point of view of virus production: with or without death of infected cells.If infected cells die ( > 0), then they form a spike of a constant size during the wave propagation.Therefore, the total virus production is also constant in this case,   () = ∫︀  0  (, )d → const.This case corresponds to some experimental data (Fig. 4, left) but not to coronavirus in the culture of cilia cells (Fig. 4, right).
Suppose now that infected cells do not die ( = 0).Then their concentration remains constant behind the wave ( =  0 ).Since the wave propagates with a constant speed, then the total number of infected cells grows linearly in time, () = , where  is the wave speed.Integrating equation (1.5) with respect to  over the interval [0, ], we obtain an approximate equation: Hence Thus,   () has asymptotically linear growth.This analysis shows that we should consider infection spreading without death of infected cells in order to describe the experimental results for coronavirus (Fig. 4, right) and with death of infected cell for poliovirus (Fig. 4, left).In the case of coronavirus, viral load decreases during first 10 h [16].Hence, we set  = 10 h.Virus diffusion rate varies in the limit from 10 −4 cm 2 /h to 10 −3 cm 2 /h.We take the maximal value  = 10 −3 cm 2 /h.We will discuss this choice below.The value of  is determined by the decay rate of   during the first stage.We get  = 0.1 ÷ 0.2 1/h.There are two free parameters  and .The rate of   growth during the second stage strongly depends on .The wave speed at the third stage depends on the product  but this dependence is weak.The initial viral load is taken from the experimental data.
Comparison of modelling and experiment.We describe this comparison in more detail in the case of coronavirus in the culture of cilia cells (Fig. 4, right).The function log   () in modelling shows the same three stages as in the experiment: exponential decay (linear in the log scale), exponential growth (linear in the log scale), linear growth (logarithmic in the log scale).The value   (0) is chosen in such a way that it coincides with the experimental value.The value   ( 1 ),  1 = 10 h is controlled by the value of .We set  = 0.1 in order to have the modelling value larger than the experimental one in order to minimize the difference between them for larger times.The growth rate at the second stage is determined by the parameter .We choose it in such a way that   ( 2 ),  2 = 20 h is larger than the experimental value in order to minimize the difference later.At the last stage, the growth rate in modelling is less than in the experiment.Since the wave speed weakly depend on parameters  and , their variation does not allow us to improve the accuracy.The wave speed is proportional to √ , but we take the maximal physically justified value of the diffusion coefficient.Therefore, we cannot use it to increase the wave speed.If we increase , then we can decrease the difference at  = 50 h but this will increase the difference at  = 20 h and  = 30 h.Thus, comparison of modelling and experiment show that growth rate in modelling at the last stage is less than in the experiment.Moreover, we can see from (4.2) that main term of log   () is log  in the 1D model.It is independent of the parameters of the model and cannot be increased for this model.We will see below that the growth rate is larger in the 2D model.Distributed delay.In the case of distributed delay in the virus replication rate, equation (1.3) in system (1.1)-(1.3) is replaced by the equation where This integral describes distributed delay with a uniform distribution in the interval [ −  − ℎ,  −  + ℎ].We recover the model with a point-wise delay in the limit of small ℎ.An example of numerical simulations of this model is shown in Figure 4 (right).The total viral load has qualitatively similar behavior as in the case of point-wise delay.However, virus reproduction begins earlier, and the total viral load becomes larger with the increase of ℎ.

Model with interferon
Infected cells produce interferon which downregulates virus replication.In this case, we complete system of equations (1.1)-(1.3)by the equation for the interferon concentration : ) The virus replication rate in equation (4.6) is inversely proportional to the interferon concentration .On the other hand, virus can slow down the interferon production [8,13].We take this effect into account in the interferon production rate in equation (4.7).This term contain time delay  2 which is different, in general, from the time delay in virus replication.Figure 5 shows an example of numerical simulations of system (4.4)-(4.7).Virus propagation is qualitatively similar to the case without interferon but the maximal virus concentration decreases.Interferon production follows the progression of infected cells.Its concentration strongly grows at the first stage of infection progression, then it propagates in space with damped oscillations.Let us note that the final value of virus concentration decreases due to the interferon production.However, the speed of infection progression, that is the speed of reaction-diffusion wave remains constant.Such behavior is specific for the monostable case for which the wave speed is determined by linearization at the virus free equilibrium (Sects. 2 and 3).The corresponding total virus load   () is shown in Figure 4 (right).

2D simulations
We present in this section the 2D modification of the system (1.1)-(1.3)by introducing the following equations for uninfected cells  , infected cells  and virus  : where (, ) = √︀  2 +  2 , and we set  = √︀  0 / in order to obtain the same initial viral load as in 1D.Let us note that we consider a smooth function  (, , 0) in order to avoid problems due to discontinuity of the solution on unstructured grid.Furthermore, we introduce a fine mesh near the radius  and a coarser mesh elsewhere (see (b) in Figs. 6, 8), and in order to not use a finer grid in all domain, we consider piecewise quadratic (P2) finite-elements instead of piecewise linear (P1) finite-elements.At the end, these simulations were done on a mesh with 79.468 triangles.We use the library PETSc8 (the Portable, Extensible Toolkit for Scientific Computation) which is a suite of data structures and routines for the scalable (parallel) solution of scientific applications modeled by partial differential equations.The computations were performed on a Macbook pro 2.2 GHz Intel Core i7, 16 GB of DDR4 2400 MHz RAM and took 00:24:16 over 6 threads.
We compute the viral load   () = ∫︀ Ω  (, , ) dd and we compare the results with the 1D results obtained by FreeFEM code with the same parameters described above, for  = 0 (Figs.6, 7) and  = 0.2 (Figs. 8, 9).In both cases, the total viral load is larger in the 2D case in comparison with the 1D case because the quantity of infected for the former grows in time due to the circular geometry.Though this difference does not seem essential in Figures 6c and 8c due to the log scale, the total viral load is about 3 times larger in the 2D case.

Discussion
Experimental studies of viral infection progression in cell culture allow the assessment of the efficacy and the mechanisms of virus penetration the host cells and its replication inside the cells, both processes being important for the determination of virus virulence.The duration of time delay in virus replication, clearly seen in the experiments, can depend on virus and host cells.This time delay can be taken into account in the optimization of anti-viral treatment.Furthermore, interferon produced by infected cells downregulates virus replication while some viruses (including coronavirus) can decrease the rate on interferon production [8,13].Such experiments represent a necessary stage for further investigation of viral infection in the human organism where it interacts with the immune system essentially increasing the complexity of this process.Mathematical modelling of infection development in cell cultures allows a better understanding of these processes and their quantitative description.
Mathematical models of viral infection in cell culture.Viral plaque growth in cell culture can be described as a reaction-diffusion wave which converts uninfected cells into infected and eventually dead cells due to virus production and diffusion.Surprisingly, there are only few works where this process is modelled, some of them with continuous models [10,22] and some other with individual based models [1,15].The advantage    of reaction-diffusion models is that they allow a theoretical investigation of infection waves, including their existence and the speed of propagation.The only delay reaction-diffusion model we are aware of differs from our model by the equation for infected cells [10].It contains the first time delay taking into account latently infected cells before they become infectious and the second time delay for the death of infected cells.This second (negative) delay term needs to be controlled to avoid the negativity of solutions.Contrary to this previous work, we consider time delay in the virus replication term and, this term being positive, it preserves the positiveness of solutions and the applicability of powerful mathematical methods based on the maximum principle and comparison theorems.We use them to prove the existence of waves.Let us note that this question is not only of pure mathematical interest because it allows us to determine the minimal wave speed, which is particularly important for the description of behavior of solutions.We discuss this question in the next paragraph.
Individual based models allow a more detailed description of intracellular regulation of virus replication and interferon production.It is also possible to consider more complete continuous models and hybrid discretecontinuous models [5].There is an increasing interest for all these modelling approaches in relation with the new coronavirus disease.
Existence and stability of waves.Existence of waves described by monotone reaction-diffusion systems is studied in the case without time delay (see [20] and the references therein).In general, the introduction of time delay essentially complicates the analysis of wave existence, as reviewed in [18].In some cases, when the maximum principle remains applicable, the wave existence can be proved by conventional methods, though it may be more technically difficult and some additional conditions can be necessitated to be imposed.As such, we prove the wave existence in the monostable case under the additional assumption that time delay is limited from above.This condition is related to the location of the essential spectrum, and we do not know whether the waves exist if it not satisfied.It is possible that oscillating waves emerge instead of monotone waves as it is the case for some other models [18].
Let us recall that the wave existence is proved by the method of lower and upper functions.It implies also the wave stability understood in a certain sense.If the decreasing solution provided by the upper function and the increasing solution provided by the lower function have the same limit, then the limiting function is unique, and all solutions between the upper and lower functions converge to the same wave.However, the uniqueness of the limiting function is not proved, so we have only a partial (one-side) stability result.Furthermore, the existence of waves for the minimal speed is proved passing to the limit of the waves with a larger speed.This method does not provide the stability of the wave with the minimal speed.However, we know for the scalar equation without time delay that the wave with the minimal speed is mostly interesting because solutions with initial conditions with a finite support converge to it.We observe this convergence numerically for the problem studied in the present work but neither local stability of the wave with the minimal speed nor global convergence to it are not yet proved.Stability of waves for delay reaction-diffusion systems requires further investigations.
Infection spreading in the tissue.Infection spreading in cell culture is accompanied by the interferon production by the infected cells, which is a part of the innate immune response.In the case of infection development in tissues of the human organism, there are also cells of innate immune response, such as macrophages and dendritic cells, and also the adaptive immune response determined by T and B lymphocytes.Both components of the immune response also involve numerous cytokines and chemokines.Altogether, the process of infection spreading becomes much more complex than in cell cultures.Some simplified models are considered in [2] for a single equation for virus concentration and in [4] for two equations for virus and immune cells.These models allow the determination of some regimes of infection spreading.The model developed in the present work provides an appropriate background for a more detailed and biologically justified investigation of viral infection development in the tissue of the organisms.This question will be studied in the future works.

A.3. 1D examples
The FreeFEM code was validated in the 1D case by the comparison with the C++ code based on implicit finite difference method (Fig. 3).We start our simulation by considering the 1D domain Ω = [0, ] and we took the following initial condition:  (, 0) =  0 = 1, (, 0) = 0 and  (, 0) =  0 for 0 ≤  <  0 and 0 otherwise.We note that we obtain similar between the two codes and a small difference for the CPU time: 11 s with FreeFEM code vs. 18 s with C++ code.

Figure 2 .
Figure 2. The distributions of virus concentration  (left) and of infected cells  (right) in numerical simulations of system (1.1)-(1.3).Different curves correspond to consecutive moments of time (every 2 time units).The values of parameters:  = 10,  = 5,  = 0,  = 0.4,  = 0.001,  = 2,  = 10,  0 = 1,  0 = 10,  0 = 0.1.The arrows in the left image show three stages of infection progression: virus concentration decay during time delay in its replication, explosive growth in its concentration when infected cells begin to reproduce new viral particles, propagation of infection wave along the cell culture.

Figure 3 .
Figure 3. Numerical simulations of system (1.1)-(1.3) in the case where infected cells die with certain rate.Concentration distributions (upper left) are different in comparison with Figure 2: the virus distribution (brown curve) and the concentration of infected cells (green curve) have a specific pulse-like form.The same virus distribution as a function of  and  is shown in the right figure.We can notice decaying oscillations of the maximal concentration.The corresponding total viral load log   () also oscillates (lower figure).The values of parameters are as follows:  = 10,  = 100,  = 1,  = 1,  = 0.001,  = 2,  = 5,  0 = 1,  0 = 10,  0 = 0.1.A small red interval on the -axis in the upper left figure shows the support of the initial viral load  0 .

Figure 6 .
Figure 6.Comparison between the results obtained in 1D (solid lines) and in 2D (dashed lines) with FreeFEM (a) normalized solution  (log scale),  and  at  = 25, (b) zoom on the center of the mesh used for the 2D case and (c) the total viral load log   () for  = 0.

Figure 7 .
Figure 7.The cut of the normalized solutions of  (a),  (b) and  (c) at  = 25 for  = 0.

Figure 8 .
Figure 8.Comparison between the results obtained in 1D (solid lines) and in 2D (dashed lines) with FreeFEM (a) normalized solution  (log scale),  and  at  = 25, (b) zoom on the center of the mesh used for the 2D case and (c) the total viral load log   () for  = 0.2.

Figure 9 .
Figure 9.The cut of the normalized solutions of  (a),  (b) and  (c) at  = 25 for  = 0.2.