key: cord-1006614-38xcv2se
authors: Djordjevic, Magdalena; Djordjevic, Marko; Ilic, Bojana; Stojku, Stefan; Salom, Igor
title: Understanding Infection Progression under Strong Control Measures through Universal COVID‐19 Growth Signatures
date: 2021-03-01
journal: Glob Chall
DOI: 10.1002/gch2.202000101
sha: a0c04832bc0c8459628ec2bc528cac8fbc09f38e
doc_id: 1006614
cord_uid: 38xcv2se

Widespread growth signatures in COVID‐19 confirmed case counts are reported, with sharp transitions between three distinct dynamical regimes (exponential, superlinear, and sublinear). Through analytical and numerical analysis, a novel framework is developed that exploits information in these signatures. An approach well known to physics is applied, where one looks for common dynamical features, independently from differences in other factors. These features and associated scaling laws are used as a powerful tool to pinpoint regions where analytical derivations are effective, get an insight into qualitative changes of the disease progression, and infer the key infection parameters. The developed framework for joint analytical and numerical analysis of empirically observed COVID‐19 growth patterns can lead to a fundamental understanding of infection progression under strong control measures, applicable to outbursts of both COVID‐19 and other infectious diseases.

COVID-19 pandemic introduced unprecedented worldwide social distancing measures. [1] While interventions such as quarantine or vaccination have been extensively studied in quantitative epidemiology, effects of social distancing are not well understood, [2] [3] [4] and when addressed, they have been studied only numerically. Unique opportunity to understand these effects has been provided by COVID-19 tracing through confirmed case counts, active cases and fatalities, in a variety of countries with different demographic and environmental conditions. [5, 6] We here show that focusing on analytical and numerical derivations in distinct epidemics growth regimes, is We implement the model deterministically, as COVID-19 count numbers are very high wherever reasonable testing capacities are employed. This makes model analytically tractable, and allows robust parameter inference through combination of analytically derived expressions and tightly constrained numerical analysis, as we show below. Our analysis is applied separately to each country, as the effect of social distancing, initial numbers of infected and exposed cases, diagnosis/detection efficiency and transmission rates may be different. However, within a given country, we do not take into account different heterogeneities−demographic, spatial, population activity, or seasonality effects. [2, 9, 10] Alternatively, global dynamical properties of the outbreak can be analyzed in a probabilistic framework employing partial differential equations in an age-structured model. [11, 12] These can readily be included in our model, but would lead to model structure which is not analytically tractable, so these extensions are left for future work.

Given this, the model equations are: (3) where N is the total population number; β-the transmission rate; σ -inverse of the latency period; γ -inverse of the infectious period; δ-inverse of the detection/diagnosis period; ε-detection efficiency; h-the recovery rate; m-the mortality rate. Social distancing is included through Equation (1) (second equation), which represents the rate at which the population moves (on average) from susceptible to protected category. The term 1 ( / ) 0 α + t t n corresponds to a sigmoidal dependence (similar to Fermi-Dirac function, in quantitative biology known as the Hill function [13] ). Time 0 t determines the halfsaturation, so that well before 0 t the social distancing is negligible, while well after 0 t the rate of transition to the protected category approaches α. Parameter n (the Hill constant) determines how rapidly the social distancing is introduced, that is, large n leads to rapid transition from OFF to ON state, and vice versa. [13] Equation (3) considers that only a fraction of the infected is diagnosed, so that εδI takes into account the diagnosis and the subsequent quarantine process.

To make the problem analytically tractable, we approximate the Hill function in the first relation of Equation (1) by unit step function, so that after 0 t the second term in Equation (1) becomes α − S and dominates over the first term, that is, ( ) ≈ α − S t e t . We checked that this approximation agrees well with full-fledged numerical simulations ( Figure 1D and Supporting Information). In all comparisons with analytical results, numerical analysis is done with the full model, allowing an independent check of both analytical derivations and employed approximations. Under this assumption, Equations (1) and (2) reduce to:

We next introduce two time regions: I) For ( ) I I t , we take ( 0) 0 = ≡ I t I , and restrict to dominant (positive) Jacobian eigenvalue, leading to the exponential regime:

By shifting

Equation (6) 

whose general solution for noninteger ν is given by:

where ( , ) ν J x represents Bessel function of the first kind, and , 1 2 C C are arbitrary constants. In our case 1 α γ σ α = + , 1

taking into account the following relation between standard and modified ( ( , ) ν I x ) Bessel functions of the first kind: [15, 16] ( , ) ( , ) ν ν = ν − I x i J ix , the general solution of Equation (6) reads:

To determine 1 C , 2 C , we use the following boundary conditions: (0) ( )

where the first derivative in region II has the following expression:

In obtaining the expression above, the following identities were frequently used: [15, 16] d ,

After derivations, where the following relation [16] 1, ,

together with sin(( 1) ) sin( ) ν π νπ ± =− and the identity relating modified Bessel function of the first and second kind [15, 16] we finally obtain a surprisingly simple result:

where ( , ) ν K x is the modified Bessel function of the second kind.

At maximum and inflection points,

, respectively. After extensive simplification of the results, this leads to (

is the basic reproduction number in the absence of social distancing: [6, 17] )

Equations (14) and (15) have to be solved numerically, but, as γ and σ are constants, we, interestingly, obtain that solutions will depend only on α . Since, for the analysis of superlinear and sublinear regimes, only the left inflection point and the maximum are important, we will further omit the second solution of Equation (15) (Equation (14) has one solution), and denote

(these two solutions are presented as upper and lower curves on Figure 2C , respectively), so that the effective reproduction numbers at inflection and maximum points ( e,i R and e,m R ) are: 

From this follows the length of superlinear regime (between inflection and maximum points):

We further Taylor expand ( ) II I t around the inflection point: 

In the superlinear regime ( ) ( )

where υ is the scaling exponent and s t marks the beginning of this regime. By Taylor expanding ( ) D t around i t , using Equations (18) and (3):

which is always larger than 1, as expected for the superlinear regime. As t i is localized toward the beginning of the regime,

. Finally, to provide analytical constrain on α , we Taylor expand ( ) I t II around the maximum: 

As

, we see that the quadratic term in Equation (20) is always negative, that is, ( ) D t curve enters sublinear regime around maximum of the infection. By fitting ( ) D t

f t t in this regime, and by using Equation (20) together with Equation (3), we obtain:

which allows to directly constrain α. We first numerically analyze outburst dynamics in the countries that continuously updated [18] three observable categories (D, A, and F). For a large majority of countries active cases were either not tracked or were not continuously updated, so the analysis is done for ten countries listed in the outline above.

In the exponential regime, the analytical closed-form solution is given by Equation (5) . From this, and the initial slope of ln( ) D curve (once the number of counts are out of the stochastic regime), β can be directly determined, while the corresponding eigenvector sets the ratio of 0 I to 0 E . The intercept of the initial exponential growth of D at 0 = t sets the product of 0 I and εδ . h and m can also be readily constrained, as from Equation (3), they depend only on integrals of the corresponding counts; here note that

Also, [17, 19, 20] 1/3 σ = day 1 − and 1/4 γ = day 1 − , characterize fundamental infectious process, which we assume not to change between different countries.

Only parameters related with the intervention measures ( , , , 0 α εδ t n ) are left to be inferred numerically, leading to tightly constrained numerical results. For this, we individually performed joint fit to all three observable quantities ( , , A D F) for each country. The errors are estimated through Monte-Carlo [21, 22] simulations, assuming that count numbers follow Poisson distribution.

Representative numerical results are shown in Figure 1 Information. In Figure 1A -C (and Supporting Information) we see a good agreement of our numerical analysis with all three classes of the case counts. In Figure 1D , we see sharp transitions between the three growth patterns indicated in the figure: i) exponential growth, observed as a straight line in log-linear plot in Figure 1E ; ii) superlinear growth, a straight line in loglog plot in Figure 1F ; iii) sublinear growth, a straight line in linear-log plot in Figure 1G .

Transition between the growth patterns can be qualitatively understood from Equation (3), and ( ) I t curve in Figure 1D . The exponential growth has to break after the inflection point of ( ) I t , that is, once the maximum of its first derivative ( ( ) ′ I t in Figure 1D ) is reached. In the superlinear regime, confirmed

), so this regime breaks once ( ) ′ I t (dashed blue curve) becomes negative. Equivalently, ( ) D t curve becomes concave (enters sublinear regime) once the maximum of the ( ) I t is reached. Note that the growth of ( ) D t can reemerge if the social distancing measures are alleviated. Our model can account for this by allowing transition from protected back to susceptible category, which is out of the scope of this study, but may improve the agreement with the data at later times (see Figure 1A -C). In addition to this numerical/intuitive understanding, we also showed that we analytically reproduce the emergence of these growth regimes (Equations (6), (14) , (15) ). Can we also analytically derive the parameters that characterize these regimes?

The exponential regime is straightforward to explain, as described above. The superlinear regime is in between the left inflection point and the maximum of ( ) I t , so that infective numbers grow, but with a decreasing rate. While the derivations are straightforward in the exponential regime, they are highly non-trivial during the subsequent subexponential (superlinear and sublinear) growth. As the superlinear regime spans the region between the left inflection point (17), with 1/α ≈ dependence, so that weak measures lead to protracted superlinear growth (see Figure 2A ). This tendency is also confirmed by independent numerical analysis in Figure 2A , where for each individual country we numerically infer α and extract the length of the superlinear regime. Therefore, the duration of the superlinear regime indicates the effectiveness of introduced social distancing.

The scaling exponent υ of the superlinear regime is given by Equation (19) , and shown in Figure 2B , where we predict that all countries are roughly in the same range of 1.2 1.5 υ < < (surprisingly, weakly dependent on α ), despite significant differences in the applied measures, demographic and environmental factors. This result is (independently from our model) confirmed from case count numbers (the slope in Figure 1F , and equivalently for other countries, see Figure 2B ).

How the effective reproduction number e R changes during this regime, that is, between the left inflection point and the maximum of ( ) I t ? e R quantifies the average number of secondary cases per infectious case, so that 1 e > R signifies disease outburst, while for 1 e < R the disease starts to be eliminated from the population. [17] The Equation (16) provides expressions for e,i R (at the inflection point) and e,m R (at the maximum). Interestingly, from Figure 2C , we observe that e,i R and e,m R do not depend on 0,free R and are, respectively, significantly larger and smaller than 1, which shows that transition from infection outburst to extinguishing happens during the superlinear growth. Consequently, the steepness of e R change over the superlinear regime significantly increases (larger change over smaller time interval, see Figure 2C ) with the measure strength.

Finally, in the sublinear regime, in a wide vicinity of ( ) I t maximum (which marks the beginning of the sublinear growth) leading non-linear term of ( ) D t is cubic ( 3 ≈ t , with negative prefactor). This is consistent with the expansion of ( ) I t around m t , which has leading negative quadratic ( 2 t ) dependence (see Equations (3) and (20)). The ratio between the prefactors in ( ) D t expansion is given by Equation (21), from which we see that α can be directly constrained, as shown in Figure 2D . For the ten countries with consistent tracking of , D A, and F, we independently numerically determined α and compared it with analytical results coming from Equation (21), obtaining an excellent agreement between our derivations and numerical results. The obtained α values should be understood as an effective epidemic containment measure-that is, estimating the true result of the introduced measures, which can be used to evaluate the practical effectiveness of the official policies.

To demonstrate how constraining α can aid numerical analysis in the cases when A is not continuously tracked, we next analyze five additional countries listed in the outline above, so that altogether our study covers majority of COVID-19 hotspots, which (at the time of this analysis) are close to saturation in confirmed counts. Furthermore, in the specific cases of UK and Italy, where we analytically obtained both very low and very constrained α (0.01 0.04 α < < ), we chose five times larger parameter span in α in the numerical analysis, to confirm that these low values are indeed preferred by the exhaustive numerical search. For example, the finally obtained α for Italy (0.033 0.005) ± and UK (0.025 0.005 ± ), together with previously obtained agreements shown in Figure 2A -C, strongly confirm that the observed growth patterns provide invaluable information for successful analysis of the infection progression data.

To further illustrate this, the synergy of analytical derivations and numerical analysis presented above enables us to, directly from the publicly available data, infer key infection parameters necessary to assess epidemics risks (provided in Table S1 , Supporting Information). We estimate these parameters by the same model/analysis, for a number of diverse countries, allowing their direct comparison. In Figure 3 , we show together case fatality rate (CFR), infected fatality rate (IFR) and infection attack rate (AR). [17, 24] CFR is the number of fatalities per confirmed cases. CFR can, in principle, be inferred directly from the data, but since different countries are in different phases of infection, we project forward the number of confirmed cases until a saturation is reached for each country, from which we calculate CFR. IFR (crucial parameter for assessing the risks for infection progression under different scenarios) is the number of fatalities per total number of infected cases, which is a genuine model estimate, due to the unknown total number of infected cases. AR (necessary for understanding the virus recurrence risk) is also determined from our model and provides an estimate of the fraction of the total population that got infected and possibly resistant.

From Figure 3 , we see that CFR takes very different values for different countries, from below 2% (New Zealand) to above 20% (France). On the other hand, IFR is consistent with a constant value (the dashed red line in the figure) of 0.3 0.4% ≈ − . In distinction to IFR, AR also takes diverse values for different countries, ranging from 1% ≈ to as high as 15% ≈ (though with large errorbars). Although diverse, these AR values are well bellow the classical herd immunity threshold of 60-70%

To summarize, we here developed a novel quantitative framework through which we showed that: i) The emergence of three distinct growth regimes in COVID-19 case counts can be reproduced both analytically and numerically. ii) Typically, a brief superlinear regime is characterized by a sharp transition from outburst to extinguishing the infection, where effective reproduction number changes from much larger to much smaller than one; more effective measures lead to shorter superlinear growth, and to a steeper change of the effective reproduction number. iii) Scaling exponent of the superlinear regime is surprisingly uniform for countries with diverse environmental and demographic factors and epidemics containment policies; this highly non-trivial empirical result is well reproduced by our model. iv) Scaling prefactors in the sublinear regime contain crucial information for analytically constraining infection progression parameters, so that they can be straightforwardly extracted through numerical analysis. Interestingly, we found that the number of COVID-19 fatalities per total number of infected is highly uniform across diverse analyzed countries, in distinction to other (highly variable) infection parameters, and about twice higher than commonly quoted for influenza (0.3-0.4% compared to 0.1-0.2%), which may be valuable for direct assessment of the epidemics risks.

While state-of-the-art approach in epidemiological modeling uses computationally highly demanding numerical simulations, the results above demonstrate a shift of paradigm toward simpler, but analytically tractable models, that can both explain common dynamical features of the system and be used for straightforward and highly constrained parameter inference. This shift is based on a novel framework that relates universal growth patterns with characteristic points of the infective curve, followed by analytical derivations in the vicinity of these points, in an approach akin to those in a number of physics problems. The framework presented here can be, in principle, further extended toward, for example, including stochastic effects or different heterogeneities such as age-structure. However, these are non-trivial tasks, and it remains to be seen to what extent the analytical results can be obtained in those more complex models. Overall, as our approach does not depend on any COVID-19 specifics, the developed framework can also be readily applied to potential outbursts of future infections.

Supporting Information is available from the Wiley Online Library or from the author. Figure 3 . CFR, IFR, and AR, inferred for countries whose abbreviations are indicated on the horizontal axis, are denoted, respectively, by blue, grey, and green dots, with errorbars indicated by corresponding bands. The dashed red horizontal line stands for IFR consistent with a mean value (indicated in the legend). Values for PRC are from ref. [23] .

WHO report

Mathematical Tools for Understanding Infectious Disease Dynamics

Modeling Infectious Diseases in Humans and Animals

Worldometer 2020. COVID-19 Coronavirus Pandemic

Introduction to Bessel Functions

Standard Mathematical Tables and Formulae

Handbook of Mathematical Functions

Numerical Recipes: The Art of Scientific Computing

This work was supported by the Ministry of Education, Science and Technological Development of the Republic of Serbia.

The authors declare no conflict of interest.

The data used in this study are openly available in Worldometer at https://www.worldometers.info/coronavirus/, reference number [7] . Parameters inferred through the analysis are available in the supplementary material of this article.