Competition of SARS-CoV-2 Variants in Cell Culture and Tissue: Wins the Fastest Viral Autowave

Replication of viruses in living tissues and cell cultures is a “number game” involving complex biological processes (cell infection, virus replication inside infected cell, cell death, viral degradation) as well as transport processes limiting virus spatial propagation. In epithelial tissues and immovable cell cultures, viral particles are basically transported via Brownian diffusion. Highly non-linear kinetics of viral replication combined with diffusion limitation lead to spatial propagation of infection as a moving front switching from zero to high local viral concentration, the behavior typical of spatially distributed excitable media. We propose a mathematical model of viral infection propagation in cell cultures and tissues under the diffusion limitation. The model is based on the reaction–diffusion equations describing the concentration of uninfected cells, exposed cells (infected but still not shedding the virus), virus-shedding cells, and free virus. We obtain the expressions for the viral replication number, which determines the condition for spatial infection progression, and for the final concentration of uninfected cells. We determine analytically the speed of spatial infection propagation and validate it numerically. We calibrate the model to recent experimental data on SARS-CoV-2 Delta and Omicron variant replication in human nasal epithelial cells. In the case of competition of two virus variants in the same cell culture, the variant with larger individual spreading speed wins the competition and eliminates another one. These results give new insights concerning the emergence of new variants and their spread in the population.


Introduction
Viruses are non-cellular pathogens which use cell machinery to self-replicate, resulting in infection having to do with the population, the organism, the local tissue, or-in the simplest case-just the cell culture ex vivo or in vitro. In any case, replication of the virus is a "number game" between complex but still local biological processes (cell infection, virus replication in an infected cell, virus release, cell death, viral degradation-see Figure 1A) and spatial processes of viral transport between cells (convection, diffusion). Without transport, infection cannot proceed. In this study, we are interested in pure diffusion control of viral replication, which can be expected if convection (in particular, blood and lymph flow) is negligible, which is the case during viral replication in epithelium in vivo and cell cultures ex vivo/in vitro. COVID-19 pandemics made this question actual, as SARS-CoV-2 virus replicates in the upper respiratory tract (URT) and then in the lower respiratory tract (LRT), and corresponding ex vivo/in vitro models aimed at investigating and comparing new strains of this virus are being developing now [1,2]. Mathematical modelling is a powerful tool for understanding complex biological phenomena [3,4], and we use a mathematical approach here to study infection propagation in space. The SARS-CoV-2 virus enters the host from air mainly through cells of the URTolfactory and respiratory nasal epithelial cells-using the ACE2 (angiotensin-converting enzyme 2) receptor and TMPRSS2 (transmembrane serine protease 2) [5,6]. ACE2 expression decreases from URT towards LRT [7], thus passage of infection to lungs takes some time. However, most of the damage this virus causes is in the lungs, where severe COVID-19 results in diffusive alveolar damage and multiple thrombosis leading to respiratory failure [8]. Other tissues (nervous, cardiovascular, etc.) are also susceptible to this virus [5]. Thus, current experimental studies of SARS-CoV-2 replication involve multiple cell lines (in addition, from different species), but cells from URT and LRT are of the greatest concern [1,2]. In particular, comparison of different strains infecting these lines may produce insights concerning their relative virulence, severity, predictions of epidemic progression, etc.
The very recent SARS-CoV-2 variants of concern are Delta and Omicron [9]. Figure 2 shows the kinetics of replication of wild-type, Delta and Omicron variants in the URT and LRT cells reported by Hui et al. [1], Peacock et al. [2]. In both studies, a very interesting phenomenon was observed: Omicron got ahead of Delta in cells located above in the respiratory tract (A, C) during the first 1-2 days, but then, on the 2nd-3rd days, Delta overtook Omicron. In contrast, in cells located below in the respiratory tract (B, D), Delta got ahead during the whole time interval 1-3 days. While the second case seems to be obvious, the first one is not, and we aimed to study it closely. In particular, we are interested in the following questions regarding Figure 2A,C: • Can these experiments be described by the mathematical model (the simplest possible one)? • Could the situations of one strain overtaking another strain after some prolonged time interval be observed qualitatively if the spatial effects are taken into account? • Could spatial effects be very important: for example, could they result in any new information that cannot be obtained in the homogeneous system?
In spatially distributed systems with highly non-linear kinetics, one can expect completely new phenomena principally impossible in homogeneous conditions: stationary patterns and travelling fronts/waves of concentration distribution. Classical examples include animal coat patterns, epidemic progression across populations, nerve action potential propagation, Belousov-Zhabotinskii reaction, calcium waves in cells, etc. [4]. A decade ago, autowave properties were found in blood coagulation in no-stirring conditions [10] after their first prediction in 1994 [11] and subsequent extensive experimental and theoretical research which used mathematical models of different complexity, from two to dozens of variables (reviewed in [12]). Recently, autowave spatial progression of viral infection in tissues was studied theoretically by Bocharov et al. [13] using rather general mathematical models of two equations (for virus and immune cell concentrations). Mathematical models of infection propagation in tissues are based on the same approaches as models of epidemic progression across populations, models of chemical and biochemical concentration waves, etc.: these are reaction-diffusion models, but frequently with time delay because virus release from an infected cell occurs some finite time after cell infection (in the case of modelling the immune response with cell proliferation, corresponding terms frequently contain time delays, too). On the one hand, using time delays simplifies model equations, but on the other hand delays are the greatest drawbacks of such models as they hamper both analytical and numerical studies. Recent studies of autowave spatial infection spreading in cell cultures were also based on delay reaction-diffusion equations [14]. There, the conditions for infection propagation in terms of viral replication number and the formulas for the propagation speed were derived involving time delay as a parameter. In the present study, we use a similar approach, but (i) replace time delay with the intermediate cell sub-population-"exposed cells"-that is, infected but still not shedding the virus, and (ii) explicitly account for the decrease of free virus concentration concerning cell infection (which was neglected in previous studies). First, we obtain the new formulation for the viral replication number (which determines the condition for infection progression in space), estimates for the speed of spatial propagation of infection, and formulas for the final concentration of intact cells and the total (spatially-integrated) viral load. Second, we calibrate the model to recent ex vivo experimental data concerning replication of SARS-CoV-2 variants Delta and Omicron in cultures of human nasal epithelial cells reported by Peacock et al. [2]. Third, we numerically study spatial infection propagation in these cultures both in no-competition and in competition assumptions. We show that in the spatially distributed conditions, the Omicron variant, despite its much lower viral replication number, steady peak and total concentrations compared to those of Delta variant, can win the competition with Delta and completely suppress it.

Model Description and Governing Equations
Let us consider the scheme shown in Figure 1A. In the spatially distributed case, the system of equations corresponding to this scheme is: Here, variables U, E, and I stand for the concentrations of uninfected cells, "exposed" cells (infected but still not producing the virus), and infected cells (producing the virus), respectively, V stands for the concentration of virus particles in the medium, D is the diffusion coefficient of virus particles. N is the total number of virus particles produced by one cell during its life, N · β is the rate of virus production by one infected cell, ρ is the number of virions needed to infect one cell (usually only one virion is needed to infect a cell, but we keep ρ instead of 1 for generality and to study its influence on the results). Note that V and N account only for active virus particles in the medium, i.e., those that are able to infect a cell and self-replicate using it, but the total number of virus particles (which can be measured, for example, by PCR) can be much higher [2,15], see Section 2.4 below. The characteristic times of processes we account for are: cell infection, τ a = 1 aU 0 ; E → I transition, τ E = 1 γ ; cell death (and viral shedding), τ I = 1 β ; viral death, τ V = 1 σ . Boundary conditions for Equation (1d) are V x (±∞) = 0. Initial conditions for Equation (1) are where Ω 0 is some small but finite region onto which the virus has been inoculated (in calculations, we smooth the initial V(x, 0) distribution in Ω 0 , see Section 2.5). When considering the homogeneous case of Equation (1), we omit the diffusion term and the x coordinate, passing to the corresponding ODE system with initial conditions (U 0 , 0, 0, V 0 ).

Steady Travelling Solution
We are particularly interested in the steady travelling wave solution for Equation (1), which is established after some transition time period. Introducing the moving frame where c is the steady (i.e., established) wave speed, we denote u(ξ) = U(x, t), p(ξ) = E(x, t), q(ξ) = I(x, t), v(ξ) = V(x, t), where u, p, q, and v are independent on the new time τ. In new coordinates, Equation (1) become where prime denotes the derivative with respect to ξ. Conditions at infinity are Below in this article, we derive the condition for the infection progression, estimates for the steady wave speed, c, and formulas for the relative concentration of uninfected cells after the front, w = u f u 0 , and for the total virus concentration, v = +∞ −∞ v(ξ) dξ, as functions of parameters.

Two-Strain Model with Competition for Cells
The scheme shown in Figure 1B describes the coexistence of two virus strains, V 1 and V 2 , competing for the uninfected cells, U. The system of equations corresponding to this scheme is: Here, indexes 1 and 2 correspond to strains 1 and 2 (Delta and Omicron in this article), respectively. Diffusion coefficients for two strains are supposed to be equal. Initial conditions for each strain are just like those in the one-strain model, with equal initial distributions.

Model Parameters and Comparison with Experimental Data
In this article, we use two sets of model parameters: • Arbitrary non-dimensional values (typically, of the order of 1) for plotting the explanatory graphs when deriving analytical formulas in the Appendices. For this set, in numerical calculations we take L = 200, U 0 = 1. • Parameters which correspond to the SARS-CoV-2 variants Delta and Omicron replicating in human nasal epithelial cultures (hNECs) in the in vitro experiment [2] (see Table 1). We estimated them from physical reasons (e.g., D), derived/estimated from experimental results and kept fixed (e.g., V 0 , τ a ), or varied under the strict limitations as described below: 1.
The diffusion coefficient D was estimated using the Stockes-Einstein formula at 300 K assuming virus diameter 100 nm [16] and water viscosity.

2.
The experimental in vitro system in [2] is considered as homogeneous, because (a) only average concentrations are presented in that article, (b) the full size of the experimental system is 6.5 mm [17], which is comparable to the width of the front in our numerical calculations (about 1-3 mm, see Figure 5 below), and (c) substantial convection should be expected during inoculation and everyday sampling. 3.
Since the virus concentration determined by RT-qPCR for the E gene in [2] (see Figure 2C,D) was approximately 1000 times greater than the PFU concentration for all available time points and cell lineages (crossed vs. filled markers), that is only 1 out of 1000 viral particles was able to effectively infect a cell and reproduce with its help, we compared the model variable V(t) with the united data for PFU mL (filled markers) and E genes mL /1000 (open markers). In particular, we estimated V 0 as E genes Initial cell concentration U 0 was limited in the range [10 4 , 10 6 ] cells mL . 5.
To avoid large differences between the corresponding parameters of Delta and Omicron variants, the rations β 1 β 2 , γ 1 γ 2 , σ 1 σ 2 and N 1 N 2 were limited to the range of from 1 10 to 10. 7.
After the fitting, all parameters were rounded up to 1 decimal digit, and those that were close to each other were set as equal. 9.
For the obtained set of kinetic parameters, in spatial numerical calculations (with diffusion) we set L = 2 cm or more to be able to track the transition of autowave to the steady propagation regime.
Parameter estimation in the homogeneous case. We used COPASI 4.34 [18] to estimate the model parameters U 0 , γ 1,2 , β 1,2 , σ 1,2 , and N 1,2 from the experimental data (Section 2.4). During these estimations, each initial value problem was solved using the LSODA method (relative tolerance 10 −6 , absolute tolerance 10 −8 ). The parameter estimation task used the Differential Evolution algorithm (1000 generations, population size = 256) to minimize the sum of absolute differences of log(V) in the model and experiment at time points 0, 24, 48, and 72 h. Constraints for the parameter values and for their ratios described in Section 2.4 were applied. To fulfil the limitations on Delta/Omicron parameter ratios, two uncoupled versions of model (1)-for Delta and for Omicron strains-were solved simultaneously. (1) and (3), we used the previously developed package [10,19] based on the numerical methods described in [20]: Störmer-Encke's method for space discretization (uniform mesh with 401 nodes, zero Neumann boundary conditions at x = 0 and x = L) and the embedded Runge-Kutta-Fehlberg method of order 2(3) with automatic step size control for integration in time (mixed local error estimation with max norm, tol = 0.001, f ac = 0.8, f acmax = 5). The activation zone was Ω 0 = [0, L 10 ], and V 0 in this zone was multiplied by 2.257 · e −( x L/20 ) 2 to reduce gradients at t = 0, x = L 10 (thus, the average value of V(x, 0) in Ω 0 equaled V 0 ). The accuracy of the numerical method was controlled in two ways: (i) by solving the reduced Nagumo equation, u t = u 2 (1 − u) − hu + u xx , followed by comparison of the steady speed with the theoretical value, c theor = [3 (1 − 4h) − 1]/2 √ 2; relative discrepancy was <0.2% in the whole range of h ∈ [0, 0.2] allowing autowave solution; (ii) by repeating some calculations of model (1) with ten times lower tolerance and two times lower spatial mesh size. For example, for Omicron parameters (Table 1) and ρ = 1 relative changes of the steady speed and of the spatial-integrated virus concentration were 0.36% and 0.38%, respectively (relative to c ≈ 0.0062 cm/h and v ≈ 0.93 · 10 6 cells/cm 2 , see Table 2).  (1) is in the steady state (ss) (U, E, I, V) ss = (U 0 , 0, 0, 0). To study the stability of this steady state, we linearized equations around it (in the homogeneous case, neglecting the diffusion term). The stability is determined by the following Jacobian:

Solution to the reaction-diffusion problems. To integrate Equations
Here, we defined the virus replication number (VRN) as R v can be thought of as the ratio of virus production and death rates at the leading edge of the front, where U ≈ U 0 . Note that R v does not depend on the time of E → I transition, τ E = 1/γ. The characteristic equation |J ss − λE| = 0 reads as This equation for λ has one positive solution if R v > 1; otherwise, all its solutions are non-positive (see Appendix A.1). The existence of a positive λ means exponential growth of any perturbation of the steady state. Thus, the condition R v > 1 determines the instability of the steady state (U 0 , 0, 0, 0) of the homogeneous system, and infection progression in the spatially distributed system.

Estimates of the Steady Wave Speed
The speed of infection propagation in the steady travelling regime, c, is the eigenvalue to system (2). Analogous to the KPP problem, solutions exist for c in some range bounded from below, so the minimum possible speed should be found. After linearizing Equation (2) at ξ → +∞, where u → u 0 , we obtained the following estimate for the minimal wave speed (see Appendix A.2): The function F(µ) is negative at 0 < µ < µ 0 and positive at µ > µ 0 , where µ 0 is a discontinuity point, and the above equation corresponds to the minimum of the positive branch F(µ). We also obtained two explicit estimates of the minimal wave speed (see Appendix A.2): and, more precisely, where We compare these estimates with each other and with the numerical solutions in Appendix A.2 ( Figure A5) for arbitrary values of parameters. After we calibrate the model to experimental data, we use Formulas (6) and (7b) for comparison with the numerical results for Delta and Omicron spatial propagation.

Equations for the Final Concentration of Intact Cells and the Total Spatial Viral Load
In spatially distributed conditions, concentrations of all components depend on spatial coordinates. In particular, in the autowave regime, virus concentration is zero before and after the autowave, and passes through the maximum in the wave region; the concentration of intact cells in the wave region (i.e., the region of nonzero virus concentration) drops from the initial to some final value (see Figure 5 below). Depending on aims and facilities, various concentration measures can be used, for example maximum and total virus concentration in the sample, final concentration of intact cells, and its ratio to the initial value. The relation between these values and the wave speed can be obtained by integration of Equation (2) over space (see Appendix A.3): Here, w = u f u 0 ≤ 1 is ratio of the final to the initial concentration of intact cells, v is total (spatially-integrated) virus concentration. If R v > 1, the transcendental Equation (8a) has two roots: 1 and w * , where 0 < w * < 1 (Figure 3). Having found numerically the second root, w * , one can determine the total spatial viral load v using Equation (8b) provided that the wave speed c is found experimentally or theoretically (say, from Equation (6)). At R v >> 1, the root w * ≈ 0, so Equation (8) simplify to According to Equations (8) and (9), three main characteristics of infection, R v , c, and v, appear to be closely linked (we consider a fixed in this study to simplify the analysis). However, they express different properties; that is why stains of the same virus may differ in any/all of them in various combinations, leading to unexpected consequences. This is shown below based on the comparison with the experimental data.

Homogeneous Case without Competition: Omicron is "Quick" and Wins the Start but Delta Can Overtake It after 1-2 Days' Lag
To study the case when Omicron wins the start but then Delta overtakes it (Figure 2A,C), we calibrated the homogeneous version of our model (Equation (1) with D = 0) to one of these experiments ( Figure 2C). Figure 4 (left) and Table 1 show results of this calibration, and Figure 4 (right) shows time dependencies of other variables (namely, concentrations of cells) in the same solutions. Obtained values of parameters suggest that Omicron is "weaker": an infected cell produces five times less intact virions, and virion lifetime is six times smaller than for Delta. However, Omicron is "quicker": its times of E → I transition and of virion release are five times smaller than for Delta. This quickness explains the fast start of Omicron compared to Delta in homogeneous conditions (Figures 2A,C and 4). Later, after the transition period (1-2 days), Delta wins in homogeneous conditions (the same figures), apparently due to a greater N value. It is interesting to study the impact of spatial effects on this "number game", because non-homogeneous conditions can be easily expected both in vivo (say, in the respiratory tract) and in vitro (say, in the plaque assay, or in other specially designed setups).  Table 1.

Spatially Distributed Case without Competition: Omicron Can Win the Race despite Low Concentration and R v
Solutions to the same equations with the same parameters but in the spatially distributed conditions are presented in Figure 5A. We show results for two strains on the common graph for comparison (but note that equations for two strains are uncoupled). The values of the steady autowave speeds and concentrations of virus are given in Table 2. Compared to Delta, Omicron has six times lower peak concentration and 26 times lower spatially integrated concentration. However, Omicron propagates 1.27 times faster: the steady speed is 0.0062 cm/h compared to 0.0049 cm/h for Delta. The kinetics of spatially integrated concentrations ( Figure 5B) resemble homogeneous kinetics, where Omicron wins the start but then Delta overtakes it (Figure 4), but concentrations actually do not matter here: Omicron ultimately wins because it conquers more space, which means a greater number of cells or tissue volume infected. Note that the R v value for Omicron is also much lower than for Delta (Table 1). Table 2 also allows one to compare the steady speeds with those estimated analytically from Equations (6) and (7b) for ρ =1 and 0, and Figure 6 shows the gradual effect of ρ on speeds in numerical calculations and in Equation (6). For ρ = 1, estimation by Equation (6) gives lower values than numerical calculations; however, at ρ = 0, estimations and numerical calculations nearly coincide. We suggest that the discrepancy at ρ = 1 is due to overestimation of the rate of virus concentration decrease upon cell infection in the linearized model which was used to derive the analytical estimation (see the Section 4).   Table 1 for other model parameters). Speed values for ρ = 1 and ρ = 0 are presented in Table 2. Figure 5C. Equations for two strains are coupled now through the concentration of uninfected cells (the model suggests that the a cell cannot be co-infected by both strains). The calculation was stopped when the Omicron wave had travelled the same distance as in the no-competition situation (which is shown in Figure 5A). Compared to that case, the concentration of the Omicron strain is the same (solid line, note the difference in the y-scale), but the Delta strain is completely suppressed (dashed line). The steady speed of Omicron strain wave propagation is the same as in the no-competition situation, 0.0062 cm/h. The spatially integrated concentration of the Delta strain is lower than for the Omicron strain during the whole time period ( Figure 5D).

Discussion
In this study, we use the simplest possible mathematical model of viral replication able to account for the delay of virus shedding (through introduction of the sub-population of "exposed" cells) and for the spatial non-uniformity of concentration distribution of all components. A typical solution in this system has autowave behaviour: being initially stimulated near one of the boundaries, the wave of virus concentration moves away from that site and finally goes into an autowave mode characterised by constant shape and velocity ( Figure 5A). Spatial profiles of concentrations of exposed and infected cells also have the wave shapes, and the profile of uninfected cell concentration has a switch-off front shape ( Figure 5A, insert). For values of parameters listed in Table 1, the width of the forward wave front is about 1-3 mm, so the characteristic distances on which the wave should be studied are centimeters (we set L ≥ 2 cm). Experimental conditions (system size, convection, mixing, viscosity) and mode of virus transport between cells should alter the front width. For example, virus diffusion in a highly viscous medium or spreading through the cell-to-cell junctions should have a much lower effective diffusion coefficient, and thus the width of the wavefront should be smaller, leading to autowave behaviour detectable in smaller distances. In contrast, mixing or small size of the system should lead to spatial uniformity of concentrations, that is, to the homogeneous system lacking any spatial effects. That is why we use the homogeneous version of our model for calibration to experimental data obtained for hNECs in the MucilAir system (Figure 4, left).
We focus on viral replication in hNECs because these cells are important for SARS-CoV-2 entrance into the host as well as for the subsequent disease transmission, and because of intriguing kinetics reported in these cells for Delta and Omicron variants ( Figure 2C). Experimental results show, and the homogeneous model is able to reproduce, that Omicron wins at the start (1st day), but in the next 1-2 days Delta overtakes it (Figure 4, left). After describing the homogeneous kinetics of Delta and Omicron replication in hNECs, we study the spatial propagation of infection in these cells in no-competition and in competition assumptions. Without competition, Omicron starts faster and its wave runs quicker, however with much lower maximal and total concentrations than Delta (Figure 5, top). When competition is allowed, Omicron ousts Delta and completely suppresses it ( Figure 5, bottom). Thus, Omicron is "weaker but quicker", it wins the competition with Delta using higher spatial speed.
These results are of great interest because qualitatively different properties of virus replication in cell culture or tissue can be distinguished: viral replication number R v , viral concentration, either maximal v max or total (spatially integrated) v, and speed of spatial propagation of infection c. These properties appear to be closely linked in our simple model, see Equations (8) and (9). However, their "roles", or interpretations, are qualitatively different. Consequently, correct understanding of these properties is necessary for designing appropriate experimental conditions and ways of processing experimental results.
From the analysis of Equation (5), we conclude that R v determines the very possibility of spatial propagation of infection. This conclusion seems unexpected because R v does not depend on either D or γ, see Equation (4). Mahiout et al. [14] obtained the same condition for autowave propagation, R v > 1, but their model included time-delay instead of "exposed" cells, and R v was also independent on the delay. This means that general results are robust concerning variations in particular model details. R v turned out to be not the parameter that determines the fate of spatial competition between strains. The same conclusion follows for the concentration of virus (both maximal and total): ability to obtain high concentration may not be important in competition ( Figure 5). It is the speed of spatial propagation that determines the virulence and the fate of the competition between strains. The fastest strain conquers all the spatial area/volume of the culture or tissue.
Future studies should clarify whether this is a general conclusion for viral strain competition in tissues besides SARS-CoV-2 variants and hNECs. Another question which should be studied is the dependence of the competition fate on the initial conditions. As new strains appear due to mutations, from the very beginning they should compete with the old strain. Initial concentrations and spatial distributions should vary depending on location in the respiratory tract, person's time from infection, etc. We suppose that the R v value may play a role in the survival of a new strain if this new strain is a small local patch surrounded by the old strain.
In this study, we were also faced with an interesting theoretical result: taking into account the relatively subtle process-the decrease of free virus concentration due to binding to the cell it infects, −ρ · aUV term in Equation (1d)-not only slightly decreases the speed of autowave propagation, but also decreases the estimated autowave speed relative to the one found numerically ( Figure 6). The first effect is obvious, as any decrease of virus concentration should weaken the autowave. The second effect might not be so simple. It is not the numerical error, see Section 2.5. Thus, it could be some property of linearization of equations consisting in the substitution of u(ξ) with u 0 , which is the basis for steady speed estimation (see Appendix A.2). Linearization increases the absolute value of the auv terms in the equations for concentrations of exposed cells and virus because u 0 > u(ξ). Moreover, u(ξ) rapidly drops in the wave front ( Figure 5); thus, this replacement could have a great effect. However, linearization does not have any valuable effect on the speed at ρ = 0, and the discrepancy increases with increasing ρ. Thus, the effect is likely due to the over-estimated rate of virus elimination in the linearized system. Indeed, from our parameters for Delta (Table 1), the rate constant for virus degradation is σ 1 = 0.01 h −1 , but the quasi-1st order rate constant for virus binding, ρaU, at ρ = 1 and U = U 0 is aU 0 = 1 τ a = 1 h −1 . Thus, at ρ > 0, the linearization greatly increases the local virus elimination rate. The resulting effect on speed estimation is not great (Figure 6), but it should be noted that this effect is not attributed to numerical issues.

Conclusions
New experimental methods of viral investigation considering spatial effects, namely involving diffusion transport limitation and aimed at detecting the speeds and velocities of moving concentration waves, should be developed to obtain new insights concerning fundamentals of infection progression in tissues as well as between individuals. Practical application of these new methods will include facilitation of strain comparison and discovering new treatment options.
Eliminating p, q and v, we obtain the 4th degree equation for λ: If c > 0, the LHS of this equation has four real roots (see Figure A2, black solid line and square markers): At zero wave speed (c = 0), Equation (A2) reduces to the equation Figure A2, dotted line). Equation (A3) has solutions λ 00  At λ = 0, the LHSs of Equations (A2) and (A3) equal −γβ(σ + ρau 0 ). For Equation (A3), this is the minimum LHS value, and it is greater than the RHS if R v > 1 (note that both the LHS and the RHS are negative). Consequently, for R v > 1 and small c, say c ∈ [0, c min ), Equation (A2) also does not have real solutions. Increasing c shifts the parabola peak right and down, and there exists a minimum value of the wave speed, c = c min , for which Equation (A2) obtains a solution at λ = λ 0 > 0 (solid line, solid circle in Figure A2). For c > c min there are two solutions (dashed line). Thus, there exists a lower bound c min for possible wave speeds in solutions to the system (A1).  Figure A2. Notches on the µ-axis denote discontinuity points µ = µ 0 , closed circles denote minimum of the positive branch.
To find c min , we use two approaches: (1) Let us denote µ = cλ in Equation (A2): Rearranging this equation, we obtain: The function F(µ) is plotted in Figure A3 (for several ρ values). This figure shows that there exists such µ 0 > 0 that F(µ) |µ>µ 0 > 0. Moreover, the positive branch of the F(µ) graph has a single minimum. Consequently, which gives the exact estimate of the minimum wave speed in the system (A1).
(2) Let us note the approximate symmetry of the 4th order parabola near its peak, i.e., between λ 3 and λ 4 ( Figure A2, black solid line in proximity to the red line). From the expression for λ 3,4 , we conclude that Substituting this λ 0 into Equation (A2), we obtain Graphical solution to this equation is plotted in Figure A4 (solid lines). LHS (black solid line) has three negative roots, RHS = y 1 is constant (red line), and we look for their intersection at c 2 > 0 (circled region). At c 2 = 0, the LHS equals y 0 = 2γβ(σ + ρau 0 ). As y 0 Thus, a positive solution c 2 to cubic Equation (A7) exists if R v > 1. We can find the 1st and 2nd order estimates to this solution if we replace the positive branch (at c 2 ≥ 0) by the tangent line at c 2 = 0 or the parabola, respectively (black dotted and dashed lines): where k 1 = [LHS] |c 2 =0 = γβ + 2(γ + β)(σ + ρau 0 ), Solving these linear and quadratic equations, we get two estimates for c 2 min = 2D · c 2 : c * * 2 min = Dk 12 where Figure A4. Solution to Equation (A7) (solid lines) and two its approximations in the c 2 ≥ 0 region giving c * 2 (dotted line) and c * * 2 (dashed line). Parameters are the same as in Figure A2. Actual differences between solutions are hardly distinguished at these parameter values, however at larger R v all differences increase (see Figure A5).
To compare these analytical results with each other and with numerical solution to the initial system (1), let us take values of parameters listed in the Figure A2 (1) has steady speed c steady = 0.37. Consequently, the estimate by Equation (A5) is quite precise, and both Equation (A8) overestimate the speed. Here, the values of parameters corresponded to R v = 1.2 which is only a little above the boundary of autowave existence, R v = 1. Comparison in a wide range of R v values is given in Figure A5 for ρ = 0 (A) and ρ = 1 (B). Numerical solutions and estimates by Equation (A5) are close; however, discrepancy does exist and increases with increasing R v . Estimates by Equation (A8b) are correct for R v from 1 to approximately 10, and estimates by Equation (A8a) are correct only for very small R v above 1. Thus, the estimate min[F(µ)] gives the most precise speed value. Note that R v is actually the parameter of autowave existence (Section 3.1), and other parameters independently stay in formulas for F(µ), c * min , and c * * min , so these ranges of R v values where estimates work well correspond to the given particular set of parameters.  Figure A2.