<figure><img src="image_1.png" /></figure>

HAL
open science

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

Alexey Tokarev, Anastasia Mozokhina, Vitaly Volpert

► To cite this version:

Alexey Tokarev, Anastasia Mozokhina, Vitaly Volpert. Competition of SARS-CoV-2 Variants in Cell Culture and Tissue: Wins the Fastest Viral Autowave. Vaccines, 2022, 10 (7), pp.995. 10.3390/vaccines10070995. [https://doi.org/10.3390/vaccines10070995](https://doi.org/10.3390/vaccines10070995)

HAL Id: hal-04134103

https://hal.science/hal-04134103v1

Submitted on 10 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Article

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

Alexey Tokarev<sup>1,2,*</sup> <sup>1 ID</sup>, Anastasia Mozokhina<sup>1 ID</sup> and Vitaly Volpert<sup>1,3</sup>

<figure><img src="image_4.png" /></figure>

check for
updates

**Citation:** Tokarev, A.; Mozokhina, A.; Volpert, V. Competition of SARS-CoV-2 Variants in Cell Culture and Tissue: Wins the Fastest Viral Autowave. *Vaccines* **2022**, *10*, 995. https://doi.org/10.3390/vaccines10070995

Academic Editors: Khalid Hattaf and Bapi Pahar

Received: 23 May 2022

Accepted: 14 June 2022

Published: 22 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_5.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** 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.

**Keywords:** viral infection; competition of viral strains; reaction-diffusion systems; autowaves; mathematical modeling and analysis

## 1. 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

<footnote n="1">¹ S.M. Nikolskii Mathematical Institute, Peoples Friendship University of Russia (RUDN University), 6 Miklukho-Maklaya Street, 117198 Moscow, Russia; mozokhina-as@rudn.ru (A.M.); volpert@math.univ-lyon1.fr (V.V.)</footnote>

<footnote n="2">² N.N. Semenov Federal Research Center for Chemical Physics RAS, 4 Kosygina Street, Building 1, 119991 Moscow, Russia</footnote>

<footnote n="3">³ Institut Camille Jordan, UMR 5208 CNRS, University Lyon 1, 69622 Villeurbanne, France</footnote>

<footnote>* Correspondence: alexey.a.tokarev@gmail.com</footnote>
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.

<figure><img src="image_6.png" /><figcaption>**Figure 1.** Schemes depicting the transitions of cells between different infected states (black arrows), virus production (red arrows), and death of cells and virus (gray arrows). U—uninfected cells; E—“exposed” cells (infected but still not shedding the virus); I—infectied cells shedding the virus; V—free virus particles. (A) One-strain model. (B) Two-strain model with competition of strains for the uninfected cells.</figcaption></figure>

The SARS-CoV-2 virus enters the host from air mainly through cells of the URT—
olfactory and respiratory nasal epithelial cells—using the ACE2 (angiotensin-converting
enzyme 2) receptor and TMPRSS2 (transmembrane serine protease 2) [5,6]. ACE2 ex-
pression 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 se-
vere 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 com-
pletely new phenomena principally impossible in homogeneous conditions: stationary
patterns and travelling fronts/waves of concentration distribution. Classical examples in-
clude 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.

<figure><img src="image_7.png" /><figcaption>**Figure 2.** Viral replication kinetics: (A,B) in bronchus and lung cells [1]; (C,D) in nasal epithelial and human lung cells [2] (data were extracted from Figure 1 in preprint of [1] and from Figure 1 in preprint of [2]). (C) Open markers show E gene concentration (crossed markers) divided by 1000 to approximate match the corresponding PFU concentration.</figcaption></figure>

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.

## 2. Methods

### 2.1. 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:

$$ \frac{dU}{dt} = -aUV, \qquad (1a) $$

$$ \frac{dE}{dt} = aUV - \gamma E, \qquad (1b) $$

$$ \frac{dI}{dt} = \gamma E - \beta I, \qquad (1c) $$

$$ \frac{dV}{dt} = -\rho \cdot aUV + N \cdot \beta I - \sigma V + D \nabla^2 V. \qquad (1d) $$

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 \cdot \beta$ is the rate of virus production by one infected cell, $\rho$ is the number of virions needed to infect one cell (usually only one virion is needed to infect a cell, but we keep $\rho$ 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, $\tau_a = \frac{1}{aU_0}$; $E \to I$ transition, $\tau_E = \frac{1}{\gamma}$; cell death (and viral shedding), $\tau_I = \frac{1}{\beta}$; viral death, $\tau_V = \frac{1}{\sigma}$. Boundary conditions for Equation (1d) are $V'_x(\pm\infty) = 0$. Initial conditions for Equation (1) are

$$ U(x, 0) = U_0, \quad E(x, 0) = 0, \quad I(x, 0) = 0, \quad V(x, 0) = \begin{cases} V_0, & \text{if } x \in \Omega_0, \\ 0, & \text{otherwise,} \end{cases} $$

where $\Omega_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 $\Omega_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$).

### 2.2. 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

$$ \zeta = x - ct, \quad \tau = t, $$

where $c$ is the steady (i.e., established) wave speed, we denote $u(\zeta) = U(x, t)$, $p(\zeta) = E(x, t)$, $q(\zeta) = I(x, t)$, $v(\zeta) = V(x, t)$, where $u, p, q$, and $v$ are independent on the new time $\tau$. In new coordinates, Equation (1) become
$$-cu' = -auv, \qquad (2a)$$

$$-cp' = auv - \gamma p, \qquad (2b)$$

$$-cq' = \gamma p - \beta q, \qquad (2c)$$

$$-cv' = -\rho \cdot auv + N \cdot \beta q - \sigma v + Dv'', \qquad (2d)$$

where prime denotes the derivative with respect to $\xi$. Conditions at infinity are

$$u(+\infty) = u_0 \equiv U_0,$$

$$u(-\infty) = u_f,$$

$$p(\pm\infty) = q(\pm\infty) = v(\pm\infty) = v'(\pm\infty) = 0.$$

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 = \frac{u_f}{u_0}$, and for the total virus concentration, $\tilde{v} = \int_{-\infty}^{+\infty} v(\xi) d\xi$, as functions of parameters.

## 2.3. 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:

$$\frac{dU}{dt} = -a_1UV_1 - a_2UV_2, \qquad (3a)$$

$$\frac{dE_1}{dt} = a_1UV_1 - \gamma_1 E_1, \qquad (3b)$$

$$\frac{dE_2}{dt} = a_2UV_2 - \gamma_2 E_2, \qquad (3c)$$

$$\frac{dI_1}{dt} = \gamma_1 E_1 - \beta_1 I_1, \qquad (3d)$$

$$\frac{dI_2}{dt} = \gamma_2 E_2 - \beta_2 I_2, \qquad (3e)$$

$$\frac{dV_1}{dt} = -\rho_1 \cdot a_1UV_1 + N_1 \cdot \beta_1 I_1 - \sigma_1 V_1 + D\nabla^2 V_1, \qquad (3f)$$

$$\frac{dV_2}{dt} = -\rho_2 \cdot a_2UV_2 + N_2 \cdot \beta_2 I_2 - \sigma_2 V_2 + D\nabla^2 V_2. \qquad (3g)$$

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.

## 2.4. 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$, $\tau_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/mL ($t = 0$)/1000.

4. Initial cell concentration $U_0$ was limited in the range [10⁴, 10⁶] cells/mL.

5. The characteristic time of infection ($\tau_a$) was assumed to be 1 h for both strains. Characteristic times if $E \to I$ transition (1/γ₁,₂) and cell death (1/β₁,₂) were limited as 2 h and 0.1 h from below, respectively. Virus death times (1/σ₁,₂) were limited in the range [0.1, 100] h.

6. To avoid large differences between the corresponding parameters of Delta and Omicron variants, the ratios β₁/β₂, γ₁/γ₂, σ₁/σ₂ and N₁/N₂ were limited to the range of from 1/10 to 10.

7. Parameters $U_0$, γ₁,₂, β₁,₂, σ₁,₂, and N₁,₂ were determined by simultaneous best-fitting of two homogeneous versions of system (1)—for Delta and for Omicron variants—to the experimental data represented in Figure 2C from [2].

8. 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.

<table><caption>**Table 1.** Model parameters for Delta and Omicron variants.</caption><thead><tr><th>Parameter</th><th>Dimension</th><th>Delta</th><th>Omicron</th></tr></thead><tbody><tr><td></td><td></td><td colspan="2" style="text-align: center;">Equal for both strains:</td></tr><tr><td>V<sub>0</sub></td><td>virions/mL</td><td></td><td>5 &times; 10<sup>3</sup></td></tr><tr><td>U<sub>0</sub></td><td>cell/mL</td><td></td><td>1 &times; 10<sup>6</sup></td></tr><tr><td>&tau;<sub>a</sub></td><td>h</td><td></td><td>1</td></tr><tr><td>a &equiv; 1/&tau;<sub>a</sub>&middot;U<sub>0</sub></td><td>cell&middot;h/mL</td><td></td><td>1 &times; 10<sup>-6</sup></td></tr><tr><td>D</td><td>cm<sup>2</sup>/h</td><td></td><td>1 &times; 10<sup>-4</sup></td></tr><tr><td>&gamma;, &beta;</td><td>h<sup>-1</sup></td><td>&gamma;<sub>1</sub> = &beta;<sub>1</sub> = 0.04</td><td>&gamma;<sub>2</sub> = &beta;<sub>2</sub> = 0.2</td></tr><tr><td>&sigma;</td><td>h<sup>-1</sup></td><td>&sigma;<sub>1</sub> = 0.01</td><td>&sigma;<sub>2</sub> = 0.06</td></tr><tr><td>N</td><td>1 &equiv; virions/cell</td><td>N<sub>1</sub> = 50</td><td>N<sub>2</sub> = 10</td></tr><tr><td>R<sub>v</sub> = N<sub>1</sub>-1/&tau;<sub>a</sub>&sigma;</td><td>1</td><td>R<sub>v,1</sub> = 4900</td><td>R<sub>v,2</sub> = 150</td></tr></tbody></table>

## 2.5. *Numerical Methods*

1. *Parameter estimation in the homogeneous case.* We used COPASI 4.34 [18] to estimate the model parameters $U_0$, γ₁,₂, β₁,₂, σ₁,₂, and N₁,₂ from the experimental data (Section 2.4). During these estimations, each initial value problem was solved using the LSODA method (relative tolerance 10⁻⁶, absolute tolerance 10⁻⁸). 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.
2. **Solution to the reaction-diffusion problems.** To integrate Equations (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, fac = 0.8, facmax = 5). The activation zone was $\Omega_0 = [0, \frac{L}{10}]$, and $V_0$ in this zone was multiplied by $2.257 \cdot e^{-(\frac{x}{L/20})^2}$ to reduce gradients at $t = 0, x = \frac{L}{10}$ (thus, the average value of $V(x, 0)$ in $\Omega_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\sqrt{(1-4h)-1}/2\sqrt{2}; 2\sqrt{2}]$; relative discrepancy was <0.2% in the whole range of $h \in [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 $\rho = 1$ relative changes of the steady speed and of the spatial-integrated virus concentration were 0.36% and 0.38%, respectively (relative to $c \approx 0.0062 \text{ cm/h}$ and $\tilde{v} \approx 0.93 \cdot 10^6 \text{ cells/cm}^2$, see Table 2).

<table><caption>**Table 2.** Autowave speed $c$ (cm/h), spatially integrated viral concentration $\tilde{v}$ (10⁶ cells/cm²) and maximal viral concentration $v_{max}$ (10⁶ cells/mL) for Delta and Omicron steady spatial autowave propagation.</caption><thead><tr><th></th><th>Delta</th><th>Omicron</th><th>Delta</th><th>Omicron</th></tr><tr><th></th><th colspan="2">ρ = 1</th><th colspan="2">ρ = 0</th></tr></thead><tbody><tr><td>Numerical c</td><td>0.0049</td><td>0.0062</td><td>0.0094</td><td>0.011</td></tr><tr><td>c<sub>min</sub> as √min[F(μ)], Equation (6)</td><td>0.0038</td><td>0.0058</td><td>0.0093</td><td>0.011</td></tr><tr><td>c from Equation (7b)</td><td>0.0069</td><td>0.0088</td><td>0.016</td><td>0.015</td></tr><tr><td>Numerical ṽ</td><td>24</td><td>0.93</td><td>46</td><td>1.88</td></tr><tr><td>ṽ = (√min[F(μ)]/a) R<sub>v</sub>, Equations (6) and (9b)</td><td>19</td><td>0.87</td><td>45</td><td>1.65</td></tr><tr><td>Numerical v<sub>max</sub></td><td>26</td><td>4.3</td><td>27</td><td>4.9</td></tr></tbody></table>

### 3. Results

#### 3.1. Virus Replication Number Provides the Condition for the Infection Progression, $R_v > 1$

At $x \to +\infty$, system (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:

$$J_{ss} = \begin{vmatrix} -\gamma & 0 & aU_0 \\ \gamma & -\beta & 0 \\ 0 & N \cdot \beta & -\rho \cdot aU_0 - \sigma \end{vmatrix} = (R_v - 1)\gamma\beta\sigma = \frac{R_v - 1}{\tau_E \tau_I \tau_V}.$$

Here, we defined the virus replication number (VRN) as

$$R_v = \frac{aU_0}{\sigma}(N - \rho) = \frac{\tau_V}{\tau_a}(N - \rho). \qquad (4)$$

$R_v$ can be thought of as the ratio of virus production and death rates at the leading edge of the front, where $U \approx U_0$. Note that $R_v$ does not depend on the time of $E \to I$ transition, $\tau_E = 1/\gamma$. The characteristic equation $|J_{ss} - \lambda E| = 0$ reads as

$$\begin{vmatrix} -\gamma - \lambda & 0 & aU_0 \\ \gamma & -\beta - \lambda & 0 \\ 0 & N \cdot \beta & -\rho \cdot aU_0 - \sigma - \lambda \end{vmatrix} = 0,$$

which gives

$$(\lambda + \gamma)(\lambda + \beta)(\lambda + \sigma + \rho \cdot aU_0) = N \cdot aU_0\gamma\beta. \qquad (5)$$
This equation for $\lambda$ has one positive solution if $R_v > 1$; otherwise, all its solutions are non-positive (see Appendix A.1). The existence of a positive $\lambda$ 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.

### 3.2. 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 $\xi \to +\infty$, where $u \to u_0$, we obtained the following estimate for the minimal wave speed (see Appendix A.2):

$$c_{min}^2 = \min_{\mu > \mu_0} F(\mu), \text{ where } F(\mu) = \frac{D\mu^2}{\mu + \sigma + \rho au_0 - \frac{\gamma\beta}{(\mu + \gamma)(\mu + \beta)} Nau_0}. \quad (6)$$

The function $F(\mu)$ is negative at $0 < \mu < \mu_0$ and positive at $\mu > \mu_0$, where $\mu_0$ is a discontinuity point, and the above equation corresponds to the minimum of the positive branch $F(\mu)$. We also obtained two explicit estimates of the minimal wave speed (see Appendix A.2):

$$c_{min}^{*2} \approx 4D \frac{R_v - 1}{\frac{1}{\sigma} + 2(\frac{1}{\gamma} + \frac{1}{\beta})(1 + \frac{\rho au_0}{\sigma})}, \quad (7a)$$

and, more precisely,

$$c_{min}^{**2} \approx Dk_{12} \left[ \sqrt{1 + \frac{8}{k_{12}} \cdot \frac{R_v - 1}{\frac{1}{\sigma} + 2(\frac{1}{\gamma} + \frac{1}{\beta})(1 + \frac{\rho au_0}{\sigma})}} - 1 \right], \quad (7b)$$

where

$$k_{12} = \frac{k_1}{k_2} = \frac{\gamma\beta + 2(\gamma + \beta)(\sigma + \rho au_0)}{\gamma + \beta + 2\sigma + 2\rho au_0}, \text{ and } \frac{\rho au_0}{\sigma} \equiv \frac{R_v}{\frac{N}{\rho} - 1}.$$

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.

### 3.3. 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):

$$\ln w = -R_v \cdot (1 - w), \quad (8a)$$

$$\tilde{v} = -\frac{c}{a} \ln w. \quad (8b)$$

Here, $w = \frac{u_f}{u_0} \le 1$ is ratio of the final to the initial concentration of intact cells, $\tilde{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 $\tilde{v}$ using Equation (8b) provided that the wave speed $c$ is found experimentally or theoretically (say, from Equation (6)). At $R_v \gg 1$, the root $w^* \approx 0$, so Equation (8) simplify to

$$ \ln w \approx -R_v, \qquad (9a) $$

$$ \tilde{v} \approx \frac{c}{a} R_v. \qquad (9b) $$

According to Equations (8) and (9), three main characteristics of infection, $R_v$, $c$, and $\tilde{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.

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Graphical solution to Equation (8a) for two $R_v$ values. Thin line—ln $w$, bold lines—$R_v \cdot (w - 1)$.</figcaption></figure>

### 3.4. 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 \to 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).
<figure><img src="image_9.png" /><figcaption>**Figure 4.** Left: Replication of the SARS-CoV-2 variants Delta (green + red) and Omicron (blue + cyan) in hNECs in the experiment [2] (markers, data are taken from Figure 2C as described in Section 2.4) and in the homogeneous version of the model (1) (lines). Right: concentrations of cells in the same solutions. Figures are screenshots from COPASI. Model parameters are listed in Table 1.</figcaption></figure>

3.5. *Spatially Distributed Case without Competition: Omicron Can Win the Race despite Low Concentration and R<sub>v</sub>*

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*<sub>v</sub> 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).
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Spatial replication and propagation of Delta (dashed lines) and Omicron (solid lines) variants without competition for cells (A,B) and with competition (C,D). (A) Two uncoupled models given by Equation (1) were solved simultaneously. (C) Equation (3) were solved. (A,C) Solutions corresponding to the same distance (1.5 mm) travelled by Omicron autowave $V_2(x)$ are shown. (B,D) Logarithms of integrated virus concentrations vs. time are shown. $\rho = 1$, all other parameters are listed in Table 1.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Steady autowave speed as a function of $\rho$, assuming that $\rho$ can vary continuously (see Table 1 for other model parameters). Speed values for $\rho = 1$ and $\rho = 0$ are presented in Table 2.</figcaption></figure>
### 3.6. Spatially Distributed Case with Competition: Omicron Can Win and Completely Suppress Delta

Solution to Equation (3) describing the competition of two virus strains is presented in 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).

## 4. 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 \ge 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) $\tilde{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, $-\rho \cdot 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(\xi)$ 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(\xi)$. Moreover, $u(\xi)$ 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 $\rho = 0$, and the discrepancy increases with increasing $\rho$. 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 $\sigma_1 = 0.01 \text{ h}^{-1}$, but the quasi-1st order rate constant for virus binding, $\rho aU$, at $\rho = 1$ and $U = U_0$ is $aU_0 = \frac{1}{\tau_0} = 1 \text{ h}^{-1}$. Thus, at $\rho > 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.

**5. 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.

**Author Contributions:** Conceptualization and methodology, A.T. and V.V.; software, investigation and visualisation, A.T.; funding acquisition and writing, A.T., A.M., V.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is supported by the Russian Science Foundation project №22-21-00830, https://rscf.ru/en/project/22-21-00830/ (accessed on 13 June 2022).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.
**Conflicts of Interest:** The authors declare no conflict of interest.

**Abbreviations**

The following abbreviations are used in this manuscript:

<table><tr><td>URT</td><td>upper respiratory tract</td></tr><tr><td>LRT</td><td>lower respiratory tract</td></tr><tr><td>hNECs</td><td>culture of human nasal epithelial cells</td></tr><tr><td>KPP or KPPF equation</td><td>Kolmogorov–Petrovskii–Piskunov–Fisher equation</td></tr></table>

**Appendix A**

Appendix A.1. Condition for Existence of a Positive Solution to Equation (5)

The condition $R_v > 1$ for the existence of a positive solution to Equation (5) can be obtained using Figure A1, where the LHS (left-hand side) of this equation is plotted. The graph intersects the y-axis at $y = \gamma\beta(\sigma + \rho aU_0)$. Thus, a positive solution exists if $\gamma\beta(\sigma + \rho aU_0) < NaU_0\gamma\beta \Leftrightarrow \sigma < (N - \rho)aU_0 \Leftrightarrow R_v = \frac{aU_0}{\sigma}(N - \rho) > 1$.

<figure><img src="image_12.png" /><figcaption>**Figure A1.** LHS of Equation (5) for $\gamma = 1.5$, $\beta = 2$, and $\sigma + \rho aU_0 = 1$.</figcaption></figure>

Appendix A.2. Estimates of the Steady Wave Speed

Our estimations of $c_{min}$ are based on linearization of Equation (2) at $\xi \to +\infty$, where $u \approx u_0$. Substituting $u(\xi)$ with $u_0$, we have:

$$ 
\begin{aligned}
-cp' &= au_0v - \gamma p = \frac{v}{\tau_a} - \frac{p}{\tau_E}, \\
-cq' &= \gamma p - \beta q = \frac{p}{\tau_E} - \frac{q}{\tau_I}, \\
-cv' &= -\rho \cdot au_0v + N \cdot \beta q - \sigma v + Dv'' = -\rho \frac{v}{\tau_a} + N \frac{q}{\tau_I} - \frac{v}{\tau_v} + Dv''.
\end{aligned}
\quad (A1) $$

The right-hand sides (RHSs) of these equations are written first using the kinetic constants and then using the characteristic times. The solution of this system is close to the solution of system (2) in the leading front of the autowave. Let us look for the solution in the form $p = p_1e^{-\lambda\xi}$, $q = q_1e^{-\lambda\xi}$, $v = v_1e^{-\lambda\xi}$, where $\lambda > 0$. Substituting these conjectures into Equation (A1), we obtain:
$$c\lambda p = au_0v - \gamma p = \frac{v}{\tau_a} - \frac{p}{\tau_E},$$

$$c\lambda q = \gamma p - \beta q = \frac{p}{\tau_E} - \frac{q}{\tau_I},$$

$$c\lambda v = -\rho \cdot au_0v + N \cdot \beta q - \sigma v + D\lambda^2 v = -\rho \frac{v}{\tau_a} + N \frac{q}{\tau_I} - \frac{v}{\tau_v} + D\lambda^2 v.$$

Eliminating $p, q$ and $v$, we obtain the 4th degree equation for $\lambda$:

$$(c\lambda + \gamma)(c\lambda + \beta)[D\lambda^2 - c\lambda - \sigma - \rho \cdot au_0] = -N \cdot au_0\gamma\beta. \quad (A2)$$

If $c > 0$, the LHS of this equation has four real roots (see Figure A2, black solid line and square markers):

$$\lambda_1 = -\frac{\gamma}{c} = -\frac{1}{c\tau_E},$$

$$\lambda_2 = -\frac{\beta}{c} = -\frac{1}{c\tau_I},$$

$$\lambda_{3,4} = \frac{1}{2D} \left[ c \mp \sqrt{c^2 + 4D(\sigma + \rho au_0)} \right] = \frac{1}{2D} \left[ c \mp \sqrt{c^2 + 4D(\tau_v^{-1} + \rho \tau_a^{-1})} \right].$$

At zero wave speed ($c = 0$), Equation (A2) reduces to the equation

$$\gamma\beta \cdot [D\lambda^2 - \sigma - \rho \cdot au_0] = -N \cdot au_0\gamma\beta \quad (A3)$$

with roots of the LHS $\lambda_{3,4}^0 = \mp\sqrt{\frac{\sigma+\rho \cdot au_0}{D}} = \mp\sqrt{\frac{\tau_v^{-1}+\rho \tau_a^{-1}}{D}}$ (Figure A2, dotted line). Equation (A3) has solutions $\lambda_{3,4}^{00} = \mp\sqrt{\frac{\sigma}{D}(1-R_v)}$, which requires $R_v \le 1$ for $\lambda_{3,4}^{00}$ to be real. Thus, for $R_v > 1$, Equation (A3) does not have real solutions (dotted line does not intersect the red one).

<figure><img src="image_13.png" /><figcaption>**Figure A2.** Solution of Equations (A2) and (A3). Black lines—LHSs, red line—RHSs. $\gamma = 1.5$, $\beta = 2$, $\sigma = 1$, $D=1$, $R_v = 1.2$, $\frac{N}{\rho} = 10$, and $c$ values are specified.</figcaption></figure>

At $\lambda = 0$, the LHSs of Equations (A2) and (A3) equal $-\gamma\beta(\sigma + \rho 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 \in [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 $\lambda = \lambda_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><img src="image_14.png" /><figcaption>**Figure A3.** Function $F(\mu)$ from Equations (A4) and (6) at various $\rho$ values. $N = 10$, other parameters are the same as in Figure A2. Notches on the $\mu$-axis denote discontinuity points $\mu = \mu_0$, closed circles denote minimum of the positive branch.</figcaption></figure>

To find $c_{min}$, we use two approaches:

(1) Let us denote $\mu = c\lambda$ in Equation (A2):

$$ (\mu + \gamma)(\mu + \beta)\left[D\frac{\mu^2}{c^2} - \mu - \sigma - \rho au_0\right] = -N \cdot au_0\gamma\beta. $$

Rearranging this equation, we obtain:

$$ c^2 = \frac{D\mu^2}{\mu + \sigma + \rho au_0 - \frac{\gamma\beta}{(\mu+\gamma)(\mu+\beta)}Nau_0} = \frac{D\mu^2}{\mu + \sigma(1 + \frac{R_v}{\frac{N}{\rho}-1}) - \frac{\gamma\beta}{(\mu+\gamma)(\mu+\beta)}\sigma\frac{R_v}{1-\frac{N}{\rho}}} = F(\mu). \quad (A4) $$

The function $F(\mu)$ is plotted in Figure A3 (for several $\rho$ values). This figure shows that there exists such $\mu_0 > 0$ that $F(\mu)|_{\mu>\mu_0} > 0$. Moreover, the positive branch of the $F(\mu)$ graph has a single minimum. Consequently,

$$ c_{min}^2 = \min_{\mu>\mu_0} F(\mu), \quad (A5) $$

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 $\lambda_3$ and $\lambda_4$ (Figure A2, black solid line in proximity to the red line). From the expression for $\lambda_{3,4}$, we conclude that

$$ \lambda_0 \approx \frac{\lambda_3 + \lambda_4}{2} = \frac{c_{min}}{2D}. \quad (A6) $$
Substituting this $\lambda_0$ into Equation (A2), we obtain

$$ \left( \frac{c_{\min}^2}{2D} + \gamma \right) \left( \frac{c_{\min}^2}{2D} + \beta \right) \left( D \frac{c_{\min}^2}{4D^2} - \frac{c_{\min}^2}{2D} - \sigma - \rho \cdot au_0 \right) = -N \cdot au_0 \gamma \beta. $$

Denoting $c_2 = \frac{c_{\min}^2}{2D} > 0$, we have

$$ (c_2 + \gamma)(c_2 + \beta)(c_2 + 2\sigma + 2\rho au_0) = 2N \cdot au_0 \gamma \beta. \quad (A7) $$

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\gamma\beta(\sigma + \rho au_0)$. As $\frac{y_0}{y_1} = \frac{1-\frac{\rho}{N}}{R_v} + \frac{\rho}{N}$ and $\frac{\rho}{N} < 1$, $y_0 < y_1$ if $R_v > 1$. 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 \ge 0$) by the tangent line at $c_2 = 0$ or the parabola, respectively (black dotted and dashed lines):

$$ \begin{aligned} y_0 + k_1 c_2^* &= y_1, \\ y_0 + k_1 c_2^{**} + k_2 c_2^{**2} &= y_1, \end{aligned} $$

where

$$ \begin{aligned} k_1 &= [\text{LHS}]'_{c_2=0} = \gamma\beta + 2(\gamma + \beta)(\sigma + \rho au_0), \\ k_2 &= \frac{1}{2}[\text{LHS}]''_{c_2=0} = \gamma + \beta + 2\sigma + 2\rho au_0. \end{aligned} $$

Solving these linear and quadratic equations, we get two estimates for $c_{\min}^2 = 2D \cdot c_2$:

$$ c_{\min}^{*2} = 4D \frac{R_v - 1}{\frac{1}{\sigma} + 2(\frac{1}{\gamma} + \frac{1}{\beta})(1 + \frac{R_v}{\frac{N}{\rho} - 1})}, \quad (A8a) $$

$$ c_{\min}^{**2} = Dk_{12} \left[ \sqrt{1 + \frac{8}{k_{12}} \cdot \frac{R_v - 1}{\frac{1}{\sigma} + 2(\frac{1}{\gamma} + \frac{1}{\beta})(1 + \frac{R_v}{\frac{N}{\rho} - 1})}} - 1 \right], \quad (A8b) $$

where

$$ k_{12} = \frac{k_1}{k_2} = \frac{\gamma\beta + 2(\gamma + \beta)\sigma(1 + \frac{R_v}{\frac{N}{\rho} - 1})}{\gamma + \beta + 2\sigma(1 + \frac{R_v}{\frac{N}{\rho} - 1})} = \frac{1 + 2\frac{\tau_I + \tau_E}{\tau_V}(1 + \frac{R_v}{\frac{N}{\rho} - 1})}{\tau_I + \tau_E + 2\frac{\tau_I \tau_E}{\tau_V}(1 + \frac{R_v}{\frac{N}{\rho} - 1})}, $$

and

$$ \frac{R_v}{\frac{N}{\rho} - 1} \equiv \frac{\rho au_0}{\sigma}. $$
<figure><img src="image_15.png" /><figcaption>**Figure A4.** Solution to Equation (A7) (solid lines) and two its approximations in the $c_2 \ge 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).</figcaption></figure>

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 caption. In that figure, the solid line (corresponding to the minimal speed for which the autowave solution in the system (A1) exists) is plotted assuming $c = 0.38$. This value was obtained from Equation (A5): $\min F(\mu) \approx 0.144$ (Figure A3, $\rho = 1$ curve), so $c_{min} \approx \sqrt{0.144} \approx 0.38$. Estimates from Equations (A8) give $c_{min}^* \approx \sqrt{0.22} \approx 0.47$ and $c_{min}^{**} \approx \sqrt{0.21} \approx 0.46$. For the same parameters, numerical solution of Equation (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 correspond 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 $\rho = 0$ (A) and $\rho = 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 $\sqrt{\min[F(\mu)]}$ 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(\mu)$, $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><img src="image_16.png" /><figcaption>**Figure A5.** Steady autowave speed (squared) as a function of $R_v$ in the numerical solution (Equation (1), squares) and estimated using Equations (A5) (circles), (A8a) (dashed lines), and (A8b) (solid lines). (A) $\rho = 0$, (B) $\rho = 1$. Other parameters are the same as in Figure A2.</figcaption></figure>
## Appendix A.3. Derivation of Equations for the Final Concentration of Intact Cells and the Total Spatial Viral Load

Equation (8) for the final concentration of intact cells and the total spatial viral load can be derived as follows. From Equation (2a),

$$ 
\begin{aligned}
-c \frac{u'}{u} &= -av, \\
c \int_{u_0}^{u_f} \frac{du}{u} &= a \int_{+\infty}^{-\infty} v(\xi) d\xi, \\
c \ln \frac{u_f}{u_0} &= -a \int_{-\infty}^{+\infty} v(\xi) d\xi = -a\tilde{v}.
\end{aligned}
$$

From combined Equations (2a) and (2b),

$$ 
\begin{aligned}
-c(u + p)' &= -\gamma p, \\
c \int_{u(-\infty)}^{u(+\infty)} du + c \int_{p(-\infty)}^{p(+\infty)} dp &= \gamma \int_{-\infty}^{+\infty} p(\xi) d\xi, \\
c(u_0 - u_f) + 0 - 0 &= \gamma \tilde{p}.
\end{aligned}
$$

From combined Equations (2a)–(2c),

$$ 
\begin{aligned}
-c(u + p + q)' &= -\beta q, \\
c \int_{u(-\infty)}^{u(+\infty)} du + c \int_{p(-\infty)}^{p(+\infty)} dp + c \int_{q(-\infty)}^{q(+\infty)} dq &= \beta \int_{-\infty}^{+\infty} q(\xi) d\xi, \\
c(u_0 - u_f) + 0 - 0 + 0 - 0 &= \beta \tilde{q}.
\end{aligned}
$$

From Equation (2d),

$$ 
\begin{aligned}
-c \int_{v(-\infty)}^{v(+\infty)} dv &= -\rho a \int_{-\infty}^{+\infty} uv d\xi + N\beta \int_{-\infty}^{+\infty} q(\xi) d\xi - \sigma \int_{-\infty}^{+\infty} v(\xi) d\xi + D \int_{v'_{\xi}(-\infty)}^{v'_{\xi}(+\infty)} dv'_{\xi}, \\
-0 + 0 &= -\rho a \tilde{u} + N\beta \tilde{q} - \sigma \tilde{v} + 0 - 0.
\end{aligned}
$$

From combined Equations (2a) and (2d),

$$ 
\begin{aligned}
-c(v - \rho u)' &= N\beta q - \sigma v + D v'', \\
-c \int_{v(-\infty)}^{v(+\infty)} dv + c\rho \int_{u(-\infty)}^{u(+\infty)} du &= N\beta \int_{-\infty}^{+\infty} q(\xi) d\xi - \sigma \int_{-\infty}^{+\infty} v(\xi) d\xi + D \int_{v'_{\xi}(-\infty)}^{v'_{\xi}(+\infty)} dv'_{\xi}, \\
-0 + 0 + c\rho(u_0 - u_f) &= N\beta \tilde{q} - \sigma \tilde{v} + 0 - 0.
\end{aligned}
$$

Altogether, we have:

$$ 
\begin{aligned}
-c \ln \frac{u_f}{u_0} &= a\tilde{v}, \\
c(u_0 - u_f) &= \gamma \tilde{p} = \beta \tilde{q}, \\
c\rho(u_0 - u_f) &= N\beta \tilde{q} - \sigma \tilde{v} = \rho a \tilde{u} \tilde{v}.
\end{aligned}
$$

From the first equation, we obtain Equation (8b). Eliminating $\tilde{q}$ and $\tilde{v}$ and taking into account Equation (4), we obtain Equation (8a).
References

1. Hui, K.P.; Ho, J.C.; Cheung, M.c.; Ng, K.c.; Ching, R.H.; Lai, K.I.; Kam, T.T.; Gu, H.; Sit, K.Y.; Hsin, M.K.; et al. SARS-CoV-2 Omicron variant replication in human bronchus and lung ex vivo. *Nature* **2022**, *603*, 715–720. [CrossRef] [PubMed]

2. Peacock, T.P.; Brown, J.C.; Zhou, J.; Thakur, N.; Newman, J.; Kugathasan, R.; Sukhova, K.; Kaforou, M.; Bailey, D.; Barclay, W.S. The SARS-CoV-2 variant, Omicron, shows rapid replication in human primary nasal epithelial cultures and efficiently uses the endosomal route of entry. *bioRxiv* (version posted January 3) 2022, preprint. [CrossRef]

3. Murray, J. *Mathematical Biology. I. An Introduction*, 3rd ed.; Springer: New York, NY, USA, 2002. [CrossRef]

4. Murray, J. *Mathematical Biology. II. Spatial Models and Biomedical Applications*, 3rd ed.; Springer: New York, NY, USA, 2003. [CrossRef]

5. Flerlage, T.; Boyd, D.F.; Meliopoulos, V.; Thomas, P.G.; Schultz-Cherry, S. Influenza virus and SARS-CoV-2: Pathogenesis and host responses in the respiratory tract. *Nat. Rev. Microbiol.* **2021**, *19*, 425–441. [CrossRef] [PubMed]

6. V'kovski, P.; Kratzel, A.; Steiner, S.; Stalder, H.; Thiel, V. Coronavirus biology and replication: Implications for SARS-CoV-2. *Nat. Rev. Microbiol.* **2021**, *19*, 155–170. [CrossRef] [PubMed]

7. Hou, Y.J.; Okuda, K.; Edwards, C.E.; Martinez, D.R.; Asakura, T.; Dinnon, K.H.; Kato, T.; Lee, R.E.; Yount, B.L.; Mascenik, T.M.; et al. SARS-CoV-2 Reverse Genetics Reveals a Variable Infection Gradient in the Respiratory Tract. *Cell* **2020**, *182*, 429–446.e14. [CrossRef] [PubMed]

8. Gąsecka, A.; Borovac, J.A.; Guerreiro, R.A.; Giustozzi, M.; Parker, W.; Caldeira, D.; Chiva-Blanch, G. Thrombotic Complications in Patients with COVID-19: Pathophysiological Mechanisms, Diagnosis, and Treatment. *Cardiovasc. Drugs Ther.* **2021**, *35*, 215–229. [CrossRef] [PubMed]

9. WHO. Tracking SARS-CoV-2 Variants. Available online: https://www.who.int/en/activities/tracking-SARS-CoV-2-variants (accessed on 13 June 2022).

10. Dashkevich, N.M.; Ovanesov, M.V.; Balandina, A.N.; Karamzin, S.S.; Shestakov, P.I.; Soshitova, N.P.; Tokarev, A.A.; Panteleev, M.A.; Ataullakhanov, F.I. Thrombin activity propagates in space during blood coagulation as an excitation wave. *Biophys. J.* **2012**, *103*, 2233–2240. [CrossRef] [PubMed]

11. Ataullakhanov, F.; Guriia, G. Spatial aspects of the dynamics of blood coagulation. I. Hypothesis. *Biophysics* **1994**, *39*, 89–96.

12. Tokarev, A.; Ratto, N.; Volpert, V. Mathematical modeling of thrombin generation and wave propagation: from simple to complex models and backwards. In *BIOMAT 2018 International Symposium on Mathematical and Computational Biology*; Springer: Cham, Switzerland, 2019. [CrossRef]

13. Bocharov, G.; Meyerhans, A.; Bessonov, N.; Trofimchuk, S.; Volpert, V. Interplay between reaction and diffusion processes in governing the dynamics of virus infections. *J. Theor. Biol.* **2018**, *457*, 221–236. [CrossRef] [PubMed]

14. Mahiout, L.; Mozokhina, A.; Tokarev, A.; Volpert, V. Virus replication and competition in a cell culture: Application to the SARS-CoV-2 variants. *Appl. Math. Lett.* **2022**, *113*, 108217. [CrossRef] [PubMed]

15. Sender, R.; Bar-On, Y.M.; Gleizer, S.; Bernshtein, B.; Flamholz, A.; Phillips, R.; Milo, R. The total number and mass of SARS-CoV-2 virions. *Proc. Natl. Acad. Sci. USA* **2021**, *118*, 1–9. [CrossRef] [PubMed]

16. Zhu, N.; Zhang, D.; Wang, W.; Li, X.; Yang, B.; Song, J.; Zhao, X.; Huang, B.; Shi, W.; Lu, R.; et al. A Novel Coronavirus from Patients with Pneumonia in China, 2019. *N. Engl. J. Med.* **2020**, *382*, 727–733. [CrossRef] [PubMed]

17. Epithelix. FAQ, Epithelix's Tissues. Available online: https://www.epithelix.com/faq (accessed on 13 June 2022).

18. Hoops, S.; Sahle, S.; Gauges, R.; Lee, C.; Pahle, J.; Simus, N.; Singhal, M.; Xu, L.; Mendes, P.; Kummer, U.; et al. COPASI—A COMplex PAthway SImulator. *Bioinformatics* **2006**, *22*, 3067–3074. [CrossRef] [PubMed]

19. Tokarev, A.A. Velocity–Amplitude Relationship in the Gray–Scott Autowave Model in Isolated Conditions. *ACS Omega* **2019**, *4*, 14430–14438. [CrossRef] [PubMed]

20. Hairer, E.; Nørsett, S.; Wanner, G. *Solving Ordinary Differential Equations I. Nonstiff Problems* (Russian Translation, 1990); Mir: Moscow, Russia; Springer: Berlin/Heidelberg, Germany; New York, NY, USA; London, UK; Paris, France; Tokyo, Japan, 1987.