Next Article in Journal
A Multiplex Quantitative Polymerase Chain Reaction for the Rapid Differential Detection of Subgroups A, B, J, and K Avian Leukosis Viruses
Previous Article in Journal
Correction: Magazine et al. Mutations and Evolution of the SARS-CoV-2 Spike Protein. Viruses 2022, 14, 640
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Mathematical Modelling of Virus Spreading in COVID-19

1
Faculty of Physical Science and Technology, Inner Mongolia University, 235 West College Road, Hohhot 010021, China
2
College of Science, Inner Mongolia University of Technology, 49 Aymin Street, Hohhot 010051, China
*
Authors to whom correspondence should be addressed.
Viruses 2023, 15(9), 1788; https://doi.org/10.3390/v15091788
Submission received: 3 August 2023 / Revised: 20 August 2023 / Accepted: 21 August 2023 / Published: 23 August 2023
(This article belongs to the Section SARS-CoV-2 and COVID-19)

Abstract

:
A mathematical model is proposed to analyze the spreading dynamics of COVID-19. By using the parameters of the model, namely the basic reproduction number (R0) and the attenuation constant (k), the daily number of infections (DNI) and the cumulative number of infections (CNI) over time (m) are deduced and shown to be in good agreement with experimental data. This model effectively addresses three key issues: (1) inferring the conditions under which virus infections die out for a specific strain given R0; (2) explaining the occurrence of second waves of infection and developing preventive measures; and (3) understanding the competitive spread of two viruses within a region and devising control strategies. The findings highlight the potential of this simple mathematical framework in comprehensively addressing these challenges. The theoretical insights derived from this model can guide the evaluation of infection wave severity and the formulation of effective strategies for controlling and mitigating epidemic outbreaks.

1. Introduction

The ongoing COVID-19 pandemic has posed several unresolved questions. Firstly, our aim is to understand the reason behind the waveform pattern observed in the daily number of infections (DNI) and to predict the scale and duration of an infection wave, including the cumulative number of infections (CNI) and the time until infection ends. Secondly, we investigate why the second wave of infection typically follows the first wave and how to predict its occurrence. Lastly, we explore how virus spreading in a region is influenced by the emergence of new strains and develop methods to control their spread. Existing epidemic models have limitations in addressing these problems comprehensively. Since the 1920s, differential equations were used to model the population distribution of disease spread, including susceptible, infected, and recovered/dead pools [1,2]. Such models could not examine important social and behavioural factors, such as the behavioural responses of individuals to policy measures, and the effect of heterogeneous social contacts on diffusion patterns. Next, since the 1990s, agent-based simulations were proposed that included some important sources of population heterogeneity and explored the structure and dynamics of transmission networks. However, none of the agent-based models are based on explicit empirical/theoretical assumptions of individual behaviour, social transmission mechanisms, and social structure constraints [3,4,5]. Recently, several network models simulating the spread of epidemics in the population were proposed, such as the Susceptible-Infectious-Susceptible model on random networks [6], and the epidemic spreading on modular networks [7]. Although these specific models have solved some aspects of the complexity of infectious diseases and obtained meaningful results, they still fail to comprehensively answer the aforementioned unresolved questions.
As the Chinese idiom goes, “the greatest truths are the simplest”. Starting from the principle of natural selection—the interaction and compromise between the virus and its host (human)—we propose a mathematical model based on insights from experimental data. The model is based on two fundamental assumptions. First, each viral strain has a basic reproduction number (R0), which represents the average number of infections caused by an initial infectious person in a completely susceptible population [8]. R0 is a measure of the transmission potential of a particular infectious disease. Second, virus infections undergo a series of events due to the presence of population immunity or intervention that makes the instantaneous reproduction number, Rt, decreasing step by step, described by an attenuation parameter k (k < 1) [9,10]. The parameter k is influenced by social contacts and policy measures, where looser contacts and measures result in higher values of k. It is worth noting that public health measures in a given region usually undergo significant changes over a longer period of time compared to the duration of a virus wave [11]. The stringency index of policy measures can be found at https://ourworldindata.org/covid-stringency-index (accessed on 13 August 2023). Therefore, it is reasonable to assume that the attenuation parameter k is approximately constant during a single wave of virus spread. Based on two parameters, R0 and k, we can deduce the general formula for the daily number of infections (DNI) and the cumulative number of infections (CNI) over m steps. The CNI is denoted by F(R0, k; m). These formulas are in good agreement with experiments and provide a generalized framework to simulate single and multiple virus infections. By introducing F1 and F2 we can also study the competitive spread of two viruses in a region, including scenarios with the introduction of a new virus strain. It is worth noting that the loosening of public health measures and/or the emergence of new viruses can lead to subsequent waves of virus transmission.

2. Materials and Methods

The data of daily virus infections of COVID-19 were taken from WHO (https://covid19.who.int/data, accessed on 12 May 2023) for the UK and from the public database (https://ourworldindata.org/coronavirus, accessed on 12 May 2023) for Hong Kong. Since Gauss discovered the least square method and successfully applied it to astronomical observation, the ordinary least square is recognized as one of the best methods for curvilinear regression. In the present study, we used the least square simulation of observational data of the virus infection to determine the parameters in each COVID-19 pandemic. The goodness of the nonlinear least square fitting (NLLSF) is tested by R2, i.e.,
R 2 = 1 i = 1 n ( y i π ^ i ) 2 i = 1 n y i y ¯ 2 ,
where yi are the observation values, y ¯ their average, and π ^ i the predictions of the model. The goodness R2 and the p-value of Prob(F-statistic) are calculated for each simulation. The F-statistic is defined as
F - s t a t i s t i c = ( i = 1 n y i y ¯ 2 i = 1 n y i π ^ i 2 ) / ( k 1 ) i = 1 n y i π ^ i 2 n k ,
(n–number of samples, ν–number of parameters) that can be calculated from the experimental data. The p-value is decided by the percentile of the Fisher’s F-distribution Fν−1,nν(z), namely,
p - v a l u e = t h e   p r o b a b i l i t y   o f   F ν 1 , n ν ( z > F - s t a t i s t i c )
The goodness R2, the p-value, and the root mean squared error (RMSE) for each simulation are given in the figures.

3. Results

3.1. Derivation of Formulas for Cumulative Number of Infections (CNI) and Daily Number of Infections (DNI)

The potential for infection of a given virus strain is determined by the basic reproduction number known as R0. For the virus to spread, it must continually re-infect the rest of the susceptible population. However, countermeasures by humans must aim to limit these subsequent infections. Essentially, this implies that the interaction and compromise between the virus and host (humans) necessitates the j-th reproduction number Rj to be smaller than the (j − 1)-th reproduction Rj−1, as expressed in the equation
R j = k R j 1 , ( k < 1 ) .
The reason Rj changes with j is twofold. Firstly, as the virus continues to spread, the susceptible population decreases, leading to a decrease in the transmission potential. Secondly, ongoing public health measures are being implemented to control the spread of the virus, further contributing to the reduction in the instantaneous reproduction number. The attenuation constant ‘k’ in this scenario is influenced by social contact rates and policy measures. After m rounds of infection spread, the total number of infections in a specific area can be represented as
F ( R 0 , k ; m ) = R 0 + R 0 R 1 + + R 0 R 1 R 2 R m 1 = R 0 + k R 0 2 + + k m ( m 1 ) / 2 R 0 m
This function F(R0, k; m) describes the cumulative number of infections (CNI), which is an increasing function of m. Assuming the total number of infected individuals to be N, solving the equation F(R0, k; m) = N will produce ‘m’, the transmission number of the virus strain. The parameter m grows with time t. Assuming each transmission takes place every q days, we have
m = ( t + t 0 ) / q
where t0 is a time shift parameter linked to the start of infection (the time of m = 1). Equations (1) and (2) provide the cumulative number of infections (CNI) for each branch of the epidemic. If various generations of different branches exist simultaneously, the total CNI is simply the sum of the CNIs of each branch. The CNI can be observed experimentally. Inserting (2) into (1), we obtain a formula of CNI vs. t represented by four parameters R0, k, q, and t0 that can be used to simulate the change of a cumulative number of infections with time in a region. As seen in Figure 1, our simulations of the UK-alpha strain (November 2020 to April 2021) and the Hong Kong-delta/omicron strain (February 2022 to April 2022) spreading are well fitted to the data obtained from COVID-19 pandemic updates. The high accuracy of the simulation shows that the assumption of the constancy of k is reasonable. It is important to note that our model not only allows for simulation of COVID-19 pandemics but also enables predictions of virus spread over longer periods by utilizing parameter data derived from a particular phase of the outbreak.
Applying Equations (1) and (2), the daily number of infections (DNI) is derived:
d F / d t = ( d F / d m ) / q ,   d F / d m = k m ( m 1 ) / 2 R 0 m .
From Equation (3), it is easy to understand why DNI first rises then falls, considering that R0 > 1 (for most infectious viruses) and k < 1. Zero DNI arrives when dF/dm = 0, which requires m to be large enough that
k ( m 1 ) / 2 R 0 < 1   o r   R 0 < ( 1 / k ) ( m 1 ) / 2
This is a condition for the end of an infection wave.

3.2. Insights from Typical Figures of CNI vs. m

Utilizing Equation (1), we plotted the change of CNIs with variables k and m for the given R0 (Figure 2). Figure 2A–E show typical cases that help in understanding the development and termination of an infection wave.
From Figure 2, one can observe that the curve of F(R0, k; m) increases with m and approaches a stable value when m > mst for each k < kth. We denote F(R0, kth; mst) as Nth. Calculating for mst = 15 in Figure 2A–E, we obtained the threshold kth and the corresponding CNI Nth for each given R0 (each virus strain). On the other hand, CNI attaining a stable value means the DNI is approaching zero. By defining Ekm = (1/kth)(mst−1)/2, Equation (4) can be written as R0 < Ekm. The threshold kth, the corresponding CNI Nth, and the parameter Ekm are listed in Table 1. We found the relation R0 < Ekm to hold well across all data.
The daily increasing number dF/dm = 0 means the virus infection is dying out. The above result shows the higher the basic reproduction number R0, the stricter the public health measure is required to be to increase Ekm to satisfy Equation (4) and terminate the wave of virus infection effectively. From Table 1, we found the cumulative infection numbers, Nth, of many virus strains (except Omicron with R0 = 18.6) are lower than 105 if an appropriate k (lower than kth) is introduced.
For the virus strain of given R0, if F(R0, k; m) (for all m) exceeds a threshold Nmax, then the spread of this strain would lead to a wave of COVID-19 infection in a region with a population larger than Nmax. In order to end this spread, the necessary condition would be
F ( R 0 ,   k ;   m ) < N m a x ( f o r   a l l   m )
Equation (5) provides a constraint on k for strains with R0. The critical value of k is denoted as kcr. Taking Nmax = 107, the values of kcr are also listed in Table 1. Note that the parameter kth is the threshold value of k required for ending the spread in the 15th generation, but the parameter kcr is that value for ending the spread in an arbitrary generation. The latter constraint is looser than the former. Therefore, kcr is larger than kth.
Virus strains with low R0 values have higher kcr and kth values. As a result, they can spread in regions with smaller populations under looser public health measures. This theoretical prediction can explain why strains like SARS-CoV in 2003 only spread in restricted regions and soon disappeared globally. On the contrary, virus strains with high R0 values, such as Omicron, have lower kcr and kth values. To satisfy Equations (4) and (5), the parameter k needs to meet very strict constraints. In cases where social management measures are not stringent enough, the number of infected individuals will quickly surpass the Nmax limit, triggering a global pandemic wave, possibly culminating in a type of coexistence between viruses and humans.
In summary, the termination of an infection wave is determined by the condition dF/dm = 0 on the CNI diagram (Figure 2), which requires a sufficiently large m and the condition in Equation (4) satisfied, i.e., R0 < Ekm = (1/kth)7. Additionally, the necessary condition for a pandemic to die out in a population of Nmax is expressed by Equation (5), imposing a limitation on k, namely k < kcr.

3.3. Prediction on the Second Wave of Pandemics

Early models based on previous pandemics such as SARS, MERS, and the 2009 H1N1 outbreak can effectively predict the occurrence of the first wave of a disease. However, their predictive power decreases when it comes to anticipating the possibility of a second wave [12]. This raises the question: why does the second wave of COVID-19 infections often follow the first wave? The present model aims to address this issue.
As shown in Figure 2, there are multiple curves (F versus m) for a given R0 value. These curves differ from each other by the parameter k. When the number of daily new infections (DNI) approaches zero, any fluctuation in k significantly influences the spread of the virus. Therefore, changes in public health measures can induce variations in virus spread along different curves. Generally, as public health measures are relaxed, the parameter k increases, causing the virus to transition from one curve of F to another steeper curve. This signifies the onset of a new wave of virus spread. For instance, in Figure 2A, the spread of the Omicron variant (with an R0 of 18.6) is plotted. Assuming the initial spread is along the curve with k equal to 0.613, a stable state is reached at m = 15. At this point, the k value increases to 0.722. In response to this change, the virus begins spreading along a new curve with k = 0.722, starting from m = 4, as the value of F(18.6, 0.722; m = 4) is equal to the original CNI value, F(18.6, 0.613; m = 15). This example explains how the second wave of virus infection occurs. However, in cases where the parameter k decreases when the DNI approaches zero, the curve of F will transition to a flatter one. This indicates that the first wave of viral infection will end soon, and no second wave will occur.
By examining Figure 2A–E, we can analyze how the occurrence of a second wave depends on the R0 value of the virus. For instance, as k changes from 0.85 to 0.9, the CNI (at m = 15) for high R0 viruses increases hundreds of times, whereas it only increases tens of times for low R0 viruses. Therefore, our model predicts that multiwave infections are more likely to occur with viruses that have high R0 values.
In the aforementioned discussions, we have assumed that no virus mutation occurs and that only one type of virus is spreading. In reality, changes in public health measures may be accompanied by virus mutations. In this case, the change in public health measures would cause the jump of F not only between different curves with a given R0 but also between different R0 values, providing more opportunities for the occurrence of a second wave during virus spread.
Another crucial point to consider is that the change in k can simultaneously induce a change in q, as the eigen-time (inherent time) m depends on q (Equation (2)). In the case of a single wave, the parameter m can be used to represent time dependence, and q can be simply set to 1 (known as normalization). However, when studying two continuous waves, the dependence on q should be clearly indicated due to the different eigen-times m. We have mentioned that the dependence of F on the change in k results in a jump from one curve to another. Meanwhile, the dependence of F on the change in q only affects the lengthening or shortening of the abscissa of the graph without altering the shape of the curve. When the public health measures change and a jump between curves occurs, the abscissa of the graph of the second wave simultaneously lengthens or shortens.
The occurrence of the second wave of infection is more likely when public health measures are relaxed. This change in the second wave is accompanied by a change in the q value. These predictions align with experimental data. For example, in the UK from May to September 2021, public health measures were relaxed, and a second wave of infection followed the first wave (Figure 1A). Similarly, in Hong Kong several months after May 2022, the looser public health measures led to a second wave occurring after the first wave (Figure 1B). (The data on the change of public health measures can be found at https://www.bsg.ox.ac.uk/research/covid-19-government-response-tracker, accessed on 13 August 2023).
In summary, the dependence of the CNI on k (given a specific R0) is determined by the jump between different curves on the graph F(R0, k; m) versus m, referred to as k-transformation. The dependence of the CNI on the duration of m is obtained by stretching the abscissa m on the graph, known as q-transformation. The k-transformation and q-transformation, occurring when the first wave is nearing its end, are the causes of a continuous second wave. The change in k is attributed to the modification of public health measures, while the modification of q is due to the change in physical time within a unit of m. The changes in k and q provide an explanation for the experimental data on continuous multiwave infections.

3.4. Cross-Spread of Two Viruses: Discriminant Function

Viral infections often involve the simultaneous presence of two or more viruses. To accurately simulate the cross-spread of two viruses, it is necessary to consider the differences in eigen-times (m1 and m2) between these viruses and their relationship with the physical time t. Consequently, additional parameters q and t0 (as given in Equation (2)) must be taken into account. By utilizing the four parameters R0, k, q, and t0, the CNI F(R0, k; m) can be expressed as:
F ( R 0 ,   k ;   m ( t ) ) = F ( t ;   a ,   b ,   c ,   d ) ,   ( a = R 0 ,   b = k ,   c = q ,   d = t 0 )
When the CNIs of two virus strains F1(t; a1, b1, c1, d1) and F2(t; a2, b2, c2, d2) intersect at tcr,
F 1 ( t c r ;   a 1 ,   b 1 ,   c 1 ,   d 1 ) = F 2 ( t c r ;   a 2 ,   b 2 ,   c 2 ,   d 2 )
and
F 1 ( t ;   a 1 ,   b 1 ,   c 1 ,   d 1 ) > F 2 ( t ;   a 2 ,   b 2 ,   c 2 ,   d 2 )   a s   t < t c r , F 1 ( t ;   a 1 ,   b 1 ,   c 1 ,   d 1 ) < F 2 ( t ;   a 2 ,   b 2 ,   c 2 ,   d 2 )   a s   t > t c r .
This represents a transition in the population of virus strains from F1 to F2. As it is challenging to directly solve the intersection equation (Equation (7)), we introduce a function:
D 21 = l o g [ ( d F 2 / d t ) / ( d F 1 / d t ) ]
which satisfies:
D 21 > 0   as   d F 2 / d t > d F 1 / d t , D 21 < 0   as   d F 2 / d t < d F 1 / d t .
By utilizing Equations (2), (3), (6) and (8), D21 can be formulated as a simple quadratic function of time
D 21 = α t 2 + β t + γ
where
α = log b 2 / 2 c 2 2 log b 1 / 2 c 1 2 ,   β = ( 2 d 2 c 2 ) log b 2 / ( 2 c 2 2 ) ( 2 d 1 c 1 ) log b 1 / ( 2 c 1 2 ) + log a 2 / c 2 log a 1 / c 1 ,   γ = d 2 ( d 2 c 2 ) log b 2 / ( 2 c 2 2 ) d 1 ( d 1 c 1 ) log b 1 / ( 2 c 1 2 ) + d 2 log a 2 / c 2 d 1 log a 1 / c 1 + log ( c 1 / c 2 ) .
In order to determine if a real root of D21 exists, we define:
= β 2 4 α γ
The value of Δ determines the existence of the real root in the quadratic form D21. This form, known as the discriminant function, serves as a tool to identify the occurrence of tcr and the domain of its existence. The prediction rules can be summarized as follows:
Rule 1: When Δ > 0, the quadratic form intersects with the t axis at ts and tm (tm > ts),
t s = β 2 α   ( for   α > 0 ) β + 2 α   ( for   α < 0 ) ,   t m = β + 2 α   ( for   α > 0 ) β 2 α   ( for   α < 0 ) .
These values partition the time into three distinct domains: t < ts in the first domain, ts < t < tm in the second domain, and t > tm in the third domain.
Rule 2: In the first domain, ta = qm0t0 is defined as the initial time where m0 >> 1, indicating m1,2 = (t + d1,2)/c1,2 >> 1. If D21 > 0, there will be no intersection of F1(t) and F2(t) in the domain between ta and ts when the initial values of F at ta satisfy F1(ta) < F2(ta), but there may be one tcr (the number of tcr is either 1 or 0) in the domain when F1(ta) > F2(ta). If D21 < 0, there will be no intersection of F1(t) and F2(t) in the domain when the initial values satisfy F1(ta) > F2(ta), but there may be one tcr (the number of tcr is either 1 or 0) when F1(ta) < F2(ta).
Rule 3: In the second domain, if D21 > 0 there will be no intersection of F1(t) and F2(t) in the domain when the F-values at t = ts satisfy F1(ts) < F2(ts), but there may be one tcr (the number of tcr is either 1 or 0) when F1(ts) > F2(ts). If D21 < 0 there will be no intersection of F1(t) and F2(t) in the domain when the F-values at t = ts satisfy F1(ts) > F2(ts), but there may be one tcr (i.e., the number of tcr is either 1 or 0) when F1(ts) < F2(ts). The F-values at t = ts are determined by F1(t) and F2(t) in the first domain.
Rule 4: In the third domain, if D21 > 0 there will be no intersection of F1(t) and F2(t) in the domain when the F-values at t = tm satisfy F1(tm) < F2(tm), but there may be one tcr (the number of tcr is either 1 or 0) when F1(tm) > F2(tm). If D21 < 0 there will be no intersection of F1(t) and F2(t) in the domain when the F-values at t = tm satisfy F1(tm) > F2(tm), but there may be one tcr (i.e., the number of tcr is either 1 or 0) when F1(tm) < F2(tm). The F-values at t = tm are determined by F1(t) and F2(t) in the second domain.
Rule 5: There can be at most one tcr in a given domain because the symbol of D21 is definite in any domain. The magnitude of the domain is an important factor to predict the occurrence of a tcr. For example, in the second domain the magnitude is tm − ts = Δ1/2/(2|α|), and the necessary condition for the occurrence of tcr is a large enough (tmts) or Δ.
Rule 6: When Δ < 0, the quadratic form D21 does not intersect with t axis. In this case, there is only one domain, and the rule is the same as that in the first domain given by Rule 2.
Figure 3 presents examples of cross-spread of two virus strains, 1 and 2. The left panel shows the discriminant function, and the right panel displays the cross-spread of the two strains. The influence of the change in parameters (as an example, we only assume the R0 value of strain 1 changes) on the intersection of two strains is shown in the figure. In Figure 3A, there is no intersection. In Figure 3B,C, there are two intersections in the second and third domain, respectively. In Figure 3D, there is one intersection in the second domain. These intersection occurrences are in agreement with the aforementioned prediction rules.
In the first domain we have introduced ta = qm0t0 as the initial time. To study the cross-spread in the time interval between m = 1 and ta where dF/dt in D21 is difficult to be defined, one should use the F(t)-ladder method as follows:
The step size of CNIs of the first few steps in ladders 1 and 2 are given by a1, b1a12, b13a13, b16a14, b110a15, b115a16, b121a17, b128a18, b136a19, etc., and a2, b2a22, b23a23, b26a24, b210a25, b215a26, b221a27, b228a28, b236a29, etc., respectively (Equation (1)). Strain 1 spreads from the 1st step a1 at t = c1d1 on ladder 1 and strain 2 spreads from the 1st step a2 at t = c2d2 on ladder 2. The CNI and arrival time t of the two strains on their respective ladders are listed in Table 2. For example, let us take c1 = 7, d1 = 15, c2 = 5, d2 = 10. The earlier arrival times, in turn, are t = c1d1 = −8, c2d2 = −5, 2c1d1 = −1, 2c2d2 = 0, 3c2d2 = 5, 3c1d1 = 6, etc. By using Table 2, one can easily calculate the CNI at each arrival time as the parameters a1, b1, a2, b2 are given, compare the CNI values on two ladders and obtain the information on the cross-spread of two strains.

3.5. Examples of the Cross-Spread of Two Viruses

The epidemics that occurred in the UK from November 2020 to February 2022 provide a clear demonstration of the cross-spread phenomenon involving multiple viruses. This process can be divided into five stages, namely: (A) the alpha epidemic, (B) the delta invasion and cross-spread of two strains, (C) the delta dominant stage, (D) the omicron invasion and cross-spread of delta and omicron, and (E) the omicron dominant stage. In this section, we will specifically focus on studying the cross-spread of two strains during stages B and D.
Firstly, let us examine the cross-spread between the alpha strain (designated as strain (1) and the delta strain (designated as strain (2) during stage B. The spread of the alpha strain during stage A has been represented in Figure 1A, where we obtained the parameter a1 = R0(1) = 3.9. The spread of the delta strain during stage C has been illustrated in Figure 4A, and we derived the parameter a2 = R0(2) = 5.1 based on this data. By using R0(1) and R0(2) as inputs, we simulated the cross-spread of the two strains, as depicted in Figure 4B. Furthermore, utilizing all the parameters, ai, bi, ci, di (i = 1,2), obtained from Figure 1A and Figure 4A,B, we constructed the discriminant function of the cross-spread and plotted the intersection of the two strains during stage B in Figure 4C. Interestingly, we discovered that one tcr occurs at t = 71 in the region where t > tm (tm = 44.5), which is consistent with the prediction outlined in Rule 4.
Similarly, we investigated the cross-spread between the delta strain (strain 1) and the omicron strain (strain 2) during stage D. The results of this analysis are shown in Figure 5A–C. From Figure 5C, we observed that a tcr appears at t = 42 in the region where t > tm (tm = 28.5), which is in agreement with the prediction made by Rule 4.
In both simulations, we assumed April 1 was the initial time for stage B, and November 11 was the initial time for stage D. Due to the uncertainty surrounding the exact timing of the delta invasion in stage B and the omicron invasion in stage D, we shifted the initial times t of stages B and D by several days. Remarkably, we found that the same values of k1, k2, q1, q2, t0(1) were obtained, and only t0(2) varied while still maintaining the invariant t + t0(2).

4. Discussion

4.1. On the Simulation of COVID-19 Cases

The traditional SIR-type epidemic models depict the exponential growth of the number of infected individuals. However, empirical data have demonstrated that COVID-19 outbreaks do not exhibit exponential growth, but rather follow a three-parameter Gompertz growth function [13,14]. A new compartment model, known as the broken link model, has been proposed in the literature to explain the mechanism of Gompertz growth [15]. However, our proposed model is logically simple. In this model, we suggest that the spread of the virus depends on four parameters: R0, k, q, and t0. R0 describes the inherent infectious ability of the virus, k represents the strictness of social management, q represents the time needed for one step of infection, and t0 is an additional parameter that aligns with the starting date of the experimental data. By incorporating these four parameters, our four-parameter simulation accurately fits all existing COVID-19 epidemic data and will contribute to the prediction of future outbreaks. Moreover, the mechanism of Gompertz growth has been elucidated by our model.

4.2. On the Mutation of the SARS-CoV-2 Virus

The SARS-CoV-2 virus continuously undergoes mutations, giving rise to new strains. This ongoing mutation process is the reason why the pandemic has persisted for more than three years. As a result of natural selection, these new mutant strains possess higher infectious ability but lower lethality rates. Although the lethality rate of the new strains may be lower, it still results in a significant number of deaths. Therefore, since mutation occurrence is inevitable and costly, the key focus should be on reducing the epidemic probability of mutants. While humans cannot prevent the virus from mutating, they can hinder the mutant strains from dominating the competition. According to six prediction rules on the cross-spreading of two strains, it is highly unlikely for the intersection point tcr to manifest under the following conditions: (1) when Δ > 0, there are three domains and in this case there will be no occurrence of tcr within a domain if the period of this domain is short enough (i.e., tsta small, Δ small, or tm large) or if the symbol of D21 within a domain does not align with the symbol of F1F2 at the initial time of this domain; (2) when Δ < 0, there is only one domain and in this case there will be no tcr if the period of virus spread is short enough or if the symbol of D21 does not align with the symbol of F1F2 at the initial time.
An example to highlight the phenomenon of no occurrence of tcr is the winter epidemic in China in 2022. Among tens of millions of SARS-CoV-2 cases in one city, no new mutant strain emerged or triumphed over the competition, with the exception of the original omicron strand. This suggests that the pandemic was effectively terminated within a short period, lasting only one month.
In summary, the strategies proposed by our model to control an epidemic have two main aspects. First, it is critical to prevent the occurrence of a second wave. Second, one should aim to avoid competition among mutant strains. In the latter scenario, minimizing the duration of virus spread is an efficient approach.

4.3. Dependence of Virus Infection Potential on Temperature and Humidity

The conformational equilibrium between the open and closed conformations of the receptor binding domain (RBD) of the spike (S) protein can be analyzed using first-principle techniques for a susceptible individual. The RBD can exist in either the open or closed position, known as the up or down conformation, respectively. The population of conformational states can be determined based on the free energy change during conformational transitions of the S protein. Let us denote the closed conformation as state A and the open conformation as state B. The Gibbs free energies of states A and B are represented by GA and GB, respectively. Generally, if GA is lower than GB, the RBD will primarily assume the inactive conformation A. Conversely, in order to initiate the infectious process, the equilibrium should shift towards the open conformation, implying that GB should be lower than GA [16]. The free energy G in a given conformation is related to the partition function Z, which can be expressed as follows:
G = k B T l n Z , Z = n e β E n ( β = 1 k B T ) .
Here En represents the energy level of a given conformation,
( E n ) A , B = V A , B + ( n + 1 / 2 ) ω A , B
VA,B represents two minima of conformational potential, respectively, and (n + 1/2)ħωA,B represents the corresponding vibrational energy around the minimum of the conformational potential. By the summation of Boltzmann factor over vibration states one has
Z A Z B = e β ( V A V B ) Y A / B , Y A / B = e 1 2 β ω B e 1 2 β ω B e 1 2 β ω A e 1 2 β ω A
Let TC represent the phase transition temperature, which can be determined by ΔG = GBGA = 0. From Equations (12) and (14) we obtain a simplified equation for TC
2 ( V B V A ) ( ω A ω B ) = c o t h ω A 2 k B T C ( a s ω A ω B ω A 1 )
and ΔG > 0 for T > TC, ΔG < 0 for T < TC. Therefore, if the environmental temperature decreases to T < TC, the conformational transition from the closed conformation to the open conformation occurs rapidly. This explains why the virus has higher transmission rates during the winter and the entrance to host cells is prioritized. Conversely, the summer season provides the most favorable conditions for virus elimination. In addition to temperature, the conformational equilibrium is influenced by humidity. The virus can be modeled as a charged sphere, and through the application of electrostatics principles to salty solutions, one can derive an expression for the potential at the surface of the charged sphere, where the dielectric constant ε is incorporated [17]. This implies that the elastic frequency ω2 should be replaced by ω2/ε. Considering water has a dielectric constant of ε = 80, the frequency parameter takes a reduced value of ω/9 in the presence of a fully salty solution rather than ω in a vacuum. Consequently, this provides a quantitative estimation of the virus’s infection potential, which is strongly dependent on humidity.
The above discussions on the susceptible individual can be extended to the population level, providing evidence that R0 is influenced by temperature and humidity. The point will be discussed in detail in the future work.

5. Conclusions

(1)
We propose a logically simple model for analyzing the spread of the COVID-19 virus. The model is based on two fundamental assumptions: each viral strain possesses a basic reproduction number R0, which quantifies its transmission potential, and the instantaneous reproduction number (Rt or Rj) decreases gradually due to factors such as population immunity and interventions, which can be represented by an attenuation constant k (k < 1).
(2)
The daily number of infections (DNI) and the cumulative number of infections (CNI) versus time (m) are deduced based on the aforementioned two assumptions. By utilizing the explicit relation m(t) where t is the physical time, our simulations with the four parameters (R0, k, q, t0) demonstrate excellent agreement with all experimental data.
(3)
Insights obtained from typical plots of CNI vs. m (i.e., typical plots of F(R0, k; m)) provide valuable information regarding the conditions required for the decline of a viral infection wave, as well as an explanation for the occurrence of continuous second waves of infection and its preventive measures.
(4)
The persistence of the SARS-CoV-2 pandemic for more than three years can be attributed to frequent mutations. We thoroughly examine the cross-spread of two strains within a region and lay a theoretical foundation for designing strategies to avoid competition among mutant strains or hindering the mutant strains from dominating the competition.

Author Contributions

Conceptualization, L.L.; validation, J.L. and L.L.; investigation, J.L.; writing—original draft preparation, L.L.; writing—review and editing, J.L.; visualization, J.L.; supervision, L.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Kermack, W.O.; McKendrick, A.G. A contribution to the mathematical theory of epidemics. Proc. R. Soc. Lond. A 1927, 115, 700–721. [Google Scholar] [CrossRef]
  2. Murray, J.D. Mathematical Biology, 2nd ed.; Springer-Verlag GmbH: Berlin/Heidelberg, Germany; New York, NY, USA, 1993; pp. 1–92. [Google Scholar] [CrossRef]
  3. Squazzoni, F.; Polhill, J.G.; Edmonds, B.; Ahrweiler, P.; Antosz, P.; Scholz, G.; Chappin, E.; Borit, M.; Verhagen, H.; Giardini, F.; et al. Computational models that matter during a global pandemic outbreak: A call to action. JASSS 2020, 23, 10. [Google Scholar] [CrossRef]
  4. Adam, D. Special report: The simulations driving the world’s response to COVID-19. Nature 2020, 580, 316–318. [Google Scholar] [CrossRef] [PubMed]
  5. Sridhar, D.; Majumder, M.S. Modelling the pandemic. BMJ 2020, 369, m1567. [Google Scholar] [CrossRef] [PubMed]
  6. Parshani, R.; Carmi, S.; Havlin, S. Epidemic threshold for the susceptible-infectious-susceptible model on random networks. Phys. Rev. Lett. 2010, 104, 258701. [Google Scholar] [CrossRef] [PubMed]
  7. Valdez, L.D.; Braunstein, L.A.; Havlin, S. Epidemic spreading on modular networks: The fear to declare a pandemic. Phys. Rev. E 2020, 101, 032309. [Google Scholar] [CrossRef] [PubMed]
  8. Wan, S.; Liu, J.; Liu, M. Progress on the basic reproduction number of SARS-CoV-2. Chin. Sci. Bull. 2020, 65, 2334–2341. [Google Scholar] [CrossRef]
  9. Bugalia, S.; Tripathi, J.P.; Wang, H. Estimating the time-dependent effective reproduction number and vaccination rate for COVID-19 in the USA and India. Math. Biosci. Eng. 2023, 20, 4673–4689. [Google Scholar] [CrossRef]
  10. Hasan, A.; Susanto, H.; Tjahjono, V.; Kusdiantara, R.; Putri, E.; Nuraini, N.; Hadisoemarto, P. A new estimation method for COVID-19 time-varying reproduction number using active cases. Sci. Rep. 2022, 12, 6675. [Google Scholar] [CrossRef]
  11. Hale, T.; Angrist, N.; Goldszmidt, R.; Kira, B.; Petherick, A.; Phillips, T.; Webster, S.; Cameron-Blake, E.; Hallas, L.; Majumdar, S.; et al. A global panel database of pandemic policies (Oxford COVID-19 Government Response Tracker). Nat. Hum. Behav. 2021, 5, 529–538. [Google Scholar] [CrossRef]
  12. Costris-Vas, C.; Schwartz, E.J.; Smith, R. Predicting COVID-19 using past pandemics as a guide: How reliable were mathematical models then, and how reliable will they be now? Math. Biosci. Eng. 2020, 17, 7502–7518. [Google Scholar] [CrossRef]
  13. Levitt, M.; Scaiewicz, A.; Zonta, F. Predicting the trajectory of any COVID19 epidemic from the best straight line. medRxiv 2020, 2020.06.26.20140814. [Google Scholar] [CrossRef]
  14. Ohnishi, A.; Namekawa, Y.; Fukui, T. Universality in COVID-19 spread in view of the Gompertz function. Prog. Theor. Exp. Phys. 2020, 12, 123J01. [Google Scholar] [CrossRef]
  15. Ikeda, Y.; Sasaki, K.; Nakano, T. A new compartment model of COVID-19 transmission: The broken-link model. Int. J. Environ. Res. Public. Health 2022, 19, 6864. [Google Scholar] [CrossRef] [PubMed]
  16. Luo, L.F.; Zuo, Y.C. Spike conformation transition in SARS-CoV-2 infection. arXiv 2021, arXiv:2009.11288. [Google Scholar]
  17. Phillips, R.; Kondev, J.; Theriot, J.; Garcia, H.G. Physical Biology of the Cell, 2nd ed.; Garland Science: New York, NY, USA, 2012; pp. 355–382. [Google Scholar]
Figure 1. COVID-19 pandemics in the UK (A) and Hong Kong (B).
Figure 1. COVID-19 pandemics in the UK (A) and Hong Kong (B).
Viruses 15 01788 g001
Figure 2. CNI functions F(R0, k; m) for several typical R0 values. (A): R0 = 18.6 (for Omicron BA.4, BA.5 in South Africa, 2022-1), (B): R0 = 9.5 (for Omicron B.1.1.529 in many countries, 2021-11), (C): R0 = 6.5 (for Delta in India, R0 = 5–8, 2020-10), (D): R0 = 4.5 (for Alpha in the UK, 2020-9, Beta in South Africa, 2020-5, Gamma in Brazil, 2020-11, R0 = 4–5), and (E): R0 = 2 (for SARS-CoV 2003, R0 = 2–3).
Figure 2. CNI functions F(R0, k; m) for several typical R0 values. (A): R0 = 18.6 (for Omicron BA.4, BA.5 in South Africa, 2022-1), (B): R0 = 9.5 (for Omicron B.1.1.529 in many countries, 2021-11), (C): R0 = 6.5 (for Delta in India, R0 = 5–8, 2020-10), (D): R0 = 4.5 (for Alpha in the UK, 2020-9, Beta in South Africa, 2020-5, Gamma in Brazil, 2020-11, R0 = 4–5), and (E): R0 = 2 (for SARS-CoV 2003, R0 = 2–3).
Viruses 15 01788 g002
Figure 3. Discriminant Function and intersection occurrence.
Figure 3. Discriminant Function and intersection occurrence.
Viruses 15 01788 g003
Figure 4. CNI simulation in the UK from November 2020 to November 2021. (A): CNI simulation of delta spread in stage C, (B): CNI simulation of alpha/delta spread in stage B, and (C): discriminant function and intersection of alpha/delta spread in stage B (left panel gives discriminant function and right panel the cross-spread of two strains). Note: CNI simulation of alpha spread in stage A from November 2020 to April 2021 has been plotted in Figure 1A.
Figure 4. CNI simulation in the UK from November 2020 to November 2021. (A): CNI simulation of delta spread in stage C, (B): CNI simulation of alpha/delta spread in stage B, and (C): discriminant function and intersection of alpha/delta spread in stage B (left panel gives discriminant function and right panel the cross-spread of two strains). Note: CNI simulation of alpha spread in stage A from November 2020 to April 2021 has been plotted in Figure 1A.
Viruses 15 01788 g004
Figure 5. CNI simulation in the UK from November 2021 to February 2022. (A): CNI simulation of omicron spread in stage E (only the data of first omicron peak are used.); (B): CNI simulation of delta/omicron spread in stage D; and (C): Discriminant function and intersection of delta/omicron spread in stage D (left panel gives discriminant function and right panel the cross-spread of two strains).
Figure 5. CNI simulation in the UK from November 2021 to February 2022. (A): CNI simulation of omicron spread in stage E (only the data of first omicron peak are used.); (B): CNI simulation of delta/omicron spread in stage D; and (C): Discriminant function and intersection of delta/omicron spread in stage D (left panel gives discriminant function and right panel the cross-spread of two strains).
Viruses 15 01788 g005
Table 1. Parameters related to virus infection dying out.
Table 1. Parameters related to virus infection dying out.
R0kth (mst = 15)Ekm = (1/kth)7Nthkcr (Nmax = 107)
18.60.61330.71050.722
9.50.7229.78∼1050.8
6.50.757.49∼1040.85
4.50.84.77∼1030.9
20.92.09∼102>0.95
Table 2. CNI and arrival time t of two strains on F(t)-ladder.
Table 2. CNI and arrival time t of two strains on F(t)-ladder.
tCNI
ci − di1ai
2ci − diai + biai2
3ci − diai + biai2 + bi3ai3
4ci − diai + biai2 + bi3ai3 + bi6ai4
5ci − diai + biai2 + bi3ai3 + bi6ai4 + bi10ai5
6ci − diai + biai2 + bi3ai3 + bi6ai4 + bi10ai5 + bi15ai6
1 i = 1, 2 refer to two strains, respectively, to save space only six steps of the F(t)-ladder are listed.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Luo, L.; Lv, J. Mathematical Modelling of Virus Spreading in COVID-19. Viruses 2023, 15, 1788. https://doi.org/10.3390/v15091788

AMA Style

Luo L, Lv J. Mathematical Modelling of Virus Spreading in COVID-19. Viruses. 2023; 15(9):1788. https://doi.org/10.3390/v15091788

Chicago/Turabian Style

Luo, Liaofu, and Jun Lv. 2023. "Mathematical Modelling of Virus Spreading in COVID-19" Viruses 15, no. 9: 1788. https://doi.org/10.3390/v15091788

APA Style

Luo, L., & Lv, J. (2023). Mathematical Modelling of Virus Spreading in COVID-19. Viruses, 15(9), 1788. https://doi.org/10.3390/v15091788

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop