Mathematical modeling of respiratory viral infection and applications to SARS-CoV-2 progression.

Viral infection in cell culture and tissue is modeled with delay reaction-diffusion equations. It is shown that progression of viral infection can be characterized by the viral replication number, time-dependent viral load, and the speed of infection spreading. These three characteristics are determined through the original model parameters including the rates of cell infection and of virus production in the infected cells. The clinical manifestations of viral infection, depending on tissue damage, correlate with the speed of infection spreading, while the infectivity of a respiratory infection depends on the viral load in the upper respiratory tract. Parameter determination from the experiments on Delta and Omicron variants allows the estimation of the infection spreading speed and viral load. Different variants of the SARS-CoV-2 infection are compared confirming that Omicron is more infectious and has less severe symptoms than Delta variant. Within the same variant, spreading speed (symptoms) correlates with viral load allowing prognosis of disease progression.

INTRODUCTION

SARS‐CoV‐2 infection spreads through the epithelial tissue from the upper respiratory tract (URT) to the bronchi and lungs. , Virus replicates in the infected cells that expel viral particles into the extracellular space mainly through the apical surface (towards the airways), , where they diffuse to the neighboring cells. As such, viral infection spreads in the tissue causing further clinical manifestations related to inflammation and tissue damage, while viral load produced in the URT determines infectivity of viral infection by means of virus‐containing droplets. Thus, infection spreading speed and viral load represent important characterizations of viral infection. In this work, we use mathematical modeling in order to determine infection spreading speed and viral load, and we apply these results to characterize different variants of SARS‐CoV‐2 infection from the point of view of their infectivity and disease severity.

Bronchial epithelium is covered with the airway surface luquid (ASL) which consists of periciliary layer (PCL) and mucus produced by goblet cells. Coordinated motion of cilia at the surface of cilia cells moves ASL towards the pharynx. Mucus from the URT also moves to the pharynx but from the other side. Mucus protects the epithelial tissue from pathogens removing them with convective flow. It can also neutralize virus due to mucin molecules. SARS‐CoV‐2 can infect goblet and cilia cells and influence mucus production and motion. Moreover, mucus becomes more viscous in the infected tissue due to debris of viral particles and dead infected cells and because more viscous mucus is produced due to infection and inflammation. Furthermore, convective motion of ASL can influence infection spreading in the tissue. Thus, mucus motion should be taken into account in modeling respiratory infections.

Infection progression in cell cultures and tissues can be modeled with reaction‐diffusion equations for the concentrations of uninfected cells, infected cells, and virus, , while the innate and adaptive immune responses can be taken into account through the concentration of interferon, antigen presenting cells, and lymphocytes (in various ODE models; see, e.g., Lee ). Analysis of these models shows that viral infection spreads in cell tissue as a reaction‐diffusion wave. , , . Such waves are mainly characterized by two parameters: the quantity of produced virus (viral load) and the speed of propagation. We determine them through the original model parameters, such as the rate of cell infection by virus and the virus replication rate in the infected cells, and we show that they are independent of each other in the sense that high (or low) viral load can be combined with high or low wave speed.

High viral load in the URT can be associated with a high rate of disease transmission. , However, low wave speed slows down disease progression in the lungs and, consequently, decreases disease severity. Recent animal studies show that, indeed, lung infection in the case of Omicron variant is low or not present at all. Moreover, cell culture and tissue experiments show that infection progression depends on cell types. , We use experimental data on time‐dependent viral load from Chanet al and Peacock et al in order to identify model parameters and to characterize infection progression. In particular, we show that Omicron variant spreads faster in the URT than Delta variant and slower in the lungs, resulting in higher infectivity, shorter incubation period, and less severe symptoms.

We model infection progression in the epithelial tissue taking into account mucus motion in a two‐phase reaction‐diffusion model with a fluid flow in a two‐dimensional strip and one‐dimensional epithelial cell layer at the boundary of this strip (Section 2). We show that this model can be reduced to a one‐dimensional two‐phase model. This reduction allows us to find analytical expressions for the viral load and wave speed (Section 3). We apply these results to characterize Delta and Omicron variants in Section 4 and discuss the results in Section 5.

INFECTION PROGRESSION IN THE EPITHELIAL TISSUE

Bronchial epithelium is covered with ASL consisting of PCL and mucus and transported by the coordinated motion of cilia at the surface of cilia cells. In the process of viral infection, viral particles produced in the epithelial cells are expelled through the apical cell surface in the extracellular space where they diffuse and they are also transported by convective flow. Hence, studying infection progression in the epithelial tissue, we need to take into account ASL motion.

2D model

We consider a 2D problem in the domain . Epithelial cell layer is considered as 1D interval located at the lower boundary , fluid is located in the whole rectangular domain, and is the width of the fluid layer (mucus and PCL). Approximation of constant layer width implies that we do not consider mucus production and accumulation. This question will be considered in the forthcoming work.

Virus production in infected cells

We consider the system of equations for the concentration of uninfected cells, the concentration of infected cells, and the virus concentration in the cell layer that includes epithelial cells and adjacent extracellular space. Next, is virus concentration in fluid at the lower boundary, and . Parameter characterizes the rate of cell infection, corresponds to the virus replication rate in the infected cells, is the death rate of infected cells, is the rate of virus death, parameters and characterize virus exchange between cell layer and fluid layer, and is virus diffusion coefficient in the cell layer. All these parameters are positive constants. A similar model without fluid was considered in the recent works. ,

Virus distribution in fluid

Virus distribution in the fluid layer is described by the diffusion‐convection equation where is the diffusion coefficient, is fluid velocity considered as a given function which can be approximated by a linear function in the PCL and constant in the mucus layer, and is the virus death rate in the fluid layer. This equation is completed by the boundary conditions: where is the cell height and initial condition for which corresponds to a virus‐containing droplet in the fluid layer.

Stationary distribution

We begin the study of problem (2.1)–(2.5) with the analysis of its stationary solutions in the cross‐section of the layer (independent of ). It has a trivial stationary solution . Neglecting death of infected cells ( ), we can find another stationary solution, independent of the space variable and satisfying the equations where , and the boundary conditions:

We note that Equations (2.1) and (2.2) for give the balance equation , where is the initial concentration of uninfected cells. This balance equation allows us to reduce these two equations to equation (2.6). Searching stationary solution of this problem, we find from Equation (2.6) and . Then we obtain a closed problem for the stationary solution (we keep for convenience the same notation):

The solution of this problem is given by the following expression: where Let us note that for the considered range of parameters (e.g., Figure 1), the variation of the function in the interval is very weak (at most 0.1 %). This property of solution will allow us to reduce the 2D problem to a 1D problem.

FIGURE 1

Virus distribution in infected cell (left) and the concentration of uninfected cells (right) in consecutive moments of time (every 1 h). The values of parameters are as follows: [Colour figure can be viewed at wileyonlinelibrary.com]

Wave propagation

Viral infection propagates in the epithelial tissue as a reaction‐diffusion wave if the influence of the boundary is negligible. From the mathematical point of view, it is a solution of problem (2.1)–(2.5) of the form considered for all real with the limits for , and at for . The values and are determined in the previous section.

Propagation of viral infection in the epithelial layer described by models (2.1)–(2.5) is shown in Figure 1 for the values of parameters determined in Ait Mahiout et al for the SARS virus in the culture of cilia cells (Sims et al. ). It is a typical reaction‐diffusion wave of infection spreading. Damped oscillations are caused by the time delay in virus replication in the infected cells. Two‐dimensional virus distribution in the fluid layer is shown in Figure 2. Its dependence on is weak (less than 1%) since the fluid layer is narrow. This allows us to reduce the 2D problem two one‐dimensional two‐layer problem.

FIGURE 2

A snapshot of virus distribution in fluid for the same values of parameters as in Figure 1 [Colour figure can be viewed at wileyonlinelibrary.com]

Reduction to 1D model

Set and integrate equation (2.4) with respect to . Then we obtain

Since the dependence of on is weak, we can use the approximation . Then the last equation can be written as follows: where is an average flow velocity. Hence, we obtain the following 1D model: where is virus concentration in cells, its concentration in fluid, and is an average fluid velocity.

ONE‐DIMENSIONAL TWO‐LAYER MODEL

Virus replication number and viral load

Experimental data and mathematical modeling show that infection spreads in cell culture as a reaction‐diffusion wave. In order to study this type of solutions, we consider system of Equations (2.15)–(2.18) on the whole axis and set , where is the wave speed. Then we obtain the following system of equations: where prime denotes the derivative with respect to the variable . We are interested in solution of this system of equations with the following limits at infinity:

In this section, we will determine and viral load setting for brevity. From Equation (3.1), we obtain

Taking a sum of (3.1) and (3.2) and integrating

Next, from (3.3), and from (3.4)

We denote the corresponding integrals by and obtain from (3.8) and (3.9): We substitute these expressions into (3.6) and (3.7) and finally obtain the equation with respect to . Here, is the virus replication number. Equation (3.10) has a solution in the interval if and only if . The total viral load satisfies the equality where is the wave speed. We can now determine .

Viral infections are usually characterized by large values of . In this case, solution of Equation (3.10) is close to 0, and . Therefore, expression (3.12) for the total viral load can be written in the following form: We will see in the next section that the wave speed depends on and , while the remaining factor in (3.13) depends only on (for other parameters fixed). Therefore, the wave speed and the total viral load can vary independently of each other, and they represent two different characterizations of viral infection.

Wave speed

Analytical formula for a simplified model

In order to obtain an analytical formula for the wave speed, consider a simplified model which does not take into account the difference in virus concentrations in the layer of epithelial cells and in the adjacent fluid layer. Using the linearization method, , we consider traveling wave solution , and we replace by its value at . Therefore, we obtain the following system of equations for : where we keep notation for the independent variable and omit for simplicity of notation. We look for a solution of this system in the form . Substituting them in (3.17) and (3.18), we obtain We should find the minimal value of for which this system of equations has a positive solution . Introducing an independent parameter and excluding and , we obtain the following equation: Hence, where is a positive solution of the equation . It is proved in Ait Mahiout et al that for , this formula gives the exact value of the minimal wave speed. In general, it gives an estimate of the minimal wave speed from below. Numerical simulations confirm that we obtain the exact value of the minimal speed also for . This formula allows us to assess the dependence of the wave speed on parameters and to obtain an approximation of the wave speed for the two‐layer problem.

Approximation of the wave speed for the two‐layer problem

Comparison with numerical simulations shows that the linearization method applied to the two‐layer problem (2.15)–(2.18) gives an estimate from below for the minimal speed and not the exact value. Moreover, this estimate can be essentially different from the minimal speed, and it cannot be used as a reliable approximation. In order to obtain a better approximation of the minimal speed in this case, consider system (2.15)–(2.18) assuming for simplicity that . Set . Suppose that . Then we obtain the following system of equations: where . We approximate by a constant and determine it from the limiting values of and at (for ): . We can now apply formula (3.19) to this system. It gives a good approximation for 1D and 2D numerical results (Figure 3, left).

FIGURE 3

Left: wave speed in numerical simulations of system (2.15)–(2.18) (curves 1,2,3), in the analytical approximation with formula (3.19) (curves 1,3), and in 2D simulations (dots). In the 1D case, numerical results coincide with the analytical approximation given by formula (3.19); thus, the corresponding approximating curves overlap. The values of parameters are as follows: ; 1. ; 2. ; 3. . Right: viral load in in system (3.14)–(3.16) as a function of time. The values of parameters: (upper curve), (middle curve), (lower curve) [Colour figure can be viewed at wileyonlinelibrary.com]

Dependence of the wave speed on parameters

From the analytical and numerical results, we can make the following conclusions about the dependence of the wave speed on parameters.

Dependence on and . Formula (3.19) shows that the wave speed depends on the product and not on the individual parameters otherwise. Numerical simulations for system (2.15)–(2.18) confirm this conclusion for the whole range of parameters. Changing and up to the three order of magnitude in such a way that their product remains the same, the wave speed changes in the limit of 2 %. Dependence of the wave speed on is shown in Figure 3. This dependence is weak, the wave speed changes at most 1.5 times when changes by four orders of magnitude.

Dependence on . New viral particles produced by infected cells are expelled into the extracellular matrix through their apical surface adjacent to the fluid layer. We will consider two different cases: either viruses diffuse only in the fluid ( ) or they also diffuse along the cell surface ( ). The difference between these two cases is not essential in the absence of fluid flow ( ), but it becomes essential for physiological values of fluid velocity (see below).

Dependence on . If virus diffusion occurs only in the fluid, then the wave speed dependence on the flow velocity is essential. For the flow velocity 10 cm/h (lower estimate of physiological values), the wave speed against the flow becomes less than 10−3 cm/h (for and other parameters as in Figure 3). This means that viral infection will not arrive to the lungs from the URT during given time interval (in average 9 days according to the clinical data). Therefore, we need to assume nonzero diffusion in the cell layer. Furthermore, it should be of the same order of magnitude 0.001 cm2/h as we take for fluid; otherwise, the wave speed remains too small. If , the wave speed against the flow does not depend on the flow velocity.

Dependence on other parameters. Dependence of the wave speed on the replication delay is quite essential (Figure 3). Its dependence on is weak and not essential in the evaluation of infection progression.

Dependence of the viral load on parameters

The value of viral load during the third stage of infection progression (wave propagation) is given by the analytical formula (3.13). If we take the value of the wave speed in numerical simulations, then we obtain exact correspondence between the analytical and numerical values of the viral load. An analytical approximation of the viral load during the previous two stages (decay due to the replication delay, explosive growth) can be constructed as a solution of delay differential equations. Dynamics of viral load in numerical simulations is shown in Figure 3 (right) for different values of the cell death coefficient .

CHARACTERIZATION OF DELTA AND OMICRON VARIANTS

We compare the results of numerical simulations of infection progression in the epithelial tissue (without fluid motion) with the experimental data from Chan et al and Peacock et al. Experimental results in human bronchial tissue show that viral load for the Omicron variant of the SARS‐CoV‐2 infection is larger than for the Delta variant during the first two days, but it becomes opposite at 72 h. Parameters of the simulations are chosen to fit the data (Figure 4).

FIGURE 4

Total viral load as a function of time in the experiments in bronchial tissue and in numerical simulations for Delta (blue line, blue crosses) and Omicron (red line, red asterisks) variants. Left: bronchial tissue. The values of parameters for Delta: , for Omicron: . The initial condition is for and 0 otherwise. Experimental data are taken from Chan et al. Parameter dimensions are the same as in Figure 3. Tissue culture infectious dose (TCID50) is a measure of virus titers. Right: lung tissue. The values of parameters for Delta: , for Omicron: (other parameters are the same). Experimental data are taken from Chan et al. Parameter dimensions are the same as in Figure 3 [Colour figure can be viewed at wileyonlinelibrary.com]

Since the values of parameters are important in order to characterize infection progression, let us briefly discuss how they are determined with the example of Figure 4 (left). The total viral load as a function of time has three consecutive stages: exponential decay before the beginning of virus replication, explosive growth in the beginning of replication, more slow growth during infection spreading. There are three steps of growth of viral load before it stabilizes at the maximal value (it can decrease or grow later depending on the value of ). The maximal viral load at 72 h determines the value of . The beginning of explosive growth depends on the value of , while the size of the first step (for a given value of ) depends on the value of . Since the value of viral load for Omicron at 24 h is not very different from the next two points, then should be at the second step, and therefore, the time delay of replication is less than 12 h. It is different for Delta, since should be at the first step, and therefore, is larger than 12 h. The best results are obtained for h and h. These values can depend on the experimental conditions, and they can be different from in vivo replication delays.

Let us note that the value of parameter reflects the fact that only minor proportion of virus particles penetrate uninfected cells. This may be related to various cell protection mechanisms , and to the experimental conditions (cf. Sims et al. ).

Analysis of the simulations allows us to conclude that Omicron starts to replicate faster than Delta in the bronchial tissue but slower in the lung tissue. The viral load is larger for Omicron during the first 2 days in the bronchial tissue and smaller in the lung tissue (all 3 days). The speed of infection spreading in the bronchial tissue is cm/h and cm/h for Omicron and Delta, respectively. In the lung tissue, they are cm/h and cm/h. Thus, Delta variant preserves its spreading speed while it essentially decreases for Omicron.

It is interesting to note that the final value of the viral load for Delta in the bronchial tissue is larger than for Omicron, but its spreading speed is lower due to larger replication time. The viral load slowly decreases during further infection spreading due to death of infected cells. The virus replication number for both variants is the same in the bronchial tissue , but it is different in the lung tissue, . We see that virus replication number is not sufficient to determine viral load and spreading speed. These three characteristics are independent of each other.

Comparison of viral load for Omicron and Delta variants in cell culture of human nasal epithelial cells and lung cells is presented in Peacock et al. Parameter fitting and modeling of these results (not shown) confirm the conclusions about spreading speeds of the two variants in epithelial and lung tissues.

DISCUSSION

Comparison of SARS‐CoV‐2 variants

Comparison of Delta and Omicron variants is carried out in this work on the basis of experimental data on viral load in cell cultures and tissues. , Different experimental conditions can influence the time‐dependent viral load, and it can also differ from in vivo values. In particular, the replication delay in these experiments (10–20 h) exceeds the estimates 6–7 h in Bridges et al and Harcourt et al. Furthermore, the corresponding spreading speed is smaller than previously estimated for the SARS virus in Ait Mahiout et al on the basis of the experiments in Sims et al. The latter is about 0.06 cm/h (1.5 cm/day) in agreement with 9 day time interval for infection spreading from the URT to the lungs. Parameter estimates are consistent with other modeling studies.

However, different experimental results converge in relative comparison of viral loads for Delta and Omicron variants in the epithelial and lung tissue. Hence, based on different experiments, we validate the conclusion that the spreading speed of Omicron variant is larger than the spreading speed of Delta variant in the epithelial tissue and smaller in the lung tissue.

Virus diffusion and convective transport

The estimates of infection spreading speed presented above allow us to evaluate the contribution of diffusive and convective transport. An important conclusion from this analysis is that virus diffusion occurs in the cell layer (cells and adjacent extracellular space) where it is not influenced by convective transport. Furthermore, infection spreading speed from URT to the lungs (against fluid flow) does not essentially depend on the flow velocity due to diffusion in the cell layer. However, infection spreading in the opposite direction can be essentially accelerated by the flow.

Symptoms and infectivity

Clinical manifestation of viral infection related to tissue damage is not determined by the viral load itself but by the part of the infected tissue, that is, by the speed of infection spreading. As we discussed above, these are two different characteristics, though both of them are determined by the original model parameters. Hence, higher viral load can be associated with the same or even smaller wave speed. Indeed, since the wave speed depends on the product , we can vary these two parameters in such a way that their product remains constants. Then viral load is a linear function of (see 3.13), and for the same speed, we can obtain any viral load. Thus, more sever or rapid symptoms can be observed for smaller viral load for different virus variants.

For a given variant, the parameters of the model are fixed, but inter‐patient variation can be determined by the virus replication rate that depends, in particular, on the interferon concentration. Increasing leads to the increase of the viral load (linearly proportional to ) and of the spreading speed ( ). Hence, there is a correlation between the value of viral load in the URT and disease severity, in agreement with clinical observations. These results can be used for the prognosis of the disease progression.

Comparing the speeds of propagation for the Delta and Omicron variants in different tissues, we can conclude that the incubation period for Omicron is about twice shorter than for Delta, and its spreading in lungs is about twice slower.

Immune response

In the innate immune response, interferon produced by infected cells slows down virus replication rate. It is implicitly taken into account through the parameters of the model (cf. Ait Mahiout et al. ). The influence of the adaptive immune response, which manifests itself on later stages of infection progression, will be considered in further works. We will restrict ourselves here by some remarks related to the previous analysis.

Denote by the virus replication number taking into account the adaptive immune response. It is obtained from replacing by , where is the rate of infected cell elimination by cytotoxic lymphocytes CTL, and by , where characterizes virus neutralization by antibodies. Adaptive immune response stops infection progression if . This condition, that is, the efficacy of the adaptive immune response in stopping infection, is not directly determined by the viral load. Both of them can be expressed through the original model parameters, but they are not necessarily proportional to each other. Indeed, it follows from (3.13) that viral load is inversely proportional to and can adopt high values for small enough. Since , then depends on , but it weakly depends on . Hence, the same action of the adaptive immune response can be obtained for very different viral loads.

Conclusions and perspectives

In this work, we have determined the viral load and infection spreading speed for a generic respiratory viral infection and applied these results to characterize different variants of SARS‐CoV‐2 infection. More detailed description of mucus motion and of the immune response will be considered in the forthcoming works.

CONFLICT OF INTEREST

This work does not have any conflicts of interest.