Role of error catastrophe in transmission ability of virus

ABSTRACT

The role played by error catastrophe is explicitly taken into account in a mathematical formulation to analyze COVID-19 data. The idea is to combine the mathematical genetics formalism of the error catastrophe of mutations in virus gene loci with the standard model of epidemics, which lacks the explicit incorporation of the effect of mutation on the spreading of viruses. We apply this formalism to the case of SARS-CoV-2 virus. We assume the universality of the error catastrophe in the process of analyzing the data. This means that some basic parameter to describe the error catastrophe is independent of which group (country or city) we deal with. Concretely, we analyze Omicron variant data from South Africa and then analyze cases from Japan using the same value of the basic parameter derived in the South Africa analysis. The excellent fit between the two sets of data, one from South Africa and the other from Japan, using the common values of genetic parameters, justifies our assumption of the universality of these parameters.

INTRODUCTION

In usual analysis of virus spreading based on the standard model of epidemiology by Kermack and McKendrick (1927), we ignore the possible role played by mutations, at least in the explicit form of the model. One of the present authors published two papers (Sugawara, 2021a, 2021b) analyzing the world-wide COVID-19 situation based on the standard model, through which he learned that the standard model does not properly incorporate the role played by mutations. Because of this, some difficulty was encountered in interpreting the result of the analysis.

A simplistic assumption of the absence of mutations in the virus genome is still made even for recently developed models in which, for example, multiple pathogen strains interact with each other through cross-immunity (see Kucharski et al., 2016 for a review) or strains interfere through the removal of individuals from the susceptible pool after an acute infection (e.g., Rohani et al., 2002). Meanwhile, the importance of the random nature of pathogen mutations and the emergence of new strains have been well recognized and studied by simulation under the assumption that a host’s risk of infection is proportional to the genetic distance between the infecting phenotype and the closest phenotype in the host’s immune history (Bedford et al., 2012). Recently, Miller et al. (2022) explored an existing agent-based model for influenza A/H3N2 that incorporates mutation effects into both the transmission rate and the antigenic space in the standard model.

We present here a way to remedy the current situation for the presence of largely deleterious but occasionally advantageous antigenic mutations by introducing into the scheme of the standard model some elements showing the vital role played by mutations before, during and after the period of extensive virus spreading. When we analyze the data for variant strains of SARS-CoV-2 such as Alpha, Beta, Gamma, Delta, Omicron BA.1 and Omicron BA.2, we notice that temporal changes of the number of infected hosts show an interesting similarity among the different variant strains. In this regard, we emphasize an outstanding role played by error catastrophe (Eigen, 1971): viruses can improve their transmission ability by mutations, but it seems that there is a maximum ability they can reach due to errors during the DNA/RNA replication processes. After they reach the maximum (error threshold), both spontaneous and forced mutations (by an enzyme such as APOBEC3 (Sadeghpour et al., 2021), as was recently proposed) continue to cause errors in the copying process of the virus multiplication and thus reducing the transmission ability.

In the following sections, we first explain the biological justification of our model in more detail. We then formulate the virus mutation process in mathematical terms, and go on to apply this mathematical scheme to actual data analysis, taking infection data from Tokyo as an example. The last section is devoted to the analysis of Omicron data from South Africa and Japan. Appendices 1 and 2 deal with the issue of symbiosis and the generality in our formalism, respectively.

BIOLOGICAL JUSTIFICATION OF THE MODEL

As stated in the Introduction, analyses of the data for variant strains of SARS-CoV-2 show that temporal changes of the number of infected hosts exhibit an interesting similarity among the different variant strains. So far, no clear explanation of this phenomenon exists within the framework of mathematical epidemiology. The data also show similar decaying patterns of the infection numbers for a particular variant once they reach the maximum in various countries where different epidemic policies are adopted, including vaccination rates. These similarities especially of the decaying patterns, which are independent of the kind of variant strains or of epidemic policies, seem to be explainable by assuming that the reason is on the side of viruses concerned.

A possible model is what we call an error catastrophe model (Eigen, 1971), which claims that errors accumulated in the form of mutations make it hard for the virus multiplication process to be smoothly performed. The infection rate of a virus such as SARS-CoV-2, which has a high rate of mutation, will be strongly influenced by the errors and its decaying patterns are largely determined by them. In the case of SARS-CoV-2, which allows APOBEC3 enzymes to directionally mutate C to U, error catastrophe may be particularly important. We do not need to treat this RNA editing by APOBEC3 distinctly from the usual mutations that occur during replication (Sadeghpour et al., 2021). It simply emphasizes the importance of error catastrophe and, therefore, it is hidden in our mathematical formulation.

The genome size of SARS-CoV-2 is about 30 kb, three times larger than that of typical RNA viruses such as influenza viruses and HIV. Despite being an RNA virus, SARS-CoV-2 accumulates mutations comparatively slowly because the virus displays high-fidelity replication by expressing the 3′-5′ exoribonuclease nsp14-ExoN (Smith et al., 2013). However, if nsp14-ExoN is disrupted by self-mutagenesis or external agency, the mutation rate may increase 15-fold (reviewed by Robson et al., 2020; Thakur et al., 2022). Nearly half of the nucleotide changes in SARS-CoV-2 are C to U transitions and they occur preferentially in both 5′U/A and 3′U/A flanking sequence contexts that are favored motifs of the human APOBEC3 protein family of cytidine deaminases. As a result, SARS-CoV-2 may mutate at a rate of 3×10⁻⁴/site/year, or 9/genome/year, consistent with a mean of 5.5–9.5 nucleotide differences between variants (Simmonds, 2020; Ratcliff and Simmonds, 2021). Such abundant and context-dependent mutations likely hinder proper replication and lower the power of infection, although by unknown molecular mechanisms, and allow maladaptive amino acid changes to accumulate, leading to suicidal explosion of the virus or error catastrophe. Even so, antigenic escape can occur and challenge pandemic control efforts. This unique paradigm of SARS-CoV-2 genome evolution is not incorporated by standard models of molecular epidemiology.

In the past, efforts have been made to incorporate pathogen evolution into the framework of epidemic models (SIR or its extension); for example, see section 3 of Kucharski et al. (2016). Compared to the reports cited there, our approach is unique in its emphasis on error catastrophe.

MATHEMATICAL FORMULATION

The purpose of this section is to combine mathematical formulation of the mutation process with the standard model of mathematical epidemiology. We start with a formulation of the mutation process and properly combine it with the standard model of epidemiology. Our target space is a well-defined population of humans such as a country or a city. We concentrate on the genes that are relevant for the transmission of the virus. This means that we focus mostly, but not exclusively, on mutations in the spike proteins.

We start by defining x_i to be the number of viruses that belong to the i-th variant. Here, “variant” means not just a prominent variant such as Delta or Omicron, but any population of viruses that has mutations distinguishing it from the original virus. Let a_i denote the rate at which variant i reproduces (a_i does not necessarily mean the infection rate of variant i), and let Q_ij denote the probability of a virus of variant i mutating to variant j. The rate of change of x_j is then given by   

$$\frac{d{x}_j}{dt}={\sum }_i{a}_i{Q}_{ij}{x}_i.$$

Suppose we group one kind of variants that have higher infection rates into type X and another kind that have lower infection rates into type Y. We respectively denote y_i, b_i and R_ij rather than x_i, a_i and Q_ij for the latter case. Then, we obtain   

$${\sum }_j\frac{d{x}_j}{dt}={\sum }_{i,j}{a}_i\left(1-Q_{ij}\right)x_i+{\sum }_{i,j}{b}_i{R}_{ij}{y}_i$$

and   

$${\sum }_j\frac{d{y}_j}{dt}={\sum }_{i,j}{b}_i\left(1-R_{ij}\right)y_i+{\sum }_{i,j}{a}_i{Q}_{ij}{x}_i.$$

To simplify the discussion, we assume that a_i, b_i, Q_ij, and R_ij are all independent of suffix “i”. We then obtain   

$$\frac{dx}{dt}=a\left(1-Q\right)x+bRy$$[3.1]

and   

$$\frac{dy}{dt}=b\left(1-R\right)y+aQx.$$[3.2]

Here, we have   

$$x={\sum }_j{x}_j,\mathrm{\hspace{0.17em} }y={\sum }_j{y}_j,\mathrm{\hspace{0.17em} }Q={\sum }_j{Q}_{ij},\mathrm{\hspace{0.17em} }R={\sum }_j{R}_{ij}.$$

In the two-component vector form, Eqs. [3.1] and [3.2] become   

$$\frac{d}{dt}\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}a\left(1-Q\right)& bR\\ aQ& b\left(1-R\right)\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)\equiv \mathrm{\Omega }\left(\begin{array}{c}x\\ y\end{array}\right).$$[3.3]

Here, 1+a is the average reproduction rate of variants in type X and 1+b is that for type Y. Both a and b can be positive or negative considering replication errors. We note that all the parameters in the matrix Ω may depend on time t. However, if the time dependence is small, we can use an adiabatic approximation by considering the constant parameters that will be given a small time dependence after solving the differential equation [3.1], [3.2] or [3.3].

The number of infected hosts I at time t of our target population will change according to the following formula:   

$$\frac{dI}{dt}={\omega }_{x}x-{\omega }_{y}y$$[3.4]

where ω_x and ω_y are constants. The eigenvalues of matrix Ω are given by   

$$det\left(\begin{array}{cc}a\left(1-Q\right)-\xi & bR\\ aQ& b\left(1-R\right)-\xi \end{array}\right)=0$$[3.5]

which gives   

$$\begin{array}{l}\xi =\\ \frac{1}{2}\left[a\left(1-Q\right)+b\left(1-R\right)\pm \sqrt{{\left\{a\left(1-Q\right)-b\left(1-R\right)\right\}}^2+4abQR}\right].\end{array}$$[3.6]

In general, ξ is a complex number:   

$$\xi ={\xi }_R\pm i{\xi }_I$$[3.7]

with   

$${\xi }_R=\frac{1}{2}\left\{a\left(1-Q\right)+b\left(1-R\right)\right\}$$[3.8a]

and   

$${\xi }_I=\frac{1}{2}\left[\sqrt{-4abQR-{\left\{a\left(1-Q\right)-b\left(1-R\right)\right\}}^2}\right].$$[3.8b]

We obtain   

$$\mathrm{\Omega }=P^{-1}\left(\begin{array}{cc}{\xi }_R+i{\xi }_I& 0\\ 0& {\xi }_R-i{\xi }_I\end{array}\right)P$$[3.9]

where P is the rows of the eigenvectors of Ω. Eq. [3.3] gives   

$$\left(\begin{array}{c}x\\ y\end{array}\right)=P^{-1}\left(\begin{array}{c}A\text{exp}\{\int ({\xi }_R+i{\xi }_I)dt\}\\ B\text{exp}\{\int ({\xi }_R-i{\xi }_I)dt\}\end{array}\right).$$[3.10]

The general real solution is given as   

$$x=C\text{exp}\left({\int }_0^t{\xi }_{R}dt\right)\text{sin}\left({\int }_0^t{\xi }_{I}dt+k_1\right)$$[3.11a]

and   

$$y=D\text{exp}\left({\int }_0^t{\xi }_{R}dt\right)\text{sin}\left({\int }_0^t{\xi }_{I}dt+k_2\right)$$[3.11b]

where C and D are constants. We may assume that the time dependence of ξ_R+iξ_I is not large and we ignore the contributions after the linear term in the power series expansion of the eigenvalues in terms of t for the later analysis of the data. Although a more precise treatment such as perturbation theory developed in particle physics is available, it will be explored elsewhere. We note that we have   

$$\begin{array}{l}tan\left(k_2-k_1\right)=\frac{{\xi }_I}{{\xi }_R-a\left(1-Q\right)}\\ =\frac{\sqrt{-4abQR-{\left\{a\left(1-Q\right)-b\left(1-R\right)\right\}}^2}}{-a\left(1-Q\right)+b\left(1-R\right)}\end{array}$$[3.12a]

and   

$$D=\left|\frac{bQ}{aR}\right|C.$$[3.12b]

Principle of error catastrophe

We see from Eqs. [3.11] that both x and y have an exponential factor $exp\left({\int }_0^t{\xi }_{R}dt\right)$. This means (due to Eq. [3.4]) that for positive ξ_R, the infection number behaves wildly, oscillation becoming larger and larger in size. For vanishing ξ_R, we have a periodical oscillation as in a normal influenza virus infection. The error catastrophe is defined to be the case when we have ξ_R < 0. A condition for this to happen for time-independent ξ_R is given by   

$${\xi }_R=\frac{1}{2}\left\{a\left(1-Q\right)+b\left(1-R\right)\right\}<0.$$[3.13]

We assume that Eq. [3.13] is valid for SARS-CoV-2 virus, at least for x > x₀ (error threshold). There will be periods when this condition is not valid, especially at the beginning of a surge of certain virus variants. The change of the infection number given in Eq. [3.4] is written as   

$$\frac{dI}{dt}={\omega }_{x}x-{\omega }_{y}y={\text{e}}^{{\xi }_{R}t}\left[\rho sin\left({\xi }_{I}t+k_1\right)-\lambda sin\left({\xi }_{I}t+k_2\right)\right]$$[3.14a]

and   

$$\rho ={\omega }_{x}C,\mathrm{\hspace{0.17em} }\lambda ={\omega }_{y}D.$$[3.14b]

Here, for simplicity, we assumed the eigenvalues to be time independent. An adiabatic approximation will allow us to include moderate time dependence of the eigenvalues and it will be used in our actual analysis of the data in later sections. To see the behavior of dI/dt for a positive or a negative ξ_R, we plot Eq. [3.14a] in Fig. 1A and 1B for two cases.

Fig. 1.

Behavior of Eq. [3.14a] for daily cases (the ordinate) against days (t) elapsed since the beginning of infection (the abscissa). (A) ξ_R = 0.2, ξ_I = 0.1, ρ =5 , λ = −0.1, k₁ = 0, k₂ = 0.83. (B) ξ_R = −0.1, ξ_I = 0.1, ρ = 1, λ = 0.1, k₁ = 1.0, k₂ = 0. (C) ξ_R = 0.2 for 0 ≤ t ≤ 17 and ξ_R = −0.1 for 17 ≤ t ≤ 39. The other parameters are common: ξ_I = 0.1, ρ = 1, λ = 0.1, k₁ = 1.0, k₂ = 0.

The values chosen in Fig. 1 correspond to the following basic parameters:   

$$Q=R=\frac{3}{4},\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }b=-a+\epsilon$$[3.15a]

with   

$${\xi }_R=\epsilon \left(1-Q\right)=\frac{\epsilon }{4}$$[3.15b]

  

$${\xi }_I=\sqrt{\frac{1}{2}\left\{{\left(a-\frac{\epsilon }{2}\right)}^2-\frac{3{\epsilon }^2}{8}\right\}}$$[3.15c]

  

$$tan\left(k_2-k_1\right)=\frac{{\xi }_I}{{\xi }_R-\frac{a}{4}}$$[3.15d]

and ε = 0.8 for case (A) in Fig. 1A and ε = −0.4 for case (B) in Fig. 1B. We also plot Eq. [3.14a] when ξ_R changes from 0.2 to −0.1 as the spreading peaks. Concretely, we choose (A) initially, and (B) immediately after the peak of spreading. We see from these plots that:

(1) There is a maximum value for the change of infection number for either case of ξ_R > 0 or ξ_R < 0

(2) Whether the positive ξ_R is retained (A) or it changes to a negative value (B), that is, whether dI/dt chooses the case depicted in Fig. 1C, depends on the actual data.

Incidentally, we observe that any city or country, except perhaps for China, has experienced multiple peaks in the COVID-19 daily infection number. It is likely that each peak is caused by a new variant, which is sometimes drastically new, such as the Delta or Omicron type, and sometimes just a small number of mutations different from the preceding variant.

Incorporating the standard model of epidemiology

Our method to combine the above formalism on mutation and the standard model of epidemiology, often called the SIR model as mentioned earlier, is straightforward. The SIR model is defined to be the following set of differential equations for S (susceptible), I (infected) and R (recovered):   

$$dS/dt=-\beta IS/N$$[3.16a]

  

$$dI/dt=\beta IS/N-\gamma I$$[3.16b]

  

$$dR/dt=\gamma I.$$[3.16c]

Here, β IS/N is the rate of decrease of the susceptible S, γ is the recovery rate and S + I + R = N. We assume that N is a constant (the total population of a body we want to discuss). On the other hand, we have Eq. [3.4], which defines the infection numbers in terms of two types of variants X and Y defined by Eqs. [3.1] and [3.2]. By regarding Eqs. [3.16a] and [3.16c] as defining S and R, the only useful equation is Eq. [3.16b]:   

$$\frac{dI}{dt}=\beta \frac{IS}{N}-\gamma I.$$[3.17]

Assuming that most of the members of the group are susceptible, that is, S/N ≈ 1, we obtain   

$$\frac{dI}{dt}=\left(\beta -\gamma \right)I.$$[3.18]

Our idea is that we additively combine epidemiologic Eq. [3.18] with genetic Eq. [3.4]. Unlike Bedford et al. (2012), an implicit assumption is that the infection rate is solely determined by pathogen variants and independent of both the host immune history and possible time-dependent antigenic mutation rates. We further assume that infection spreads not only from the officially recognized infected in proportion to I, but from the officially unrecognized infected and/or partly as zoonosis. We then have the equation   

$$\frac{dI}{dt}=\mathrm{\Gamma }I+{\omega }_{x}x-{\omega }_{y}y$$[3.19]

with   

$$\mathrm{\Gamma }=\beta -\gamma .$$[3.20]

The daily infection number is defined as   

$$\mathrm{D}\mathrm{a}\mathrm{i}\mathrm{l}\mathrm{y}\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\equiv \frac{dI}{dt}.$$[3.21]

The condition that the daily infection number is non-negative is given by $\frac{dI}{dt}\ge 0$. This gives I ≥ (−ω_xx + ω_yy)/Γ, which can be satisfied as long as the accumulated daily infection is large enough. However, sometimes it may happen that the inequality is not satisfied. In this case, we can fit the data of this period separately from the other periods.

We solve linear differential Eq. [3.19] for I using the ordinary method of variation of parameters. We assume I in the form of I(t) = f(t)Γe^Γt, substitute it into both sides, and obtain the result of f(t):   

$$\begin{array}{l}\frac{dI}{dt}=\frac{df\left(t\right)}{dt}\mathrm{\Gamma }{e}^{\mathrm{\Gamma }t}+\mathrm{\Gamma }I=\mathrm{\Gamma }I+{\omega }_{x}x-{\omega }_{y}y\\ \mathrm{\hspace{1em}\hspace{0.17em}}=\mathrm{\Gamma }{e}^{\mathrm{\Gamma }t}\left[f\left(0\right)+{\int }_0^t\left({\omega }_{x}x-{\omega }_{y}y\right)e^{-\mathrm{\Gamma }\tau }d\tau \right]+{\omega }_{x}x-{\omega }_{y}y\\ \mathrm{\hspace{1em}\hspace{0.17em}}=\mathrm{\Gamma }{e}^{\mathrm{\Gamma }t}[f\left(0\right)+\rho \left\{H\left(t+\frac{k_1}{{\xi }_I}\right)-H\left(\frac{k_1}{{\xi }_I}\right)\right\}\\ \mathrm{\hspace{1em}\hspace{0.17em}}-\lambda \left\{H\left(t+\frac{k_2}{{\xi }_I}\right)-H\left(\frac{k_2}{{\xi }_I}\right)\right\}]\\ \mathrm{\hspace{1em}\hspace{0.17em}}+{\omega }_{x}x-{\omega }_{y}y\end{array}$$[3.22]

where we used Eq. [3.14a] and defined H(x) as   

$$H\left(x\right)=\frac{e^{({\xi }_R-\mathrm{\Gamma })x}\left(-{\xi }_{I}cos\left[{\xi }_{I}x\right]+({\xi }_R-\mathrm{\Gamma })\text{sin}\left[{\xi }_{I}x\right]\right)}{{\left({\xi }_R-\mathrm{\Gamma }\right)}^2+{\xi }_I{}^2}.$$[3.23]

We note in Eq. [3.22] that the multiplicative factor f(t) reflects the interacting effect of mutations during the transmission in a population (Appendix 2). Even though the two terms in the right-hand side of Eq. [3.19] are assumed to be additive, the differential equation mixes their contributions.

ANALYSIS OF ACTUAL DATA

In this section, we analyze some data based on Eq. [3.22] with Eq. [3.23]. We see that the formula has the generic form:   

$$\frac{dI}{dt}=c_1{e}^{\mathrm{\Gamma }t}+e^{{\xi }_{R}t}\left\{c_{2}sin\left({\xi }_{I}t\right)+c_3\text{cos}\left({\xi }_{I}t\right)\right\}$$[4.1]

where c_i (i = 1, 2, 3) are certain constants. We note that ξ_R and ξ_I may be linear functions of t as discussed in the previous model section. We see that the first term on the right-hand side represents the contribution of non-genetic factors such as vaccination, physical distancing and PCR testing. The second term represents the genetic effect. Here and subsequently, we use these terminologies as far as the generic form is concerned. It is to be noted that the generic form merely singles out the two seemingly independent factors. Because of this, when Eq. [4.1] is used for overall data fitting, it does not allow us to assess relative contributions of genetic and non-genetic factors.

Tokyo 5th wave

We analyze the Tokyo 5th wave data to see how well our new way of analysis works, especially in comparison with the standard epidemiology model in Eqs. [3.16]. The Tokyo data on daily infection numbers in the 5th wave period are given in https://www.google.com/search?q=tokyo+data+on+covid-19&oq=tokyo+data&aqs=chrome.1.69i57j35i39j0i512l3j0i30l5.9499j0j15&sourceid=chrome&ie=UTF-8.   

$$\begin{array}{l}\mathrm{D}\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{T}\mathrm{o}\mathrm{k}\mathrm{y}\mathrm{o}5\mathrm{t}\mathrm{h}\mathrm{w}\mathrm{a}\mathrm{v}\mathrm{e}=\{453,619,317,716,\\ 896,830,1008,1359,3177,2195,4498,4989,4377,\\ 4392,4227,3168,968,1273,763,253,299,200\}.\end{array}$$[4.2]

Data [4.2] show 22 daily cases in units of 5 days from June 18 to October 1, 2021, so that the entire period spans 110 days (Figs. 2, 3).

Fig. 2.

Tokyo daily infection numbers from June 18 to October 1, 2021 (ordinate) against days (t) elapsed since June 18, 2021 (abscissa).

Fig. 3.

Daily infection numbers in Tokyo (ordinate) against days (t) since June 18, 2021 (abscissa): (A) during the initial 40 days of the 5th wave expressed in Eq. [4.3], (B) during the middle 25 days expressed in Eq. [4.4], (C) during the later stage of the Tokyo 5th wave as expressed in Eq. [4.5], and (D) during the entire Tokyo 5th wave period fit by Eq. [4.9].

We fit the data with Eq. [4.1] using the nonlinear fitting of Mathematica. We found that we cannot fit the data for the entire period of the 5th wave if we retain the time-independent ξ_R. This is reasonable because the error threshold (ξ_R = 0) cannot be reached at the start of the spreading; it is achieved sometime during the spreading process. We use two strategies to fit the data:

(1) We split the data into three periods and fit the data in each period separately.

(2) We assume a simple form of time dependence for ξ_R. The advantage of this method is that we can calculate from the data when the error threshold is reached.

We start with the method (1). We divide the period of the Tokyo 5th wave into three periods: (A) first 40 days, (B) middle 25 days, and (C) later 45 days. The result of the fitting is as follows:

(A) period:   

$$\begin{array}{l}\frac{dI}{dt}=0.00018e^{0.406t}-133.7e^{0.039t}\text{sin}\left[0.016t\right]\\ \mathrm{\hspace{1em}\hspace{0.5em}}+419.2e^{0.039t}cos\left[0.016t\right].\end{array}$$[4.3]

The curve [4.3] is depicted in Fig. 3A. Note that ξ_R = 0.039 is positive, indicating that the error threshold has not yet been reached.

(B) period:   

$$\begin{array}{l}\frac{dI}{dt}=34.49e^{-0.044t}-23.68e^{0.0998t}sin\left[0.0196t\right]\\ \mathrm{\hspace{1em}\hspace{0.5em}}+87.13e^{0.0998t}cos\left[0.0196t\right].\end{array}$$[4.4]

Note that ξ_R = 0.0998 is even larger than the value in the initial period and still below the threshold (Fig. 3B).

(C) period:   

$$\begin{array}{l}\frac{dI}{dt}=3550e^{-0.0026t}+3171e^{-0.0093t}sin\left[0.056t\right]\\ \mathrm{\hspace{1em}\hspace{0.5em}}-5652e^{-0.0093t}cos\left[0.056t\right].\end{array}$$[4.5]

We see that ξ_R = −0.0093 is negative in this period, satisfying the error catastrophe condition. Eq. [4.5] is depicted in Fig. 3C.

We learned above that the value of ξ_R is dependent on which period we consider. This suggests that we may assume some simple form of time dependence for ξ_R and also for Γ in Eq. [4.1] (see the discussion following Eqs. [3.3] and [3.11], as well as Appendix 1). The form we assume is the following:   

$${\xi }_R={\xi }_{R0}\left(1-At\right)$$[4.6]

  

$$\mathrm{\Gamma }={\mathrm{\Gamma }}_0\left(1-Bt\right).$$[4.7]

We then fit the entire data with a single function:   

$$c_1{e}^{{\mathrm{\Gamma }}_0\left(1-Bt\right)t}+e^{{\xi }_{R0}\left(1-At\right)t}\left\{c_{2}sin\left[{\xi }_{I}t\right]+c_3\text{cos}\left[{\xi }_{I}t\right]\right\}.$$[4.8]

The result of the fitting is given as:   

$$\begin{array}{l}\frac{dI}{dt}=453e^{0.024\left(1-0.0099t\right)t}+3.77e^{0.216\left(1-0.0078t\right)t}sin\left[0.041t\right]\\ \mathrm{\hspace{1em}\hspace{0.5em}}-2.32e^{0.216\left(1-0.0078t\right)t}cos\left[0.041t\right].\end{array}$$[4.9]

We see from Eq. [4.9] that ξ_R changes from positive to negative during the 5th wave period. We can calculate the time of crossing the error threshold from:   

$${\xi }_R=0.216\left(1-0.0078t_t\right)=0$$[4.10]

namely, t_t = 128 days after the initiation of the Delta variant spreading. Although this value of t_t indicates the later decay phase of the spreading, we see from Eq. [4.9] that other factors and even the non-genetic term contribute to the precise determination of the maximum position of the infection numbers. Fitting of the entire period by Eq. [4.9] is shown in Fig. 3D. The non-genetic term takes the maximum value at t_n = $\frac{1}{2B}$ ≈ 51 and the genetic term does so at t_g = $\frac{1}{2A}$ ≈ 64, half of t_t = 128. The fitting describes clearly the battleground of virus vs. human immune system. Initially, the virus is winning, but eventually it reaches the maximum spreading power. The APOBEC3 deaminases involved in the human innate immune system (Robson et al., 2020; Simmonds, 2020; Ratcliff and Simmonds, 2021) continually force mutation of the viral RNA, which cannot avoid more and more errors so that the spreading power starts to shrink.

ANALYSIS OF DATA FROM SOUTH AFRICA AND JAPAN

We learned in the previous section that it is possible to fit the data of an entire period of a certain wave if we take into account the time dependence of ξ_R and Γ (infection rate β minus recovery rate γ in Eq. [3.20]). We adopt this strategy rather than dividing a period into several sub-periods.

The Omicron data were still insufficient when we were writing this paper, and the following analysis should be taken as preliminary rather than final. The Omicron variant was first found in South Africa. We therefore begin our analysis of the Omicron variant in South Africa. The data from November 29, 2021 to January 13, 2022 are given in https://www.worldometers.info/coronavirus/country/south-africa/.   

$$\begin{array}{l}\mathrm{D}\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{S}\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{h}\mathrm{A}\mathrm{f}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{a}=\{2273,4393,8561,15535,\\ 16055,16366,10000,6381,13147,19842,22383,19017,\\ 17153,37875,13288,20000,26389,24785,20713,\\ 16080,15465,8511,15423,21098,21156,13000,14824,\\ 5603,3778,6000,9020,12978,11754,9793,6000,\\ 3076,8078,11106,9858,9259,7759,4482,\\ 2409,4500,6760,5917\}.\end{array}$$[5.1]

The data of the entire 46 cases are depicted in Fig. 4A. We fit the data with a function that is exactly the same as the one in Eq. [4.8], which is used to explain the entire period of the Tokyo 5th wave. The result of the fitting gives:   

$$\begin{array}{l}\frac{dI}{dt}=78.71e^{0.38\left(1-0.0156t\right)t}+207.4e^{0.2015\left(1-0.0207t\right)t}sin\left[0.072t\right]\\ \mathrm{\hspace{1em}\hspace{0.5em}}+4548e^{0.2015\left(1-0.0207t\right)t}cos\left[0.072t\right].\end{array}$$[5.2]

Fig. 4.

(A) South Africa data for daily infections (ordinate) against days (t) starting from November 29, 2021 (abscissa). (B) South Africa Omicron data (the same as (A)) and the nonlinear fitting by Eq. [5.2]. (C) Extension of (B) to a much later time.

We see from Eq. [5.2] that   

$${\xi }_R=0.2015\left(1-0.0207t\right).$$[5.3]

Eq. [5.3] satisfies the error catastrophe condition of ξ_R < 0 after t_t = 48 days. This value of 48 days is much smaller than the Tokyo 5th wave (presumed to be the Delta type) value of t_t = 128 days, indicating that Omicron infection is more powerful and in turn it reaches the error threshold much faster than the Delta type. We show the curve described by Eq. [5.2] in Fig. 4B. We extend the nonlinear fitting in Fig. 4B to a later time as in Fig. 4C. We see from Figs. 4B, 4C that Omicron variant infection has a much steeper increase but also a sharp decrease due to the earlier arrival at the error threshold. Again, this transition happens at t_g = 24, which is earlier than t_t = 48. The non-genetic term takes the maximum value at t_n = 32 (Figs. 4B, 4C).

Universality of the parameters ξ_R and ξ_I

Before fitting the Japan data, we discuss the nature of the parameters ξ_R and ξ_I. As shown above, the time to reach the error threshold is calculated using the time dependence of ξ_R. We obtained the value t_t = 128 days for the Delta variant using the Tokyo data and t_t = 48 days for the Omicron variant using the South Africa data. We assume that these values depend on the variants we consider, and not or very little on any other external effect such as vaccination and physical distancing. The reason is that ξ_R and ξ_I are genetic parameters and not epidemic parameters. They are more or less inherent in the RNA of the virus and can also depend on the genetically determined human immune system. The other parameters in Eq. [4.1] can depend on vaccination, physical distancing, PCR testing followed by segregation, etc.. The above argument allows us to use the values of ξ_R and ξ_I obtained from the South Africa data to analyze the Omicron data from other regions, although there may be some difference in the values of ξ_R and ξ_I due to population structure, such as how much of particular genomic segments of Neanderthals has introgressed into the DNA of the group being dealt with (Zeberg and Pääbo, 2020a, 2020b, 2021). Likewise, it is likely that virus epidemics and host adaptive responses have been limited to local populations. In the case of coronavirus, interactions between virus and human genes started more than 20,000 years ago but have left genomic signatures only in East Asians (Souilmi et al., 2021). In the following analysis of the Japan data, we use the South Africa ξ_R and ξ_I for simplicity.

The validity or the usefulness of the universality of ξ_R and ξ_I can be checked by analysis using an index such as the AIC. Along the same lines, we can also consider (calculate the AIC of) the limiting case of all the genetic parts of the parameters vanishing. These investigations are left to our future research.

Japan’s case: the 6th wave

We analyze the Japan data for the daily infection number starting from December 30, 2021 to January 20, 2022 when the Omicron variant was starting to dominate. We use the same ξ_R and ξ_I obtained in the South Africa analysis. Adding the numbering of dates in the first column of the data list, we take the daily cases for the 6th wave of Japan from https://www.worldometers.info/coronavirus//country/japan/.   

$$\begin{array}{l}\mathrm{D}\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{f}\mathrm{o}\mathrm{r}6\mathrm{t}\mathrm{h}\mathrm{w}\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{o}\mathrm{f}\mathrm{J}\mathrm{a}\mathrm{p}\mathrm{a}\mathrm{n}=\{\left\{1,466\right\},\left\{2,490\right\},\\ \left\{3,463\right\},\left\{4,504\right\},\left\{5,600\right\},\left\{6,783\right\},\left\{7,1256\right\},\\ \left\{8,2506\right\},\left\{9,4194\right\},\left\{10,5983\right\},\left\{11,7930\right\},\\ \left\{12,8144\right\},\left\{13,6394\right\},\left\{14,9000\right\},\left\{15,12243\right\},\\ \left\{16,17940\right\},\left\{17,21241\right\},\left\{18,22707\right\},\left\{19,26881\right\},\\ \left\{20,22000\right\},\left\{21,29862\right\},\left\{22,39841\right\}\}.\end{array}$$[5.4]

We fit this with the following function, which is exactly the same as the one used in the South Africa analysis with an untested assumption that even the non-genetic parameter (B in Eq. [4.7]) is universal:   

$$\begin{array}{l}\frac{dI}{dt}=c_1{e}^{{\mathrm{\Gamma }}_0\left(1-0.0156\text{t}\right)\text{t}}+c_2{\text{e}}^{0.2015\left(1-0.0207t\right)t}sin\left[0.072t\right]\\ \mathrm{\hspace{1em}\hspace{0.5em}}+c_3{\text{e}}^{0.2015\left(1-0.0207t\right)t}cos\left[0.072t\right].\end{array}$$[5.5]

The result of the fitting gives   

$$\begin{array}{l}\frac{dI}{dt}=0.00072e^{1.158\left(1-0.0156t\right)t}+2116{\text{e}}^{0.2015\left(1-0.0207t\right)t}\text{sin}\left[0.072t\right]\\ \mathrm{\hspace{1em}\hspace{0.5em}}+496.1{\text{e}}^{0.2015\left(1-0.0207t\right)t}cos\left[0.072t\right].\end{array}$$[5.6]

Data [5.4] and Eq. [5.6] are plotted in Fig. 5A. From Fig. 5B, we see that the Japan’s daily infection number rises to approximately 100,000 and then starts decreasing. Our analysis of the Japan data from December 30, 2021 to January 20, 2022 in Data [5.4] correctly predicted the maximum daily infection number to be approximately 100,000. However, the actual speed of the decline from the peak number is much lower than predicted by the analysis. This is due to the invalidity of the untested assumption that even the non-genetic parameters are universal.

Fig. 5.

(A) Japan daily infection numbers (ordinate) and the nonlinear fitting by Eq. [5.6] against days (t) elapsed since December 30, 2021 (abscissa). (B) Extension of (A) to a later time.

Analysis of later data from Japan

We investigate the situation by analyzing the data from January 1 to March 31, 2022 (every 3 days). Here, we invalidate the untested assumption that the non-genetic parameters are universal. Using the same format as in the case of the 6th wave of Japan, we take the data for the infection number from https://www.worldometers.info/coronavirus//country/japan/. As in Data [5.4], the numbering of dates is presented in the first column of the following list, but here in units of 3 days.   

$$\begin{array}{l}\mathrm{L}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{o}\mathrm{f}\mathrm{J}\mathrm{a}\mathrm{p}\mathrm{a}\mathrm{n}=\{\left\{1,2505\right\},\left\{2,8144\right\},\left\{3,6829\right\},\\ \left\{4,21241\right\},\left\{5,23351\right\},\left\{6,44638\right\},\left\{7,47176\right\},\\ \left\{8,69736\right\},\left\{9,82159\right\},\left\{10,93388\right\}\left\{11,94431\right\},\\ \left\{12,7170\right\},\left\{13,10097\right\},\left\{14,80234\right\},\left\{15,95603\right\},\\ \left\{16,77135\right\},\left\{17,66373\right\},\left\{18,71254\right\},\left\{19,55243\right\},\\ \left\{20,70958\right\},\left\{21,55417\right\},\left\{22,62747\right\},\left\{23,52002\right\},\\ \left\{24,50781\right\},\left\{25,49210\right\},\left\{26,29657\right\},\left\{27,47937\right\},\\ \left\{28,32916\right\},\left\{29,53544\right\}\}.\end{array}$$[5.7]

We fit the data with the following function:   

$$\begin{array}{l}\frac{dI}{d\left(3\tau \right)}=c_1{\text{e}}^{{\mathrm{\Gamma }}_0\left(1-B\times 3\tau \right)3\tau }+c_2{\text{e}}^{0.2015\left(1-0.0207\times 3\tau \right)3\tau }\text{sin}\left[0.072\times 3\tau \right]\\ \mathrm{\hspace{1em}\hspace{1em}}+c_3{\text{e}}^{0.2015\left(1-0.0207\times 3\tau \right)3\tau }cos\left[0.072\times 3\tau \right].\end{array}$$[5.8]

The meaning of both sides in Eq. [5.8] is the following:

(1) Compared to the previous fitting curve given in Eq. [5.5], we see that we use 3τ rather than t to designate the time variable (here we have data points of every 3 days as in Data [5.7] and τ ranges from 1 to 29).

(2) Otherwise, the only difference is the insertion of the parameter B in the first term. This is the most important deviation from Eq. [5.5] where the untested assumption of the non-genetic term was made.

(3) The universality of the mutation contribution is maintained: the second and third terms are the same as in Eq. [5.5].

The result of the fitting gives:   

$$\begin{array}{l}\frac{dI}{d\left(3\tau \right)}=6882{\text{e}}^{0.2431\left(1-0.0269\tau \right)\tau }+2886{\text{e}}^{0.6046\left(1-0.0621\tau \right)\tau }\text{sin}\left[0.216\tau \right]\\ \mathrm{\hspace{1em}\hspace{1em}}-3106{\text{e}}^{0.6046\left(1-0.0621\tau \right)\tau }cos\left[0.216\tau \right].\end{array}$$[5.9]

This nonlinear fitting is shown in Fig. 6A. We extend the function in Eq. [5.9] to τ = 50 (t = 150 days after January 1, 2022, i.e., the end of May). It is depicted in Fig. 6B. The genetic and non-genetic terms take their respective maximum values at τ_g = 8.05 and τ_n = 18.6. In the daily unit, these correspond to t_g = 24 and t_n = 56, respectively. Similarly, the error threshold occurs at t_t = 48 days after January 1, the same as in the case of South Africa.

Fig. 6.

(A) Nonlinear fitting of the Japan data on daily infections by Eq. [5.9] (ordinate) against τ (in units of 3 days) from January 1 to March 31 (abscissa, corresponding to the numbering of dates in the first column in Data [5.7]). (B) Extension of (A) from early January to the end of May.

Tentative summary of the analysis of the later data from Japan

(1) As one can see by comparing Eqs. [5.6] and [5.9], the difference between the previous analysis and the present one is mostly in the first term representing the contribution of non-genetic factors such as the physical distancing effect. The previous best fitting resulted in minimizing the first term, whereas the current fitting has a much larger contribution from this term.

(2) The exponential factor of the first term in the current analysis [exp{0.2431(1−0.0269τ)τ}] in Eq. [5.9] behaves in approximately the same way as that in the previous analysis [exp{1.158(1−0.0156t)t}] in Eq. [5.6], indicating that the structure of non-genetic factors was correctly taken into account in both analyses. The difference is in the overall factor of the first term (0.00072 vs. 6882).

(3) As one can see from Fig. 6B, there will be no 7th wave in Japan if our analyses are correct. Unfortunately, however, it appears that before Omicron BA. 2 completely disappeared, it had gradually been replaced by a new Omicron variant strain, BA. 5, which caused the 7th wave in Japan in the summer of 2022.

ACKNOWLEDGMENTS

One of the authors (H. S.) would like to express his sincere gratitude to Professor Toshimichi Ikemura for his valuable discussions and guidance. We thank two anonymous reviewers for their comments and constructive criticisms, and Dr. Ian Smith for his editorial supervision.

REFERENCES

•  Bedford,  T.,  Rambaut,  A., and  Pascual,  M. (2012) Canalization of the evolutionary trajectory of the human influenza virus. BMC Biol. 10, 38.
•  Eigen,  M. (1971) Selforganization of matter and the evolution of biological macromolecules. Naturwissenschaften 58, 465–523.
•  Kermack,  W. O., and  McKendrick,  A. G. (1927) A contribution to the mathematical theory of epidemics. Proc. R. Soc. Lond. A 115, 700–721.
•  Kucharski,  A. J.,  Andreasen,  V., and  Gog,  J. R. (2016) Capturing the dynamics of pathogens with many strains. J. Math. Biol. 72, 1–24.
•  Miller,  J. K.,  Elenberg,  K., and  Dubrawski,  A. (2022) Forecasting emergence of COVID-19 variants of concern. PLoS One 17, e0264198.
•  Ratcliff,  J., and  Simmonds,  P. (2021) Potential APOBEC-mediated RNA editing of the genomes of SARS-CoV-2 and other coronaviruses and its impact on their longer term evolution. Virology 556, 62–72.
•  Robson,  F.,  Khan,  K. S.,  Le,  T. K.,  Paris,  C.,  Demirbag,  S.,  Barfuss,  P.,  Rocchi,  P., and  Ng,  W.-L. (2020) Coronavirus RNA proofreading: molecular basis and therapeutic targeting. Mol. Cell 80, 1136–1138.
•  Rohani,  P.,  Green,  C. J.,  Mantilla-Beniers,  N. B., and  Grenfell,  B. T. (2002) Ecological interference between fatal diseases. Nature 422, 885–888.
•  Sadeghpour,  S.,  Khodaee,  S.,  Rahnama,  M.,  Rahimi,  H., and  Ebrahimi,  D. (2021) Human APOBEC3 variations and viral infection. Viruses 13, 1366.
•  Simmonds,  P. (2020) Rampant C→U hypermutation in the genomes of SARS-CoV-2 and other coronaviruses: causes and consequences for their short- and long-term evolutionary trajectories. mSphere 5, e00408-20.
•  Smith,  E. C.,  Blanc,  H.,  Vignuzzi,  M., and  Denison,  M. R. (2013) Coronaviruses lacking exoribonuclease activity are susceptible to lethal mutagenesis: evidence for proofreading and potential therapeutics. PLoS Pathog. 9, e1003565.
•  Souilmi,  Y.,  Lauterbur,  M. E.,  Tobler,  R.,  Huber,  C. D.,  Johar,  A. S.,  Moradi,  S. V.,  Johnston,  W. A.,  Krogan,  N. J.,  Alexandrov,  K., and  Enard,  D. (2021) An ancient viral epidemic involving host coronavirus interacting genes more than 20,000 years ago in East Asia. Curr. Biol. 31, 3504–3514.e9.
•  Sugawara,  H. (2021a) On the effectiveness of the search and find method to suppress spread of SARS-CoV-2. Proc. Jpn. Acad. Ser. B Phys. Biol. Sci. 97, 22–49.
•  Sugawara,  H. (2021b) Effect of vaccination against COVID-19 spreading. Proc. Jpn. Acad. Ser. Phys Biol Sci. 97, 543–558.
•  Thakur,  V.,  Bhola,  S.,  Thakur,  P.,  Patel,  S. K. S.,  Kulshrestha,  S.,  Ratho,  R. K., and  Kumar,  P. (2022) Waves and variants of SARS-CoV-2: understanding the causes and effect of the COVID-19 catastrophe. Infection 50, 309–325.
•  Zeberg,  H., and  Pääbo,  S. (2020a) The major genetic risk factor for severe COVID-19 is inherited from Neanderthals. Nature 587, 610–612.
•  Zeberg,  H., and  Pääbo,  S. (2020b) The MERS-CoV receptor gene is among COVID-19 risk factors inherited from Neandertals. bioRxiv. doi: https://doi.org/10.1101/2020.12.11.422139
•  Zeberg,  H., and  Pääbo,  S. (2021) A genomic region associated with protection against severe COVID-19 is inherited from Neanderthals. Proc. Natl. Acad. Sci. USA 118, e2026309118.