Viral Infection Spreading in Cell Culture with Intracellular Regulation

Abstract

Virus plaque assays are conventionally used for the assessment of viral infections, including their virulence, and vaccine efficacy. These experiments can be modeled with reaction–diffusion equations, allowing the estimation of the speed of infection spread (related to virus virulence) and viral load (related to virus infectivity). In this work, we develop a multiscale model of infection progression that combines macroscopic characterization of virus plaque growth in cell culture with a reference model of intracellular virus replication. We determine the infection spreading speed and viral load in a model for the extracellular dynamics and the kinetics of the abundance of intracellular viral genomes and proteins. In particular, the spatial infection spreading speed increases if the rate of virus entry into the target cell increases, while the viral load can either increase or decrease depending on other model parameters. The reduction in the model under a quasi-steady state assumption for some intracellular reactions allows us to derive a family of reduced models and to compare the reference model with the previous model for the concentration of uninfected cells, infected cells, and total virus concentration. Overall, the combination of different scales in reaction–diffusion models opens up new perspectives on virus plaque growth models and their applications.

1. Introduction

Virus plaque assays are a conventional experimental test used for the assessment of viral infections [1]. Plaque size, morphology, clearness, and distribution characterize virus progression and evolution [2]. In particular, plaque size caused by cell necrosis or apoptosis in cell culture correlates with virus virulence [3,4,5] and depends on the viral cell-to-cell transmission rate and on its ability to avoid the immune response mediated by interferon produced by infected cells [6]. Plaque assays are also used for the estimation of virus concentration in a multiplicity of infection tests under the assumption that each plaque corresponds to a single virion (plaque forming unit), though this may not be the case for some viruses [7]. Recently, plaque assays were used to assess different variants of SARS-CoV-2 infection [8].

Virus plaque growth is conditioned by the production of new virus particles in infected cells using their genetic machinery. Plaque assay tests show that the process of plaque growth in time and space depends on its spatial spreading, determined by virus diffusion in the extracellular matrix [9]. Virus plaque growth can be studied with different modeling approaches. A reaction–diffusion model of plaque growth is considered in [10] in the case of reversible host infection. The plaque growth rate is determined by the method of linearization. Numerical simulations of viral plaque growth described by a reaction–diffusion model with time delay are presented in [11]. Individual-based models of viral infection spreading in cell cultures are developed in [12,13,14]. One- and two-dimensional models of plaque growth are compared in [15]. Though the 2D model gives a slightly larger value of the viral load, the main features of plaque growth are the same.

It should be noted that theoretical models of virus plaque assays deal with the concentrations of infected cells but not with the intracellular regulation of viral replication. On the other hand, the models that address intracellular regulation do not take into account the spatial spread of the infection. Therefore, the investigation of the influence of different stages of intracellular virus replication on plaque growth requires the development of new models addressing intracellular regulation and spatial spreading at the same time. We develop such a model in this work.

A systematic investigation of viral infection spreading in cell culture was started in [15], where a reaction–diffusion system with time delay for the concentrations of uninfected cells, infected cells, and virus was used to find the spreading speed and viral load. Different variants of the SARS-CoV-2 infection were assessed in [16], where it was shown that in the competition of two viruses in cell culture, the virus with larger individual spreading speed persists and eliminates another one. These results corroborate the experimental data [8]. The progression of viral infections in the respiratory tract can also be described as a reaction–diffusion wave, but its propagation is influenced by the airway liquids (mucus and periciliary fluid) [17]. Mathematical questions of the existence of reaction–diffusion waves describing infection spreading in its interaction with the immune response are studied in [18,19].

This cycle of works describes infection spreading at the level of cell culture or tissue without taking into account the intracellular regulation of virus replication. This macroscopic approach is sufficiently simple that it allows the description of the experimental data, and there are only a few parameters that can be readily estimated. This approach can be combined with a more detailed description of intracellular regulation using the multiscale modeling method. In the case of immovable cells with fixed positions, we can pass from cells to the intracellular concentrations, considered as functions of space and of time. Thus, instead of the concentrations of cells, we consider concentrations of intracellular substances, implicitly associating them with each other. Intracellular proteins do not diffuse if we neglect direct cell-to-cell transport through cell junctions. On the other hand, extracellular factors (e.g., interferon) and viruses diffuse and propagate in the culture, resulting in infection spreading.

In this work, we study infection spreading in cell culture, taking into account a simplified description of intracellular virus replication. The main stages of this process are shown in Figure 1, starting with the extracellular virions entering the cell and releasing into the cytoplasm its proteins and genetic material. The latter participate in virus replication, together with some intracellular molecules. New virions are released back into the extracellular matrix. All steps of intracellular virus replication, including translation, transcription, and assembly are described here as a single reaction process for viral genomes R and required proteins A

$$R+A\to R+V_r.$$

Here, the depletion of R is neglected, since a single intracellular viral genome induces the production of multiple new viral genomes.

Figure 1. Basic model considered in this work (a) with extracellular virus $V_e$, intracellular virus $V_i$, viral genetic material R, and intracellular proteins A participating in virus replication and production of new virus particles $V_r$. In the reduced model (b), two intermediate stages are omitted.

Figure 1. Basic model considered in this work (a) with extracellular virus $V_e$, intracellular virus $V_i$, viral genetic material R, and intracellular proteins A participating in virus replication and production of new virus particles $V_r$. In the reduced model (b), two intermediate stages are omitted.

The equation for the extracellular virus concentration $V_e$,

$$\frac{\partial V_e}{\partial t}=D\frac{{\partial }^2{V}_e}{\partial x^2}-q_1{V}_e+q_2{V}_r-{\sigma }_1{V}_e$$

(1)

takes into account extracellular virus diffusion. Its concentration decreases due to its entry from the extracellular space into cells ($q_1{V}_e$) and increases due to new virus virions produced in the infected cells ($q_2{V}_r)$. The last term on the right-hand side of this equation describes virus death. The intracellular virus concentration $V_i$ is determined by its entry from the extracellular space followed by its uncoating:

$$\frac{\partial V_i}{\partial t}=q_1{V}_e-k_1{V}_i.$$

(2)

Viral genetic material (DNA or RNA) is released after uncoating, and it can also be eliminated due to various protective mechanisms:

$$\frac{\partial R}{\partial t}=k_1{V}_i-{\sigma }_{2}R.$$

(3)

Virus replication is considered a single-stage process determined by the concentrations of viral genetic material R and intracellular proteins A:

$$\frac{\partial V_r}{\partial t}=k_{3}AR-q_2{V}_r.$$

(4)

The last term in this equation characterizes the degradation or elimination of new viral particles. Finally, intracellular concentrations A participating in virus replication are described by the equation

$$\frac{\partial A}{\partial t}=-k_{4}AR.$$

(5)

Let us note that constants $k_3$ and $k_4$ can be different since multiple intracellular molecules are used in the production of a new virus particle. All concentrations in Equations (1)–(5) depend on the space variable x and on time t. This system of equations will be considered on the whole axis in the theoretical analysis and in a bounded interval in numerical simulations with Neuman (no-flux) boundary conditions and some non-negative initial conditions. It can be verified by conventional methods that the solution to this problem is bounded and non-negative.

More detailed models of intracellular regulation with some intermediate steps in virus replication can be introduced and considered under quasi-stationary approximation. This approach justifies models (1)–(5), allowing the analysis presented below.

In this study, a family of related mathematical models is presented to describe and analyze the 1D diffusion-mediated spatio-temporal kinetics of viral infection spreading in cell culture. The models differ in the details of the description of the intracellular viral life cycle. In Section 2, we derive analytical expressions for the basic viral reproduction number, the virus propagation wave speed, and the total viral load. Numerical simulations are presented that specify the quantitative dependencies of these basic characteristics on the model parameters. In Section 3, two models of reduced complexity compared to the reference one are derived and analyzed. Validation of the kinetics of virus replication in infected cells is performed in Section 4. The study concludes with a Discussion in Section 5.

2. Viral Load and Spreading Speed

2.1. Virus Replication Number

Consider the reaction part of the system (1)–(5), that is, the model without diffusion in the first equation:

$$\frac{d{V}_e}{dt}=q_2{V}_r-q_1{V}_e-{\sigma }_1{V}_e,$$

(6)

$$\frac{d{V}_i}{dt}=q_1{V}_e-k_1{V}_i,$$

(7)

$$\frac{dR}{dt}=k_1{V}_i-{\sigma }_{2}R,$$

(8)

$$\frac{d{V}_r}{dt}=k_{3}AR-q_2{V}_r,$$

(9)

$$\frac{dA}{dt}=-k_{4}AR.$$

(10)

The solution to this system with a positive initial condition is bounded and positive. Stationary points of this system are $V_r=0$, $V_e=0$, $V_i=0$, $R=0$, and any A.

The system linearized about stationary point ($A=A_0,0,0,0,0)$ is similar to the previous one where A in Equations (9) and (10) is replaced by $A_0$. We consider the corresponding matrix:

$$\mathcal{A}=\left(\begin{array}{cccc}-q_1-{\sigma }_1& 0& 0& q_2\\ q_1& -k_1& 0& 0\\ 0& k_1& -{\sigma }_2& 0\\ 0& 0& k_3{A}_0& -q_2\end{array}\right).$$

Its characteristic polynomial has positive coefficients. According to Descartes’ rule of signs, it has a positive root if and only if the last coefficient is negative:

$$\begin{array}{c}det\left(\mathcal{A}\right)=(q_1+{\sigma }_1)k_1{\sigma }_2{q}_2-q_2{q}_1{k}_1{k}_3{A}_0<0.\end{array}$$

(11)

We introduce the viral replication number:

$$R_v=\frac{q_1{k}_3{A}_0}{{\sigma }_2(q_1+{\sigma }_1)}.$$

This number represents the ratio of virus replication and virus elimination rates [16,17,20]. If $R_v>1$, then inequality (11) holds and has a positive eigenvalue, such that its solution grows with time. The virus replication number represents the ratio of virus production and elimination rates. Virus concentration grows if the virus replication number is larger than 1.

2.2. Viral Load in the Case of Uniform Virus Distribution

Consider the previous model (6)–(10) under the assumptions

$$A\left(0\right)=A_0,V_e\left(0\right)=V_i\left(0\right)=V_r\left(0\right)=R\left(0\right)=0,$$

(12)

$$A(\infty )=A_f,V_e(\infty )=V_i(\infty )=V_r(\infty )=R(\infty )=0.$$

(13)

Integrating Equations (6) and (7) from 0 to ∞, we obtain

$$q_{2}I\left(V_r\right)=(q_1+{\sigma }_1)I\left(V_e\right),q_{1}I\left(V_e\right)=k_{1}I\left(V_i\right),$$

(14)

where

$$I\left(V_r\right)={\int }_0^{\infty }{V}_r\left(t\right)dt,I\left(V_e\right)={\int }_0^{\infty }{V}_e\left(t\right)dt.$$

Taking a linear combination of Equations (8) and (10) and integrating, we have

$$-\frac{k_2}{k_4}\left(A_f-A_0\right)=k_{1}I\left(V_i\right)-{\sigma }_{2}I\left(R\right).$$

(15)

From Equations (9) and (10)

$$\frac{k_3}{k_4}(A_f-A_0)=-q_{2}I\left(V_r\right).$$

(16)

Finally, dividing Equation (10) and integrating, we obtain

$$ln\left(\frac{A_f}{A_0}\right)=-k_{3}I\left(R\right).$$

(17)

Expressing the integrals from equalities (14)–(17), we obtain the following equation

$$R_v(\omega -1)=ln\omega$$

(18)

with respect to $\omega =A_f/A_0$. Equation (18) has a solution $\omega \in (0,1)$ if and only if $R_v>1$. This case corresponds to infection progression. The total quantity of replicated virus, $I\left(V_r\right)$, can now be determined from Equation (16):

$$I\left(V_r\right)=\frac{k_3{A}_0}{k_4{q}_2}(1-\omega )\approx \frac{k_3{A}_0}{k_4{q}_2}.$$

(19)

The last approximation is applicable for large $R_v$ since $\omega \ll 1$. The total quantity of extracellular virus is given by the following expression:

$$I\left(V_e\right)\approx \frac{k_3{A}_0}{k_4(q_1+{\sigma }_1)}.$$

(20)

2.3. Viral Load in the Wave

Consider now system (1)–(5) with diffusion. As is known from the previous works (e.g., in [15,18]), infection spreads in cell culture as a reaction–diffusion wave (Figure 2).

Consider the travelling wave solution $V_e(x,t)=v(x-ct)$, $R(x,t)=u(x-ct)$, $A(x,t)=w(x-ct)$, $V_i(x,t)=y(x-ct)$, $V_r(x,t)=z(x-ct)$ with the limits at infinity:

$$v(\pm \infty )=u(\pm \infty )=y(\pm \infty )=z(\pm \infty )=0,w(-\infty )=w_f,w(+\infty )=w_0,$$

(21)

where $w_f$ is the unknown final concentration of the intracellular components. Thus, the system takes the form:

$$D{v}^{″}+c{v}^{\prime }+q_{2}z-q_{1}v-{\sigma }_{1}v=0,$$

(22)

$$c{y}^{\prime }+q_{1}v-k_{1}y=0,$$

(23)

$$c{u}^{\prime }+k_{1}y-{\sigma }_{2}u=0,$$

(24)

$$c{z}^{\prime }+k_{3}wu-q_{2}z=0,$$

(25)

$$c{w}^{\prime }-k_{4}wu=0,$$

(26)

where prime designates the derivative on $\xi =x-ct$. Integrating Equations (22)–(25) over $\xi \in (-\infty ,+\infty )$, separating variables, and integrating the last Equation (26), and using the boundary conditions, we obtain the following system of equations:

$$q_{2}J\left(z\right)=(q_1+{\sigma }_1)J\left(v\right),$$

(27)

$$q_{1}J\left(v\right)=k_{1}J\left(y\right),$$

(28)

$$k_{1}J\left(y\right)={\sigma }_{2}J\left(u\right),$$

(29)

$$cln\frac{w_f}{w_0}=-k_{4}J\left(u\right),$$

(30)

$$\frac{c{k}_3}{k_4}(w_0-w_f)=q_{2}J\left(z\right).$$

(31)

Here, we use notation $J\left(f\right)={\int }_{-\infty }^{+\infty }fd\xi$. Excluding the integrals from this system, we obtain the equation

$$R_v(\tilde{w}-1)=ln\tilde{w},$$

(32)

where $\tilde{w}=w_f/w_0$, and $R_v=q_1{k}_3{w}_0/\left({\sigma }_2(q_1+{\sigma }_1)\right)$ is the same as (11) in reaction system, where $A_0$ is substituted by $w_0$. This equation has a non-zero solution in $\tilde{w}\in (0,1)$ for $R_v>1$. In this case, the viral load is determined by the equality

$$J\left(v\right)={\int }_{-\infty }^{+\infty }vd\xi =-\frac{c{\sigma }_2}{q_1{k}_4}ln\tilde{w}$$

(33)

or

$$J\left(v\right)=\frac{c{w}_0}{q_1+{\sigma }_1}\frac{k_3}{k_4}(1-\tilde{w}).$$

(34)

If $\tilde{w}\ll 0$, then $ln\tilde{w}\approx -R_v$ and

$$J\left(v\right)\approx \frac{c{w}_0{k}_3}{(q_1+{\sigma }_1)k_4},$$

(35)

where c is the wave speed, which will be determined below.

2.4. Wave Speed

To determine the wave speed, we use the linearization method, in which we replace the value of w in the system (22)–(26) by its value $w_0$ at infinity:

$$D{v}^{″}+c{v}^{\prime }+q_{2}z-q_{1}v-{\sigma }_{1}v=0,$$

(36)

$$c{y}^{\prime }+q_{1}v-k_{1}y=0,$$

(37)

$$c{u}^{\prime }+k_{1}y-{\sigma }_{2}u=0,$$

(38)

$$c{z}^{\prime }+k_3{w}_{0}u-q_{2}z=0.$$

(39)

We search for the solution in the form $v\left(\xi \right)=p{e}^{-\lambda \xi }$, $u\left(\xi \right)=q{e}^{-\lambda \xi }$, $y\left(\xi \right)=m{e}^{-\lambda \xi }$, $z\left(\xi \right)=l{e}^{-\lambda \xi }$. Substituting this solution into the system (36)–(39), we obtain the following characteristic equation:

$$\left|\begin{array}{cccc}D{\lambda }^2-c\lambda -(q_1+{\sigma }_1)& 0& 0& q_2\\ q_1& -c\lambda -k_1& 0& 0\\ 0& k_1& -c\lambda -{\sigma }_2& 0\\ 0& 0& k_3{w}_0& -c\lambda -q_2\end{array}\right|=0,$$

(40)

or

$$(D{\lambda }^2-c\lambda -q_1-{\sigma }_1)(c\lambda +k_1)(c\lambda +{\sigma }_2)(c\lambda +q_2)=-q_1{q}_2{k}_1{k}_3{w}_0.$$

(41)

We denote $\mu =c\lambda$ and express $c^2$, and thus we obtain

$$c^2=\underset{\mu >{\mu }_0}{min}\frac{D{\mu }^2(\mu +k_1)(\mu +{\sigma }_2)(\mu +q_2)}{(\mu +q_1+{\sigma }_1)(\mu +k_1)(\mu +{\sigma }_2)(\mu +q_2)-q_1{q}_2{k}_1{k}_3{w}_0},$$

(42)

where ${\mu }_0>0$ is the value for which the denominator in the right-hand side of the last formula vanishes.

2.5. Numerical Simulations

The system of Equations (1)–(5) has been solved numerically using a finite-difference method with the first order of time approximation and the second order of space approximation. Numerical code is implemented in the C++ language with the MS Visual Studio translator. The CPU time required to calculate one time step on a mesh containing 10,000 nodes was about 1 millisecond for system (1)–(5) and about twice less for system (46)–(48).

An example of numerical simulation is shown in Figure 2 with concentration distributions in space (left) and the total viral load as a function of time (right). Let us note that the wave speed equals $0.059$ cm/h in numerical simulation and $0.058$ cm/h by Formula (42). The total viral load calculated with Formula (35) equals 5829 copy/cm${}^2$, which is in good agreement with the numerical results shown in Figure 2b.

The dependence of the wave speed and viral load on model parameters are shown in Figure 3 and Figure 4. The total viral load increases with the increase in parameters $k_1$, $q_2$, and $k_3$. The total viral load inversely depends on the parameters ${\sigma }_1$, ${\sigma }_2$, and $k_4$. The abrupt change in the dependence on ${\sigma }_2$ is determined by the transition to another regime: when ${\sigma }_2$ is big enough, $R_v$ becomes less than 1, and the infection does not develop, i.e., $J\left(V_e\right)=0$.

The total viral load can grow or decay depending on the parameter $q_1$. Figure 3 shows this dependence for two sets of parameters. The model predicts that the wave speed increases with the values of parameters $k_1$, $k_3$ and decreases for the parameters ${\sigma }_1$, ${\sigma }_2$ similarly to the total viral load (Figure 3a). The parameter $q_2$ has almost no influence on the wave speed or $k_4$. The second set of parameters does not qualitatively alter these dependencies (Figure 3b), but the curves change their relative positions, and this results in the observed qualitative changes in dependencies for the total viral load.

Table 1. Parameter estimation for the model (1)–(5). Ranges for the parameters $k_3$ and $k_4$ are calculated using respective ranges for parameters from ([21, Table 2, p. 9]) and Formula (84).

Parameter	Corresponding Parameter from [21]	Value from ([21], Table 2)
${\sigma }_1$	$d_V$	0.12 h${}^{-1}$, range: (0.06, 3.5), tuned to (0.06, 0.2)
$q_1$	$k_{fuse}$	0.5 h${}^{-1}$, range: (0.33, 1)
$q_2$	$k_{release}$	8 h${}^{-1}$ range: (8, 7200)
$k_1$	$k_{uncoat}$	0.5 h${}^{-1}$, range: (0.33, 1)
${\sigma }_2$	$d_{gRNA}$	0.2 h${}^{-1}$, range: (0.069, 0.69), tuned to (0.069, 0.4)
$k_3$	Equation (84)	$0.003$, range: ($1.05\times {10}^{-10}$, $3.7\times {10}^6$),
		range using tuned values: ($8.6\times {10}^{-9}$, $7.5\times {10}^5$)
$k_4$	Equation (84)	$6.093$, range: ($2.35\times {10}^{-7}$, $4.2\times {10}^9$),
		range using tuned values: ($1.9\times {10}^{-5}$, $8.4\times {10}^8$)

Table 1. Parameter estimation for the model (1)–(5). Ranges for the parameters $k_3$ and $k_4$ are calculated using respective ranges for parameters from ([21, Table 2, p. 9]) and Formula (84).

Parameter	Corresponding Parameter from [21]	Value from ([21], Table 2)
${\sigma }_1$	$d_V$	0.12 h${}^{-1}$, range: (0.06, 3.5), tuned to (0.06, 0.2)
$q_1$	$k_{fuse}$	0.5 h${}^{-1}$, range: (0.33, 1)
$q_2$	$k_{release}$	8 h${}^{-1}$ range: (8, 7200)
$k_1$	$k_{uncoat}$	0.5 h${}^{-1}$, range: (0.33, 1)
${\sigma }_2$	$d_{gRNA}$	0.2 h${}^{-1}$, range: (0.069, 0.69), tuned to (0.069, 0.4)
$k_3$	Equation (84)	$0.003$, range: ($1.05\times {10}^{-10}$, $3.7\times {10}^6$),
		range using tuned values: ($8.6\times {10}^{-9}$, $7.5\times {10}^5$)
$k_4$	Equation (84)	$6.093$, range: ($2.35\times {10}^{-7}$, $4.2\times {10}^9$),
		range using tuned values: ($1.9\times {10}^{-5}$, $8.4\times {10}^8$)

3. Simplified Models

3.1. Model Reduction

In this section, we study a simplified model obtained by reducing system (1)–(5) to a model of three equations. Suppose that parameters $q_1$, $k_1$, $k_3$, and $q_2$ are big enough, and $q_1=q_{10}\times N$, $k_1=k_{10}\times N$, $k_3=k_{30}\times N$, $q_2=q_{20}\times N$, where N is a large parameter. Then, the second equation has the form:

$$\frac{1}{N}\frac{\partial V_i}{\partial t}=q_{10}{V}_e-k_{10}{V}_i.$$

(43)

If $N\to \infty$, we formally obtain:

$$V_i=\frac{q_{10}}{k_{10}}{V}_e.$$

(44)

Similarly, from Equation (4),

$$V_r=\frac{k_{30}}{q_{20}}AR.$$

(45)

Substituting these expressions into (1)–(5), we obtain the following system of three equations:

$$\frac{\partial V}{\partial t}=D\frac{{\partial }^{2}V}{\partial x^2}+k_{1}AR-k_{2}V-{\sigma }_{1}V,$$

(46)

$$\frac{\partial R}{\partial t}=k_{2}V-{\sigma }_{2}R,$$

(47)

$$\frac{\partial A}{\partial t}=-k_{1}AR,$$

(48)

where we use, for convenience, the following notation: $V_e\to V,k_3=k_4\to k_1,q_1\to k_2$. We obtained the system containing three variables: extracellular virus V, intracellular viral proteins R, and cell proteins A participating in the production of new virions (Figure 1b).

3.2. Virus Replication Number

The corresponding ODE system has the form:

$$\frac{dV}{dt}=k_{1}AR-k_{2}V-{\sigma }_{1}V,$$

(49)

$$\frac{dR}{dt}=k_{2}V-{\sigma }_{2}R,$$

(50)

$$\frac{dA}{dt}=-k_{1}AR.$$

(51)

The stationary point of this system of equations can be easily determined. From Equation (51), we conclude that $AR=0$. Therefore, from Equation (49), $V=0$, and from Equation (50), $R=0$.

Consider a linear approximation where $A=A_0>0$ in Equation (49). Then, linear system of Equations (49) and (50) has a positive eigenvalue if and only if $R_v>1$, where

$$R_v=\frac{k_1{k}_2{A}_0}{(k_2+{\sigma }_1){\sigma }_2}$$

is the virus replication number. If this condition is satisfied, then solutions V and R of this system grow over time.

3.3. Viral Load

Let us complete system (49)–(51) with the initial conditions

$$A\left(0\right)=A_0,V\left(0\right)=V_0,R\left(0\right)=0.$$

(52)

It can be easily verified that the solution satisfies $V\left(t\right)\to 0$, $R\left(t\right)\to 0$, $A\left(t\right)\to A_f$ as $t\to \infty$, where $A_f\ge 0$ is some constant. Indeed, since $V^{\prime }+R^{\prime }+A^{\prime }=-{\sigma }_{2}R$, and all variables are non-negative for $t\ge 0$, then $R\left(t\right)\to 0$. We conclude from Equation (49) that $V\left(t\right)\to 0$. Since $A^{\prime }\left(t\right)\le 0$, then $A\left(t\right)$ has a limit as $t\to \infty$.

Taking a sum of Equations (49) and (51) and integrating from 0 to ∞, we obtain

$$A_f-A_0-V_0=-(k_2+{\sigma }_1)I\left(V\right),$$

(53)

where $I\left(V\right)={\int }_0^{\infty }V\left(t\right)dt$. Integrating Equation (50), we conclude that

$$k_{2}I\left(V\right)={\sigma }_{2}I\left(R\right).$$

(54)

Finally, dividing Equation (51) and integrating, we obtain

$$ln\left(\frac{A_f}{A_0}\right)=-k_{1}I\left(R\right).$$

(55)

Excluding the integrals from equalities (53)–(55), we obtain the equation

$$ln\omega =R_v(\omega -1-{\omega }_0)$$

with respect to $\omega =A_f/A_0$, where ${\omega }_0=V_0/A_0$. Assuming that ${\omega }_0\ll 1$, we have the final form of this equation

$$ln\omega =R_v(\omega -1).$$

(56)

It has a unique solution $\omega \in (0,1)$ if and only if $R_v>1$. If this condition is satisfied, then there is an outbreak of virus production. Otherwise, if $R_v<1$, its concentration converges to 0 and $\omega =1$.

The total viral load $V_T=I\left(V\right)$ can now be found as follows:

$$V_T=\frac{A_0}{k_2+{\sigma }_1}(1-\omega )\approx \frac{A_0}{k_2+{\sigma }_1}.$$

(57)

The last approximation holds for $\omega \ll 1$, which holds for sufficiently large $R_v$.

3.4. Spatial Infection Spreading

We now consider spatial distribution of viral infection in a cell culture described by the system (46)–(48). We consider this system on the whole real axis and look for a traveling wave solution, $V(x,t)=v(x-ct),R(x,t)=u(x-ct),A(x,t)=w(x-ct)$, where c is the wave speed. Substituting these functions in system (46)–(48), we obtain the following system:

$$D{v}^{″}+c{v}^{\prime }+k_{1}uw-k_{2}v-{\sigma }_{1}v=0,$$

(58)

$$c{u}^{\prime }+k_{2}v-{\sigma }_{2}u=0,$$

(59)

$$c{w}^{\prime }-k_{1}uw=0.$$

(60)

This system is considered on the whole axis with the limits at infinity:

$$v(\pm \infty )=u(\pm \infty )=0,w(-\infty )=w_f,w(\infty )=w_0.$$

(61)

Here, $w_0$ is a given constant, and $w_f$ is unknown. We will determine it below.

3.4.1. Final Amounts of Cell Proteins and Viral Load

We proceed as before for the ODE system. Taking a sum of Equations (58) and (60) and integrating over the whole axis, we obtain

$$c(w_0-w_f)=(k_2+{\sigma }_1)J\left(v\right),$$

(62)

where $J\left(v\right)={\int }_{-\infty }^{\infty }v\left(x\right)dx$. Next, from Equation (59),

$$k_{2}J\left(v\right)={\sigma }_{2}J\left(u\right).$$

(63)

Finally, from Equation (60),

$$cln\left(\frac{w_0}{w_f}\right)=k_{1}J\left(u\right).$$

(64)

Excluding the integrals from Equations (62)–(64), we obtain the equation

$$ln\tilde{\omega }=R_v(\tilde{\omega }-1)$$

(65)

with respect to $\tilde{\omega }=w_f/w_0$, where $R_v=k_1{k}_2{w}_0/\left((k_2+{\sigma }_1){\sigma }_2\right)$ is the same parameters as before with $A_0$ replaced by $w_0$. The viral load in the wave

$$J\left(v\right)=c\frac{w_0}{k_2+{\sigma }_1}$$

depends on the wave speed. We will determine it below.

3.4.2. Wave Speed

In order to determine the wave speed in problem (58)–(61), we use the linearization method where the function $w\left(\xi \right)$ is replaced by its value at infinity. Then, we obtain the system of two linear equations:

$$D{v}^{″}+c{v}^{\prime }+k_{1}u{w}_0-k_{2}v-{\sigma }_{1}v=0,$$

(66)

$$c{u}^{\prime }+k_{2}v-{\sigma }_{2}u=0.$$

(67)

The minimal wave speed should be determined as the minimal value of c for which this system has a positive solution decaying at infinity. We look for its solution in the form $v\left(\xi \right)=p{e}^{-\lambda \xi }$, $u\left(\xi \right)=q{e}^{-\lambda \xi }$. Substituting these functions in Equations (66) and (67), we obtain the characteristic equation

$$(D{\lambda }^2-c\lambda -k_2-{\sigma }_1)(c\lambda +{\sigma }_2)+k_1{k}_2{w}_0=0.$$

Denote $\mu =\lambda c$. Then, from the last equation, we obtain

$$c^2=F\left(\mu \right)\equiv \frac{D{\mu }^2(\mu +{\sigma }_2)}{(\mu +{\sigma }_2)(\mu +k_2+{\sigma }_1)-k_1{k}_2{w}_0}.$$

The minimal wave speed is given by the equality

$$c_0^2=\underset{\mu >{\mu }_0}{min}F\left(\mu \right),$$

(68)

where ${\mu }_0>0$ is the value for which the denominator of the function $F\left(\mu \right)$ vanishes.

We simulate the full system of five equations with increased values of constants $q_1$, $k_1$, $k_3$, and $q_2$ ($q_1=10$, $k_1=10$, $k_3=k_4=10$, $q_2=10$, where units and other parameters are the same as in Figure 2). Such an increase allows the model to be reduced and allows a comparison of the results for the full and reduced systems. In the simulations, the total viral load equals $0.04$ copy/cm${}^2$, and the wave speed equals $0.09$ cm/h. Using the Formulas (35) and (42), we obtain that the total viral load is $0.039$ copy/cm${}^2$ and the wave speed is $0.096$ cm/h, while Formula (68) (with $k_1=10$, $k_2=10$, ${\sigma }_1=0.01$, ${\sigma }_2=0.36$, $w_0=4$, $D=0.001$) gives wave speed $0.18$ cm/h.

3.5. Second Simplified Model

Consider system (1)–(5) and reduce it to the three-equation model under the assumption that another set of reactions is fast. As before, from Equation (3) we obtain:

$$R=\frac{k_{10}}{{\sigma }_{20}}{V}_i=\frac{q_{10}}{{\sigma }_{20}}{V}_e,$$

(69)

where ${\sigma }_2={\sigma }_{20}\times N$. Substitute this expression into the remaining equations:

$$\frac{\partial A}{\partial t}=-\frac{k_4{q}_{10}}{{\sigma }_{20}}A{V}_e,$$

(70)

$$\frac{\partial V_r}{\partial t}=\frac{k_3{q}_{10}}{{\sigma }_{20}}A{V}_e-q_2{V}_r,$$

(71)

$$\frac{\partial V_e}{\partial t}=D\frac{{\partial }^2{V}_e}{\partial x^2}+q_2{V}_r-(q_1+{\sigma }_1)V_e.$$

(72)

We introduce notation

$$U={\kappa }_{1}A,I={\kappa }_2{V}_r,V=V_e.$$

(73)

$$a\sim \frac{k_4{q}_1}{{\sigma }_2}=\frac{k_3{q}_1{\kappa }_2}{{\sigma }_2{\kappa }_1},\beta \sim q_2,b\sim \frac{q_2}{{\kappa }_2},\sigma \sim (q_1+{\sigma }_1),U_0\sim {\kappa }_1{A}_0,$$

(74)

and obtain the system:

$$\frac{\partial U}{\partial t}=-aUV,\frac{\partial I}{\partial t}=aUV-\beta I,\frac{\partial V}{\partial t}=D\frac{{\partial }^{2}V}{\partial x^2}+bI-\sigma V.$$

(75)

This system of equations for the concentrations of uninfected cells U, infected cells I, and virus V was considered previously in [15], taking into account the time delay. In the current model, A corresponds to the concentration of cell proteins participating in virus replication (see, e.g., [22]), $V_r$ corresponds to the concentration of viral copies inside the infected cells after the replication, and $V_e$ corresponds to virus concentration in the extracellular space. Therefore, we can reduce the multiscale model with intracellular regulation to the macroscopic model for cell concentrations. Their respective spreading speeds and viral loads are in agreement with each other.

We find the total viral load in the five-equation system with increased $q_1$, $k_1$, and ${\sigma }_2$, which correspond to this reduced system. In the simulations, the viral load is 106 copy/cm${}^2$ and the wave speed is $0.116$ cm/h. The viral load and the wave speed obtained with (35) and (42) are 107 copy/cm${}^2$ and $0.118$ cm/h, respectively, and the viral load and the wave speed obtained by corresponding formulas from [15] are 221 copy/cm${}^2$ and $0.24$ cm/h, respectively. Further increases in coefficients $q_1$, $k_1$, and ${\sigma }_2$ results in a better correspondence between full and reduced models (

Table 2. The results of comparison of the total viral load and the wave speed obtained in the numerical simulation of the full system with the results obtained by analytical formulas from [15] for the second reduced model. The increase in coefficients $q_1$, $k_1$, and ${\sigma }_2$ leads to better correspondence. Parameters $q_1$, $k_1$, and ${\sigma }_2$ are measured in h${}^{-1}$, $J\left(V_e\right)$ is measured in copy/cm${}^2$, and c is in cm/h. Values of other parameters are $D=0.01$, ${\sigma }_1=0.1$, $k_3=0.1$, $k_4=0.1$, $q_2=0.1$, $A_0=1000$, with the same units as in Figure 2. Corresponding values of the parameters of system (75) are $a=0.1$ (h × copy)${}^{-1}$, $b=0.1$ copy/(h × cell), $\beta =0.1$ h${}^{-1}$, $U_0=1000$ cell/mL, and ${\kappa }_1={\kappa }_2=1.0$.

Value of ${\mathit{q}}_1={\mathit{k}}_1={\mathit{\sigma }}_2$	Numerical $\mathit{J}\left({\mathit{V}}_{\mathit{e}}\right)$	Numerical c	$\mathit{J}\left({\mathit{V}}_{\mathit{e}}\right)$ by [15]	c by [15]
1	106	0.116	221	0.24
2	54	0.11	99	0.21
5	17	0.09	25	0.13
10	5.1	0.05	6.2	0.06

).

4. Model Validation

We showed in the previous section that model (1)–(5) can be reduced to the model of three Equations (75) for the concentrations of uninfected cells, infected cells, and virus. The latter is used to describe the experimental results for various virus types, including Delta and Omicron variants of the SARS-CoV-2 infection [15,16]. Therefore, the formulated multiscale model can also be used to describe these data.

In this section, we validate the model by comparing it with kinetics of the SARS-CoV-2 life cycle presented in [21], where a detailed model of theintracellular reactions is considered, and the parameters for this model are estimated using the appropriate virological and molecular biology data. As this model contains more intermediate stages, and it does not take into account the space variable, a direct comparison is impossible. Therefore, we adapt the model and consider system (1)–(5) under the homogeneous condition, with an additional assumption that $q_2=0$, since original and replicated virions are considered in [21] separately. Thus, we obtain the following system

$$\frac{\mathrm{d}{V}_e}{\mathrm{d}t}=-(q_1+{\sigma }_1)V_e,\frac{\mathrm{d}{V}_i}{\mathrm{d}t}=q_1{V}_e-k_1{V}_i,$$

(76)

$$\frac{\mathrm{d}R}{\mathrm{d}t}=k_1{V}_i-{\sigma }_{2}R,\frac{\mathrm{d}{V}_r}{\mathrm{d}t}=k_{3}AR,\frac{\mathrm{d}A}{\mathrm{d}t}=-k_{4}AR.$$

(77)

System (76), (77) can be solved directly:

$$V_e=V_e\left(0\right)\times e^{-(q_1+{\sigma }_1)t},$$

$$V_i=C_1{e}^{-k_{1}t}+B_1{e}^{-(q_1+{\sigma }_1)t},B_1=\frac{q_1{V}_e\left(0\right)}{k_1-(q_1+{\sigma }_1)},C_1=V_i\left(0\right)-B_1,$$

$$R=C_2{e}^{-{\sigma }_{2}t}+B_2{e}^{-k_{1}t}+B_3{e}^{-(q_1+{\sigma }_1)t},$$

$$B_2=\frac{k_1{C}_1}{{\sigma }_2-k_1},B_3=\frac{k_1{B}_1}{{\sigma }_2-(q_1+{\sigma }_1)},C_2=R\left(0\right)-B_2-B_3,$$

$$A=A\left(0\right)\times exp\left(\frac{k_4{C}_2}{{\sigma }_2}{e}^{-{\sigma }_{2}t}+\frac{k_4{B}_2}{k_1}{e}^{-k_{1}t}+\frac{k_4{B}_3}{q_1+{\sigma }_1}{e}^{-(q1+{\sigma }_1)t}\right),$$

$$V_r=V_r\left(0\right)+k_3{\int }_0^{t}AR\mathrm{d}t,$$

(78)

where the last integral can be found in elementary functions.

On the other hand, we reduce the system from [21] to the following five equations (see Appendix A) and obtain the following model:

$$\frac{\mathrm{d}\left[V_{free}\right]}{\mathrm{d}t}=k_{release}\left[V_{assembled}\right]-\left(k_{fuse}+d_V\right)\left[V_{free}\right],$$

(79)

$$\frac{\mathrm{d}\left[V_{endosome}\right]}{\mathrm{d}t}=\frac{k_{fuse}{k}_{bind}}{k_{diss}}\left[V_{free}\right]-\left(k_{uncoat}+d_{endosome}\right)\left[V_{endosome}\right],$$

(80)

$$\frac{\mathrm{d}\left[gRN{A}_{(+)}\right]}{\mathrm{d}t}=k_{uncoat}\left[V_{endosome}\right]-d_{gRNA}\left[gRN{A}_{(+)}\right],$$

(81)

$$\frac{\mathrm{d}\left[V_{assembled}\right]}{\mathrm{d}t}=k_3\left[SP\right]{\left[gRN{A}_{(+)}\right]}^3-(k_{release}+d_{assembled})\left[V_{assembled}\right],$$

(82)

$$\frac{\mathrm{d}\left[SP\right]}{\mathrm{d}t}=-k_4\left[SP\right]{\left[gRN{A}_{(+)}\right]}^3,$$

(83)

where

$$k_3=\frac{k_{assemb}}{k_{V_{rel}}{n}_{SP}}W,k_4=\frac{k_{assemb}}{K_{V_{rel}}}W,$$

(84)

$$W=\frac{k_{transl}{f}_N}{d_{N-gRNA}{n}_N}\frac{k_{t{r}_{(+)}}{k}_{t{r}_{(-)}}}{d_{gRNA}{K}_{NSP}^2{d}_{gRN{A}_{(-)}}}{\left(\frac{k_{transl}{f}_{ORF1}}{d_{NSP}}\right)}^2.$$

(85)

Comparing this model with the system (1)–(5), we obtain following correspondences for variables:

$$V_e\sim \left[V_{free}\right],V_i\sim \left[V_{endosome}\right],R\sim \left[gRN{A}_{(+)}\right],V_r\sim \left[V_{assembled}\right],A\sim \left[SP\right]$$

and parameters:

$$q_1\sim k_{fuse},{\sigma }_1\sim d_V,k_1\sim (k_{uncoat}+d_{endosome})\approx k_{uncoat},{\sigma }_2\sim d_{gRNA},q_2\sim k_{release}.$$

The values of parameters and their ranges for the model (1)–(5) are listed in Table 1. The results of the calculation of this model and its comparison with the results from ([21], Figure 2) are presented in Figure 5. The curves are similar, but there is no time delay.

After the comparison with results from [21], we obtain the ranges of parameters for the five-equation model (Table 1), and we establish the nature of the dependence of the viral load on the parameters in this model. These dependencies are shown in Figure 3 and Figure 4.

Let us also note that the simplified model (75) is similar to the model studied previously in [15], where the values of parameters were determined by the comparison with the experimental data. This model reduction provides an additional validation of the parameter choice.

Figure 5. Comparison of the kinetics of SARS-CoV-2 replication predicted by model (1)–(5) (right) with results from ([21], Figure 2) (left). The values of parameters are estimated from the model reduction (Table 1): $k_1=0.5$, $q_1=0.5$, ${\sigma }_1=0.12$, ${\sigma }_2=0.2$, $k_3=0.003$, $q_2=0$, $k_4=6.093$, and the initial conditions are $V_e\left(0\right)=1$, $A\left(0\right)=1$, $V_i\left(0\right)=0$, $R\left(0\right)=0$, $V_r\left(0\right)=0$. From top to bottom: $\left[V_{free}\right]$ corresponds to $V_e$, $\left[gRN{A}_{(+)}\right]$ corresponds to R, and $\left[V_{assembled}\right]$ corresponds to $V_r$. The difference in maximal values of $\left[V_{assembled}\right]$ and $V_r$ is explained by biases due to numerous transformation of the parameters required for reduction in obtaining coefficient $k_3$, which influences the amplitude of $V_r$ accordingly (78), and the absence of a decrease in the end part of the plot is explained by the fact that there is no mortality in the model (1)–(5) for $V_r$ (and $q_2=0$ for these calculations).

Figure 5. Comparison of the kinetics of SARS-CoV-2 replication predicted by model (1)–(5) (right) with results from ([21], Figure 2) (left). The values of parameters are estimated from the model reduction (Table 1): $k_1=0.5$, $q_1=0.5$, ${\sigma }_1=0.12$, ${\sigma }_2=0.2$, $k_3=0.003$, $q_2=0$, $k_4=6.093$, and the initial conditions are $V_e\left(0\right)=1$, $A\left(0\right)=1$, $V_i\left(0\right)=0$, $R\left(0\right)=0$, $V_r\left(0\right)=0$. From top to bottom: $\left[V_{free}\right]$ corresponds to $V_e$, $\left[gRN{A}_{(+)}\right]$ corresponds to R, and $\left[V_{assembled}\right]$ corresponds to $V_r$. The difference in maximal values of $\left[V_{assembled}\right]$ and $V_r$ is explained by biases due to numerous transformation of the parameters required for reduction in obtaining coefficient $k_3$, which influences the amplitude of $V_r$ accordingly (78), and the absence of a decrease in the end part of the plot is explained by the fact that there is no mortality in the model (1)–(5) for $V_r$ (and $q_2=0$ for these calculations).

5. Discussion

5.1. Different Models of Infection Progression

The basic model of infection progression in cell culture represents a reaction–diffusion system of equations with time delay for the concentrations of uninfected cells, infected cells, and viruses [15]. It allows the determination of the main characteristics of infection progression, such as the spreading speed and viral load. This model admits numerous developments, including virus mutation [20], competition [16], and mucus motion in the respiratory tract [17].

All these models treat cell culture as a continuous medium where infected cells produce new viral particles, both of them modeled through their concentrations. Neglecting cell division and motion, as is usually the case in virus plaque assays, we can associate the space variable x not to the local quantity of cells but to the concentrations of some intracellular substances. This approach is not applicable in the case of cell division or motion, which would not allow the description of intracellular concentrations with kinetic equations.

In this work, we combine the macroscopic model of infection spread with a model of intracellular regulation. We deliberately consider a simplified model problem with a one-stage virus replication process. We show how this model can be reduced to the previous macroscopic models. This two-scale approach opens new perspectives in the modeling of viral infections.

Infection spreading in cell culture includes two processes: virus replication inside infected cells and virus motion (diffusion) between cells [9]. The model of intracellular virus replication does not depend on the space variable. It is validated by a comparison with a more complete model of intracellular regulation in [21]. The model with diffusion and space variables is validated by the comparison with the previous model in [20], which was justified by the experimental data.

Let us note that experiments on virus plaque assays are carried out in a monolayer of cells so the corresponding modeling problem is two-dimensional. The comparison of 2D and 1D models was performed in [15]. It was shown that the 1D model gives a good approximation of the 2D model. In the results presented in Figure 2a for the concentration of infected cells, the interval where this concentration is close to 0 can be interpreted as plaque radius.

5.2. Complete and Reduced Models

5.2.1. Comparison with the First Reduced Model

The basic characteristics of reaction–diffusion models, and in particular, of the current model, are the viral replication number $R_v$, total viral load $J\left(v\right)$, and wave speed c. We will compare their values with the corresponding values for the reduced model (46)–(48).

Let us recall that the virus replication number determines the condition of infection progression. In both models, $R_v$ inversely depends on the extracellular virus degradation rate ${\sigma }_1$ and on the degradation rate ${\sigma }_2$. It has the same dependence on the rate of virus penetration into the cell $q_1$ (similar to $k_2$ in the reduced model). It is also proportional to the virus replication rate $k_3$ (similar to $k_1$ in the reduced model). Thus, the qualitative structure of the viral replication number remains the same for both models.

Next, we obtain Equation (32) for $\tilde{w}$, which is the ratio of intracellular components after infection progression to the initial quantity of them. This equation coincides with Equation (65) for the model from Section 3; however, the virus replication number is defined with (11) now. For the total viral load $J\left(v\right)$, which is the integral of the extracellular amount of virus over space, the Equations (34) and (35) coincide with the corresponding equations from the Section 3, irrespective of the new formula for $R_v$. Thus, a more detailed model gives the same estimation for the total viral load as a coarser one.

Comparing the Formula (42) for the wave speed to the (68), we conclude that the numerator now has three factors with degradation coefficients for the viral genomes (DNA or RNA) ${\sigma }_2$, and also for the virus replication rate $q_2$ and for the intracellular viral genome production coefficient $k_1$. In the denominator, the first term is essentially the same, and the last term contains production coefficients not only for extracellular viruses $q_2$ and $k_1$, but also for intracellular virus $q_1$ and virus replication $k_3$. As for the viral replication number $R_v$, the wave speed has the same biological interpretation as in the reduced model, while the formula is more detailed for the complete model.

5.2.2. Comparison with the Second Reduced Model

In the previous work [16], we considered the model of virus propagation in cell culture in terms of cell populations, i.e., at the cell level; it contains the quantities of uninfected cells, infected cells, and viruses. The approach to the model construction suggested in this work is principally different. In this work, the virus quantity and their productions are considered irrespective of types of virus and target cells. This allows us to take into account intracellular reactions, receptors, and features of virus replication, and thus, we obtain more flexible but still quite simple models.

In the model from [16], the virus replication number $R_v$ directly depends on the rate of cell infection a, a rate of virus production b, and the initial quantity of uninfected cells $u_0$, and it inversely depends on the death of infected cells $\beta$ and virus degradation $\sigma$. In the model from Section 3, the $R_v$ directly depends on the rate of virus replication $k_1$, the rate of virus DNA (RNA) production $k_2$, and initial quantity of intracellular components $A_0$, and it inversely depends on the degradation rate ${\sigma }_2$, and the virus consumption, including virus degradation $k_2+{\sigma }_1$. If we correlate the indicated coefficients in the order in which they are listed, we will obtain the same formulas for $R_v$ in both models. Such a correlation is natural since the virus DNA (RNA) is produced in infected cells, so its production and degradation correspond to the production and death of infected cells. Virus replication in the infected cells uses viral proteins located inside the cells. Although the intracellular components of the model from Section 3 refer to infected cells, they act as a resource required for the virus to propagate, just like uninfected cells in the model from [16]. Though the meaning of the coefficients is similar, the quantitative values of the coefficients in both models are different. This should be taken into account for the application of the model and for the comparison of the results. Both models under consideration have the same Equation (32) for $\tilde{w}$; however, the formulas for the total viral load are different. The formulas for the wave speed have the same structure as already been analyzed above.

5.3. Biological Interpretations

The virulence of infection correlates with the size of viral plaques [3,4,5], while this size is determined by the speed of infection spread. Hence, the wave speed in the reaction–diffusion models of viral infection has a direct biological interpretation by means of infection virulence. The analysis presented in this work for the two-scale model of infection progression in cell culture allows us to determine the spreading speed through the parameters of virus transport (diffusion) and of intracellular regulation.

Furthermore, the infectivity of respiratory viral infections, that is, the rate of their transmission from infected to uninfected individuals, is proportional to the viral load in the upper respiratory tract. Therefore, we can use the formulated model to determine how virus infectivity depends on the parameters of the model.

As is usually the case for reaction–diffusion models, the spreading speed is proportional to the square root $\sqrt{D}$ of the diffusion coefficient. Since virions diffuse as inert particles, their diffusion coefficient depends on their size and on the properties of the medium with values of about 10 μm${}^2$/s determined in biophysical experiments.

The parameters of intracellular regulation in the two-scale model include the rates of virus entrance to the cell $q_1$ and its exit from cell $q_2$. The first one is particularly important, and it was largely discussed for the SARS-CoV-2 infection. The spike protein in the new coronavirus undergoes numerous mutations and leads to the emergence of new variants, some of them being more virulent than the original variants. The formula for the wave speed (42) shows that it increases as a function of $q_1$ (Figure 4), but this increase is quite weak.

On the other hand, viral load determined by formula (35) is proportional to the wave speed c and inversely proportional to $q_1$. Figure 3 shows that viral load can increase or decrease for larger values of $q_1$ depending on other parameters. Thus, mutations in the spike protein increasing the virus penetration rate into the host cell lead to the increase in virus virulence but not necessarily an increase in viral load. Note that a comprehensive analysis of the viral load kinetics and infection severity requires consideration of innate and adaptive immune responses (see, for example, [23]).

Overall, the developed family of mathematical models, differing in their complexity, provide a solid analytical basis for further applications in the analysis of the pathogenesis of viral infections via multi-scale and hybrid modeling approaches.

6. Conclusions

Viral infection propagates in cell culture as a reaction–diffusion wave [15,16,20]. In this work, we studied the influence of intracellular regulation of virus replication on infection progression. We determined viral load and spreading speed depending on the rate of virus penetration inside cells and on the rate of virus replication. These results allow the comparison of different SARS-CoV-2 variants for which mutations in the spike protein can influence the entrance rate. Let us recall that viral load in the upper respiratory tract determines the infectivity of respiratory viral infections, that is, the rate of disease transmission from infected to uninfected individuals. On the other hand, the speed of infection spreading determines the size of viral plaques, which correlates with the severity of symptoms. Therefore, we determine how the virus entrance and replication rates influence the infectivity and severity of the disease. The model developed in this work is validated by the comparison with the model of intracellular virus replication [21] and with the model of infection spreading without intracellular regulation [15].