The dynamics of entropy in the COVID-19 outbreaks

With the unfolding of the COVID-19 pandemic, mathematical modelling of epidemics has been perceived and used as a central element in understanding, predicting, and governing the pandemic event. However, soon it became clear that long-term predictions were extremely challenging to address. In addition, it is still unclear which metric shall be used for a global description of the evolution of the outbreaks. Yet a robust modelling of pandemic dynamics and a consistent choice of the transmission metric is crucial for an in-depth understanding of the macroscopic phenomenology and better-informed mitigation strategies. In this study, we propose a Markovian stochastic framework designed for describing the evolution of entropy during the COVID-19 pandemic together with the instantaneous reproductive ratio. Then, we introduce and use entropy-based metrics of global transmission to measure the impact and the temporal evolution of a pandemic event. In the formulation of the model, the temporal evolution of the outbreak is modelled by an equation governing the probability distribution that describes a nonlinear Markov process of a statistically averaged individual, leading to a clear physical interpretation. The time-dependent parameters are formulated by adaptive basis functions, leading to a parsimonious representation. In addition, we provide a full Bayesian inversion scheme for calibration together with a coherent strategy to address data unreliability. The time evolution of the entropy rate, the absolute change in the system entropy, and the instantaneous reproductive ratio are natural and transparent outputs of this framework. The framework has the appealing property of being applicable to any compartmental epidemic model. As an illustration, we apply the proposed approach to a simple modification of the susceptible–exposed–infected–removed model. Applying the model to the Hubei region, South Korean, Italian, Spanish, German, and French COVID-19 datasets, we discover significant difference in the absolute change of entropy but highly regular trends for both the entropy evolution and the instantaneous reproductive ratio.

tion, it is still unclear which metric shall be used for a global description of the evolution of the outbreaks. Yet a robust modelling of pandemic dynamics and a consistent choice of the transmission metric is crucial for an in-depth understanding of the macroscopic phenomenology and better-informed mitigation strategies. In this study, we propose a Markovian stochastic framework designed for describing the evolution of entropy during the COVID-19 pandemic together with the instantaneous reproductive ratio. Then, we introduce and use entropy-based metrics of global transmission to measure the impact and the temporal evolution of a pandemic event. In the formulation of the model, the temporal evolution of the outbreak is modelled by an equation governing the probability distribution that describes a nonlinear Markov process of a statistically averaged individual, leading to a clear physical interpretation. The time-dependent parameters are formulated by adaptive basis functions, leading to a parsimonious representation. In addition, we provide a full Bayesian inversion scheme for calibration together with a coherent strategy to address data unreliability. The time evolution of the entropy rate, the absolute change in the system entropy, and the instantaneous reproductive ratio are natural and transparent outputs of this framework. The framework has the appealing property of being applicable to any compartmental epidemic model. As an illustration, we apply the proposed approach to a simple modification of the susceptible-exposed-infected-removed model. Applying the model to the Hubei region, South

Introduction
Coronaviruses are one of the most significant threats to human society [1][2][3][4][5][6]. Limited to short outbreaks in the recent past [7][8][9], their pandemic-level potential was well known [10,11], yet most countries proved unprepared to cope with the so-called coronavirus infectious disease of 2019 . Revealed in the Hubei province, China, the novel coronavirus has spread all over the world. China responded with massive containment measures starting at the end of January 2020, which limited further contamination on the mainland [8,12]. In Europe, most individual states have responded with similar containment measures. However, there has been a lack of common European action. Strict or soft containment measures have been applied with different timeframes and specialized to individual health and socio-cultural systems, showing very different pandemic evolutions. At the time of writing, the main episode of COVID-19 is (in general) under control in China, South Korea, and continental Europe [11], despite the possibility of multiple waves. On the contrary, North America and South America are still in the middle of the pandemic, and a clear picture of the evolution of events is not possible yet.
The amount of data available allows various modelling techniques to be tested more robustly than in previous epidemics. However, no model (from the physicsbased to the purely data-driven) has been or is able to predict the long-term evolution of the pandemic accurately. (Conversely, short-term predictions are possible with some degree of accuracy [13,14].) There are several reasons behind this long-term unpredictability; an incomplete list includes the partial understanding of the phenomenon, the (many, or even infinite) missing variables, the high sensitivity of the model to parameters, the incomplete/inaccurate data acquisition scheme, and the lack of uniform measurement methods.
However, a profound reason that makes any long-term prediction difficult is the presence of endogenous variables (a well-known problem in social sciences [15]). The endogenous variables may involve local policies, socio-cultural aspects, human behaviours, and communication strategies, and they are typically difficult to model and measure. Since an epidemic evolves as a result of the interplay between the "natural evolution" of the disease and society/human interventions, a robust and generalizable microscopic model with complete characterizations of endogenous variables is challenging to build. Given this, in this study, we use macroscopic phenomenology-based modelling to gain insight into the epidemic dynamics. 1 Therefore, here, the goal is not to provide long-term numerical predictions, although the proposed modelling technique can be used for extrapolation.
Among various modelling options, susceptibleinfected-removed (SIR) types of compartmental models have gained wide popularity due to their simplicity and straightforwardness in interpreting the macroscopic phenomenology. A significant amount of SIRtype model based studies have already been carried out to investigate the transmission properties of the COVID-19, and an incomplete list includes [8,12,[16][17][18][19][20]. The spectrum of complexity of these models is broad. They can range from a minimum number of compartments (which offers a better generalization) to a large number of compartments (which offers a better local description). They can be deterministic (i.e. counting the deterministic number of individuals for each compartment) or stochastic (i.e. defining a joint probability measure of the number of individuals for each compartment). They can have different data acquisition schemes (from a simple frequentist analysis of the single parameters to a complete Bayesian inversion scheme). Finally, they can simply macroscopically describe the pandemic evolution of a given location (top-down approach) or include a spatial topological description (including mobility) and/or a different degree of spreading among individuals by including adjacent matrices (bottom-up approach).
Given a dataset, these models can be calibrated and offer new insights into the evolution of the pandemic. For example, they can shed light on how the pandemic developed by measuring the reproductive ratio (constant or time varying) and finally estimating the effectiveness of containment measures. This metric also allows for a comparison between different regions but does not provide a quantitative measure of the impact of the spread. On the other hand, the evolution of the number of infected and deaths provides a means of direct impact; however, they lack objectivity as they are strongly influenced by the different populations of the regions, the measurement strategies, and the unreliability of the data. Furthermore, they are not global metrics as they do not provide an objective and robust way to unify them into a single (scalar) measure. Therefore, there is a research gap on how to provide a macroscopic model-metric pair to compare different regions' performance and get new insights into various outbreaks.
This study aims to fill this gap by proposing a macroscopic stochastic model equipped with a global transmission metric based on entropy. In this context, the entropy evolution of the process is a metric that describes the degree of disorder (i.e. of impact) of an epidemic. This metric allows for an objective comparison between regions and provides a global measure of both the evolution and the impact of COVID-19 outbreaks. In particular, we propose a compartmental stochastic model that has the following characteristics. (i) Stochastic: the model describes a statistically averaged individual by a nonlinear Markov process with compartmental epidemic states. (ii) Time-dependent: the model parameters are decomposed onto generic basis functions (of time). (iii) Parsimonious: instead of conventional orthogonal basis functions (e.g. orthogonal polynomials, Fourier/wavelet series) the adaptive basis functions are adopted to achieve a representation with minimum number of basis functions. (iv) Bayesian: the time-dependent parameters are assumed to be random and are calibrated by full Bayesian inversion. Furthermore, we equip the model with a metric based on entropy which has the following characteristics. (i) Meaningful: the metric provides a physical and transparent measure of the COVID-19 impact in a given region; moreover, it is by definition the time integral of the entropy rate, which represents the temporal evolution of the epidemic. (ii) Global: the metric provides a global and average description of the pandemic event. (iii) Consistent: the metric is not influenced by the number of individuals and can be used objectively to compare different regions. (iv) Robust the metric is associated with an error that is a direct output of the Bayesian inversion scheme used to calibrate the stochastic model. Finally, to have a reliable description of the events, we provide robust strategies to fill in missing information and to correct the numerous inconsistencies on the current datasets.
The paper is organized as follows. First, we develop general concepts of the proposed epidemic model, including governing equation, time-dependent parameterization, and Bayesian model calibration (Sect. 2). Second, we introduce the entropy-based metric in Sect. 3. Third, we apply the proposed approach to formulate a SEIR compartmental model for modelling the temporal evolution of COVID-19 (Sect. 4). Next, we apply the proposed approach to real datasets to the following regions: Hubei (China), South Korea, Italy, Spain, Germany, and France (Sect. 5). Finally, we conclude the study by identifying the limitations, conclusions, and future research directions.

The stochastic epidemic model
In the literature, the term "stochastic compartmental model" can refer to different formulations (see e.g. [21] for a review) with distinct underlying assumptions on the source of uncertainty. For instance, the noisedriven stochastic model is formulated by: (i) introducing additive noise process into the deterministic compartmental model; (ii) translating the noise into diffusion of probability distribution; and (iii) obtaining an equation of probability distribution (e.g. Fokker-Plank equation). Clearly, in the noise-driven model the source of uncertainty is the additive noise. An alternative and more popular stochastic formulation is the event-driven model, which can be summarized as a direct stochastic simulation of the deterministic model. Specifically, in the event-driven model the deterministic rate matrix is used to define the transition probabil- denotes the population in a compartment and ΔX the intra-state increment. With the transition probability, a direct stochastic simulation (via e.g. Gillespie's Direct Method [22,23]) would yield a random scenario of the epidemic. The proposed model can also be classified as an event-driven approach in the sense that the source of uncertainty is also the aleatory variability of transitions between epidemic states. However, instead of a stochastic simulation without a governing equation of probability distribution, the proposed model strictly follows an equation of probability distribution which describes a nonlinear Markov process. Consequently, the proposed model possesses clear physical interpretations within the mathematical framework of nonlinear Markov process theory.
Compartmental models with time-varying parameters have been widely studied in the literature [24][25][26][27][28]. A fundamental question to be addressed in timedependent models is the trade-off between over-fitting and under-fitting, or equivalently, model bias versus model variance. In the extreme scenario, a pointwise kernel-based parameterization may lead to an almost exact calibration on epidemic observations, yet the explanatory/extrapolation capability would be minimized, and the model variance would be maximized. In an over-parameterized model, the nonlocal trend/structure, which is crucial to characterize/understand the epidemic dynamics, can hardly be identified. In this study, we attempt to discover nonlocal structures from the epidemic dataset using a parsimonious formulation with adaptive basis functions. Moreover, since the model calibration is formulated in a Bayesian framework, likelihood-based model selection, e.g. using the Bayesian information criterion (BIC) [29], can be conveniently applied to specify the number of adaptive basis functions.

The original deterministic model
Consider a generic compartmental epidemic model with a fixed 2 total population N and a classification of the population into M compartments X = [X 1 , . . . , X M ] . The compartmental epidemic model describes the temporal evolution of the state vector X, where every component of X is by definition nonnegative and X is subjected to the conservation law For an infinitesimal incremental Δt, we study the following master equation of state vector.
where I is the identity matrix and H(X(t), t) is a problem-specific rate matrix (infinitesimal propagator). Equation (1) is equipped with the assumption that 2 A fixed total population indicates a closed system, i.e. (i) the vital dynamics (natural birth/death) is neglected, assuming that the course of the epidemic is relatively short; (ii) the immigration and emigration are neglected, assuming that N is sufficiently large and the course of immigration/emigration is relatively slow. the evolution of X(t) is smooth, 3 i.e. without jumps. Setting Δt → 0, Eq. (1) leads to Similar to mechanics, the rate matrix H(X(t), t) governs the dynamics of X(t). To preserve the conservation law X(t) 1 = N , we must have H(X(t), t) 1 = 0, where 1 is a vector of ones and 0 the null vector. Particularly, if X(t) eventually attains a stationary state X * defined as and define H * as We obtain the stationarity condition where 0 is a column vector of zeros.

Probabilistic reformulation
Given a deterministic H(X(t), t), Eq. (2) describes a deterministic trajectory of X(t). Since variabilities inevitably exist in the specification of H(X(t), t) or/and the initial condition, the solution X(t) becomes a multivariate stochastic process. However, the aforementioned "randomization" is regarded as epistemic with respect to the model Eq. (2). 4 This section focuses on a more fundamental (aleatory) probabilistic reformulation of Eq. (2). Adopting a frequentist point of view on probability, consider a normalization of X(t) by where X (t, N ) is used to highlight that the compartmental population depends on N , and P(t) can be interpreted as the marginal probability distribution of a discrete-state continuous-time stochastic process. Observe that despite P(t) being equivalent to proportions in a deterministic model, the probabilistic individualistic interpretation leads to a fully stochastic dynamic interpretation of the problem. The underlying state associated with P(t) is an epidemic state of a statistically averaged individual. Analogous to Eqs.
(1) and (2), we obtain and where Q( P(t), t) is a rate matrix analogous to H(X(t), t) in the deterministic model. Since Q( P(t), t) explicitly depends on P(t), Eq. (8) describes a nonlinear Markov process [30]. The conservation of probability is guaranteed by Q(t) 1 = 0. Equation (7) provides a straightforward strategy for sampling random realizations of the process. In particular, for a fixed initial condition P(t 0 ), the solution

P(t) is deterministic, and Q( P(t), t) can be regarded as Q(t) with P(t) being a time-dependent parameter of Q(t). The resulting tangent non-homogeneous
Markov process has the following transient stochastic matrix In line with the macroscopic description [Eq. (6)], the initial condition P(t 0 ) and the rate matrix Q( P(t), t) are by definition exactly the same for all N individuals. This assumption corresponds ad verbum to fix a constant average number of contacts and other interaction parameters between persons per unit time. Given this, there is an implicit assumption of statistical independence among the N individuals. In an adiabatic system, this is equivalent to letting N particles following N -independent Brownian motions. Therefore, this macro-description is emerging from the micro-behaviour of individuals interacting according N -independent Brownian motions, and the virus is spreading according to a simple diffusive process. 5 Consequently, a macro-random scenario of an epidemic can be obtained via simulating N -independent and identically distributed processes from Eq. (8). 5 At the micro-level, the spreading of the virus can be better described through a stochastic branching process. The simple diffusion process can be regarded as the result of a coarse graining.
In contrast to this macro-description, one could adopt a topological structure of the interactions between different individuals. This is generally done by including an adjacency operator which accounts for the different structure of the interactions among individuals (e.g. including mobility information, or considering the presence of superspreaders). As a consequence, each individual (or group of individuals) has a different average number of contacts and different interaction parameters. This leads to a heterogeneous compartmental model, which is inevitable dependent on a specific geographical area or social system. In this study, we focus on the general transmission trend of large regions, so that the trends can be more easily extrapolated and interpreted. Therefore, the simple macro-description is adopted. It is a specific choice which leads to a novel entropy-based measure to macroscopically compare the epidemic scenarios in different regions.

Time-dependent parameter model
We assume the "correct" model of Q( P(t), t) cannot be discovered, and Q( P(t), t) is replaced by a parametric model with a set of parameters α(t), i.e.

Let α(t) represent an arbitrary component of α(t).
A generic approach to parameterizing α(t) is to consider an expansion of the following form where w i are coordinates of basis functions ψ i (t). A popular choice for the basis function is the orthogonal polynomials, e.g. Legendre/Hermite/Laguerre/Chebyshev polynomials. An issue with orthogonal polynomial basis is that it may require high-order terms to represent a complex function, and consequently, this leads to over-fitting and implausible extrapolations. A powerful alternative is to use adaptive basis functions with the form where w i are parameters of the adaptive basis. The benefit of using Eq. (12) instead of Eq. (11) is that a parsimonious representation can be formulated, at the cost of introducing additional parameters in bases. An attractive choice for the adaptive basis ψ i (t, w i ) is the sigmoid function, i.e.
The theoretical justification of using Eqs. (13) in (12) is the universal approximation theorem [31], and the resulting parametric function is in fact a feed-forward neural network with a single hidden layer. In addition, the initial condition of Eq. (8) is unknown, and we parameterize where β m are nonnegative and subjected to the linear constraint M−1 m=1 β m ∈ [0, 1]. Note that β is time independent in the sense that the starting time point can be fixed. Therefore, the full parameter set of the epidemic model is written as θ := w, w , β .

Model calibration
The goal of model calibration is to find the optimal θ using real observation. We let D denote the dataset of observations collected for an epidemic up to some reference time point. The dataset D is composed by discrete measures on the number of persons in each observable compartment (e.g. infected, recovered, and dead), and D is a matrix of dimension M o × T , where M o denotes the number of observable compartments and T denotes the number of observed unit time (e.g. days).
The likelihood function L(D|θ) measures the probability of observing D given the model specified by θ . Using Bayes rule on θ, we have where π(θ |D) is the posterior distribution of θ conditional on the observed dataset D and π(θ ) is the prior distribution of θ . The major challenge of using Eq. (15) in practice is to sample from the posterior, and typically, this can be handled by advanced Markov chain Monte Carlo methods. The likelihood L(D|θ) may depend on both the observation error and the inherent variability of the epidemic model. Even by setting the observation error to zero, for any specified θ the prediction from the model is still random. If the accumulated numbers are of interest, e.g. the total number of recovered, for a large population size the variability in the prediction is expected to be small. Specifically, in a multinomial model the marginal coefficient of variation is proportional to 1/ √ N P m (t). However, at the same time, the model prediction can be extremely sensitive to θ , and an almost negligible perturbation due to the (albeit small) randomness of θ may lead to noticeably different predictions. Therefore, the Bayesian analysis is meaningful with or without the observation error.
To formulate the likelihood function, we first denote an individual (directed) random walk among various states as a Boolean operator 6 The vectors y (n) j (of dimension M o × 1) represent the state of the n person at time t j . Therefore, the components are all zero with exception of the current state, which takes the value of one. The joint probability density function of Y (n) , denoted by f (Y (n) |θ), is readily available from the governing equation Eq. (8). Next, we note that the observation D represents a collective scenario of the N -independent (under the assumptions of the macromodel) Markov processes y (n) j . Therefore, by brute force, the (observation-errorfree) likelihood function has the following form (16) where 1(·) is an indicator function and -fold summation. This brute force summation contains impossible paths that, however, are naturally excluded by the indicator function. Observe that this likelihood is fundamentally different from the deterministic compartmental models based on proportions rather than individual probabilities. Moreover, it is also different from the classical binomial (and related) likelihood approaches (used in direct Gillespie's methods). Equation (16) is clearly computationally intractable.
To formulate a computationally tractable likelihood function, we use the Markovian property and rewrite L(D|θ ) as where D(t j ) denotes the observation at time point t j . The first term L(D(t 0 )|θ) can be easily computed from a multinomial distribution with probability vector P(t 0 ). The specific expression of L(D(t j )|D(t j−1 ); θ ) varies with the adopted epidemic model, yet it is typically in the multinomial form. All the ingredients to compute L(D(t j )|D(t j−1 ); θ ) are included in the marginal distribution P(t), and the stochastic matrix S(t j , t j+1 |θ) is expressed as Observe that due to the discretization of a continuoustime Markov process into a discrete-time Markov process, the matrix S(t j , t j+1 |θ) is "less sparse" than Q(t|θ). For example, in a finite time interval, the impossible event 2 → 4 in matrix Q may have a finite probability of occurring in matrix S (through e.g. 2 → 3 → 4). In fact, Eq. (18) can be interpreted as the result of applying Eq. (9) infinite times within the integration interval. For a simple illustration of concept in constructing the likelihood function, we consider a two-state system where state 1 can either move to state 2 or stay still, while state 2 can only stay still. We assume D(t j−1 ) records [100, 50] in occupations of states 1 and 2 (for a total of 150 Markov chains), and D(t j ) records [90, 60]. Given the aforementioned transition structure, we know 10 out of 100 chains at t j moves from state 1 to state 2. Therefore, the likelihood L(D(t j )|D(t j−1 ); θ ) is simply the binomial 100 10 P 10 1→2 P 90 1→1 , where the transition probability P i→ j can be directly read from Eq. (18). One may not be able to observe the populations in all compartments; in this case, the total probability theorem can be used to integrate the unobservable states out (see Sect. 4 for an example).

Addressing data unreliability
The likelihood function introduced above only considers the inherent stochastic variability of the model. In reality, on top of the inherent stochastic variability, the underlying errors/uncertainties of a reported dataset involve multiple alternative sources. A rigorous way to treat such unreliability of reported data is to introduce a distribution assumption on the error , and the likelihood function can be written as where L(D|θ , ) is the likelihood with a specified error, π( ) is the probability distribution of the error, and Ω represents the feasible domain of the error. Note that in general represents a set of discretized stochastic processes. Apart from the technical challenge of integrating the high-dimensional Eq. (19), the major challenge of incorporating the error is the specification of π( ). Clearly, an assumption on π( ) would reshape the likelihood function towards the shape of π( ), and an inappropriate assumption would generate artificial and even misleading transmission properties. Therefore, we adopt an indirect path to incorporate the unreliability of reported data. Specifically, we apply a kernel function κ(·) to the original error-free likelihood function, i.e.
The kernel function is selected to "flatten" the likelihood function so that the unreliability in the reported data can be, to some extent, captured. In this study, we consider an exponential kernel, and Eq. (20) is rewritten aŝ where n > 1 is a scaling factor. Clearly, if n = 1, L(D|θ ) is identical to the original likelihood L(D|θ), and if n → ∞,L(D|θ ) approaches uniform. Instead of a specific error distribution, in practice it is more likely to have a crude idea on the possible magnitude of the errors in the reported dataset. For a further simplification, we focus on the errors in the infected cases, since the causal structure of infected and recovered/dead would let the errors in infected eventually flow into recovered/dead. Therefore, the question left is to relate "the magnitude of errors in infected cases" to the n in Eq. (21). It turns out, as a consequence of a sequence of qualitative reasoning, a reasonable choice of n is to let where Δ infected represents the maximum increment of infected and Δ represents the possible error in the maximum increment of infected. Note that Eq. (22) is proposed as a crude guidance for setting the magnitude of n . The reasoning of Eq. (22) is described as follows. (i) In Eq. (21), if L(D|θ) is Gaussian, the effect of applying 1/n is to introduce a scaling factor of n to the covariance of L(D|θ ). (ii) The likelihood L(D|θ) is a product of multinomial kernels (see Eq. (17)), which can be approximated by Gaussian with the maximum variance (of the infected compartment) in the size of Δ infected . (iii) Equation (22) is obtained as one assumes the scaled variance (scaled by factor n ) has a similar magnitude as Δ 2 . 7 For example, if one has a crude idea that the error of infected can be 30% of the reported infected, using Eq. (22) one could set n ∝ 0.09Δ infected .

Entropy as a global transmission metric
In the literature, only a few studies have investigated the application of entropy in epidemics modelling [32][33][34]. Moreover, in the previous works the motivation and formulation of entropy, as well as the adopted epidemics modelling framework, are entirely different from the current study. The key feature of the proposed stochastic model is that entropy-based transmission measures can be naturally developed. Specifically, for a discretized time grid {t j , j = 1, . . . , T } and stochastic matrix S(t j , t j+1 ) (Eq. (18)), we consider the Shannon entropy rate expressed as For . Recall that the marginal distribution P(t) and the stochastic matrix S(t j , t j+1 ) vary with the initial condition P(t 0 ). Therefore, Eq. (23) and H(t 0 ) should be averaged over the posterior distribution of the initial condition (obtained from the Bayesian analysis). In evaluation of Eq. (23), the convention 0 log 0 ≡ 0 is adopted.
In a homogeneous Markov process, the entropy rate is constant, and one has the important theoretical result 7 This assumption implies if Δ = √ Δ infected , n is 1. In other words, if the error is in the size of the standard deviation of multinomial distribution, one cannot tell it is error or inherent stochastic variability. H(t 1 |t 0 ). In the proposed epidemic model, the Markov process is nonlinear and non-homogeneous. Therefore, the evolution of the entropy rate H(t j |t j−1 ) within a specified duration should be considered, and they characterize the evolution of the degree of disorder.
Using the Markovian property of the epidemic model in conjunction with the additive property of entropy, 8 the entropy H(t 0 , t 1 , . . . , t T ) has the concise form The entropy H(t 0 , t 1 , . . . , t T ) is a scalar, and it provides a global measure on the total degree of disorder for an epidemic scenario. An important feature (shared by the reproductive ratio) of the entropy rate and the total entropy is that they are quantitatively comparable across different regions. This is because the entropy-based measures are associated with the statistically averaged individual, which is similar to measuring the mean-field approximation of the complex epidemic dynamics system. Qualitative speaking, a large H(t 0 , t 1 , . . . , t T ) may be contributed by: (i) a large pulse-like H(t j |t j−1 ), i.e. the entropy rate reaches high values but stays (in high values) for a short period; (ii) a moderate flat H(t j |t j−1 ), i.e. the entropy rate evolves with moderate values for a long period. In an epidemic scenario, a large pulse-like evolution of the entropy rate implies that the virus reaches a significant proportion of population but damped out (through the accumulation of recovered/dead) fast, and a flat evolution implies that the epidemic spreads in a moderate severe state for a long time. To quantitatively analyse whether the entropy rate evolution is pulse-like or flat, we introduce a concentration measure to H(t j |t j−1 ). Specifically, we again adopt the concept of Shannon entropy such that the concentration measure of H(t j |t j−1 ) is defined as the inverse of the Shannon entropy of the normalized H(t j |t j−1 ), i.e.
whereH(t j |t j−1 )) = H(t j |t j−1 ))/H(t 0 , t 1 , . . . , t T ) is the normalized H(t j |t j−1 ). Note that the total entropy Fig. 1 Diagram of the modified SEIR model. The transition from susceptible to exposed involves the term α 1 (t)P 2 (t), indicating the exposed is contagious H(t 0 , t 1 , . . . , t T ) appears in Eq. (25) as the normalizing constant of H(t j |t j−1 ) (when H(t j |t j−1 ) is normalized into a probability mass function). Also note that the Shannon entropy instead of the variance-based measures is adopted since a large variance does not necessarily reflect a large dispersion (e.g. a mixture model with highly concentrated component densities could produce a large variance).
In this paper, we propose the entropy rate H(t j |t j−1 ), the entropy H(t 0 , t 1 , . . . , t T ), and the concentration factor C(H) as complements to the conventional reproductive ratio. Appendix A illustrates various attractive features of the entropy-based measures. In practice, instead of computing entropy-based measures for the original distribution vector P and the stochastic matrix S, one may need to reshape P and S to obtain measures of different scales. For example, a typical epidemic model may involve the states of recovered and dead. Naturally, one would prefer the scenario of "a large recovery probability and a small death probability" over "a large death probability and a small recovery probability." However, Eq. (23) or Eq. (24) does not differentiate between recovery and death, and the aforementioned two scenarios can have exactly the same entropy (rate). The conventional reproductive ratio measure has the same issue. To let the entropybased measures incorporate the concept of "high recovery probability is preferable over high death probability", one could reshape the distribution vector P and the stochastic matrix S by merging the recovery state with the infected state. Consequently, the entropy (rate) of the reshaped system would diminish the contribution from the recovered state and highlight the contribution from the dead state (see Appendix A for an example).

Application to COVID-19
In the light of the general framework introduces in Sects. 2 and 3, this section introduces a simple modification of the SEIR model with the exposed being also contagious.

Modified SEIR
The modified SEIR has a five-dimensional probability state vector P(t) described as follows: -P 1 (t): the (instantaneous probability of being) susceptible. -P 2 (t): the exposed.

The rate matrix Q( P(t), t) is written as
where α(t) = [α 1 (t), . . . , α 5 (t)] are non-negative parameters to be calibrated. Note that here for generality we write every parameter as time dependent; however, in practice it is typically sufficient to set only a few of them as time dependent. Figure 1 illustrates the flow between compartments of the modified SEIR model.

Likelihood function
The likelihood function can be derived as a simple application of the concepts introduced in Sect. 2.4. First, we introduce a compound state, denoted by 1 ∨ 2, to represent the state of being in either susceptible or exposed. The most important property of the compound state 1 ∨ 2 is that it is an observable, i.e. if an individual is not at the state of infected nor at the state of recovered/dead, it is in the compound state. We let P 1∨2→m (t j , t j+1 ) represent the transition probability from the compound state at t j to other states m, m = 3 (infected), 4 (recovered), 5 (dead) at t j+1 . Using the total probability theorem, p 1∨2→m (t j , t j+1 ) can be expressed as where P 1 (t j ) and P 2 (t j ) are solution of Eq. (8), and P 1→m and P 2→m can be obtained from Eq. (18). Next, we arrange the dataset vector D(t j ) in the form D(t j ) = D 1∨2 (t j ), D 3 (t j ), D 4 (t j ), D 5 (t j ) to, respectively, represent the instantaneous number of compound state, instantaneous number of infected, accumulative number of recovered, and accumulative number of dead. Let ΔD(t j , t j+1 ) := |D(t j+1 ) − D(t j )| represent the absolute difference between two consecutive dataset vectors. Before introducing the likelihood function, we introduce an additional assumption that the ΔD 1∨2 (t j , t j+1 ) number of Markov chains all transit to the state 3 (the infected). This assumption can be always (made) correct since: (i) if t j is sufficiently close to t j+1 , naturally one cannot jump to the recovered/dead state from susceptible/exposed; (ii) if t j+1 − t j is large, one could re-mesh the timescale and perform interpolation on the dataset, so that t j can always be close to t j+1 by construction. Finally, the conditional likelihood function L(D(t j+1 )|D(t j ); θ ) can be written as where (29) and In Eqs. (29) and (30), the notations are simplified to drop t j , t j+1 and θ . To avoid possible ambiguity, the simplification rules are: t j+1 ), P m ≡ P m (t j |θ ), and P m→m ≡ P m→m (t j , t j+1 |θ). Substituting Eqs. (28), (29) and (30) into (17), one obtains the complete likelihood function.

Transmission measures
To obtain the entropy-based measures, we reshape the five-dimensional vector P into [P 1∨2 , P 3∨4 , P 5 ] , and the corresponding stochastic matrix S is also reshaped (via the total probability theorem) accordingly. The reason to consider the compound state 1 ∨ 2 is because as a whole the state 1 ∨ 2 is an observable, therefore the possible errors in identifying the exposed can be marginalized out. The reason to consider the compound state 3 ∨ 4 is discussed in Sect. 3, and as a result, the concept "high recovery probability is preferable over high death probability" is correctly incorporated. In addition to the entropy-based measures, using the next-generation matrix approach [35], the instantaneous reproductive ratio of the modified SEIR model can be defined as Note that the instantaneous reproductive ratio R 0 (t) can be understood as the basic reproductive ratio of a tangent model defined as a constant model with parameters α equal to the instantaneous parameters α(t) at the reference time point t .

Modelling and computational details
As aforementioned, the whole parameter set θ not only involves parameters of the rate matrix Q ( P(t), t), i.e. α(t) (represented by w, w of basis functions), but also parameters to represent the initial state, i.e. β.
In the modelling practice, except α 3 , which is related to the mean incubation period, we calibrate all the other parameters (including the initial conditions) with Bayesian analysis. The mean incubation period, which is 1/α 3 in the model, is reported in various previous studies [11,36], and typically, it is around 5 and in the range of [3,7] days. Therefore, we set 1/α 3 as an epistemic random variable within [3,7]. The time-dependent parameters are modelled with sigmoid basis functions. The number of function basis for each parameter is determined in an additive manner. Specifically, we start with constant α and iteratively increase the number of basis functions until the variation in likelihood function value (Eq. (17)) or BIC index becomes small. The Gibbs sampling with a uniform proposal distribution for each component of θ is adopted to sample from the posterior distribution. The step size of the Gibbs sampling is adaptively tuned using the acceptance rate of the Markov chain [37,38]. The seed samples for the Gibbs sampler are selected in the neighbourhood of the posterior mode. This is obtained by sequential Monte Carlo method [39,40] combined with deterministic trust region optimization [41,42].

Datasets
For the studied regions, the time series of the populations of infected, recovered, and dead during January to May 2020 are used in model calibration. The data are collected from WHO, European CDC and Chinese CDC [11,43,44]. The regions considered in this study include: Hubei province, South Korea, Italy, Germany, Spain, and France. 9 For each region, the population size N is fixed to the most recent value reported by Worldometer [45]. We choose these countries/regions because they have the same order of population size (this is irrelevant to the entropy-based measure which is N independent), they applied different containment strategies, and they represent different cultures. Moreover, at the time of writing of this article the peak of the epidemic waves is passed. A complete and thorough analysis of a large number of regions is out of the scope 9 We planned to include also UK; however, the data on the recovered patients are not available yet (unfortunately). of the current study. In fact, here we focus primarily on the model and metric definition and their use.

Data correction
Due to abrupt counting policy changes and various corrections, the COVID-19 datasets for Hubei, Spain, and France not only violate the smoothness assumption 10 of the proposed modelling framework, but also contradict the fundamental fact that the accumulative number can only be non-decreasing. Therefore, the datasets must be corrected. It is obvious that the cluster/jump of data has a missing information, which is the (correct) time of occurrence. To obtain a consistent dataset we fill this missing information by using the expected time of occurrence with respect to the distribution of the previous events. Since the dataset is recorded daily, marginally they form a multinomial distribution along the discrete time axis. It follows that the missing time information can be filled by using the daily expected number of events. Specifically, let t J represent the time point when the jump/drop happens (for a specified compartment), and let ΔD J represent the magnitude of the data jump/drop. We perform a postprocess of the dataset expressed as follows: where t i = t 0 , t 1 , . . . , t J , and D(t i ) represents the cumulative 11 number at t i . Note that ΔD J could be negative. For an illustration of the correction, Fig. 2 shows the raw and the corrected datasets for Hubei province.

The overall epidemic dynamics of various regions
After performing model calibrations on datasets of various regions, we present a comparison analysis based on various transmission measures. Figure 3 shows the evolution of the entropy rate of COVID-19 outbreaks for each of the regions considered. This graph represents the time evolution of the degree of disorder (in Fig. 2 Raw and corrected datasets of Hubei province. There are two policy changes regarding the dataset: (i) in February 12, 2020, the diagnosis criterion was temporarily relaxed, and as a result, there is an artificial jump in the number of infected; (ii) in April 17, 2020, the cumulative number of infected and dead is altered by a constant jump, and the cumulative number of recovered is altered by a constant drop. The jumps/drops are marked by a rectangular in the figure. The populations of infected, recovered, and dead are corrected using Eq. (32). Note that for infected the correction is made on cumulative numbers, and then, the instantaneous infected is obtained by subtracting the accumulative recovered and dead terms of infections and deaths) introduced by the virus in an average statistical individual of the region. This graph reflects features of the daily evolution of infection and recovered/deaths, but it is fundamentally different from the evolution of each compartment. In fact, it has the key property of being objective and comparable between regions. Interestingly, the evolution of the entropy rate has a similar form for each region, but a significant difference in the magnitude of the disorder. In particular, the cumulative integral of the entropy rate represents the change of entropy in the system and, therefore, the total impact in a region. In Fig. 4-top panel, we report this impact measure for each of the regions considered. Based on this metric, Spain was the most affected region despite the epidemic wave hit the country later than Italy. On the opposite side, South Korea is the country with the least change in entropy, highlighting an effective combination of policies and cultural habits that limited the impact of the epidemic. This is probably due to the experience gained during the recent 2015 Middle East Respiratory Syndrome coronavirus (MERS-CoV) outbreak [46]. Also, Hubei's reaction, with extreme containment measures, has overall limited the impact of the epidemic. Germany has the smallest total entropy among studied European countries.
Interestingly, the peak of entropy rate for Spain, Italy, and Germany occurred in about the same period but with a different left tail behaviour (i.e. in the growing phase). On the other hand, the behaviour of the right tail (i.e. the descent phase) is similar, showing a fatter and longer tail. A similar asymmetry can also be observed in Hubei and South Korea. A deviation from this "classic" behaviour is represented by Hubei, which does not show this long tail behaviour but has a rather compact and almost symmetric shape. A surprising result is shown in Fig. 4-bottom panel. Although the impact in each country is significantly different, the concentration factor is similar to support the fact that the evolution of COVID-19 is similar for all outbreaks. The Hubei region is slightly deviating from this trend, showing a higher concentration factor corroborating the lack of a right fat tail and, therefore, showing a higher prevalence as an impulse. Figure 5 shows a comparison of the instantaneous reproductive ratio and death rate, together with the date of lockdown in each region. One can infer that the lockdown reduced R 0 (t) effectively. However, surprisingly, the most effective decrease has been observed in South Korea where no national lockdown has been imple-mented, but only local containment measures, and massive early-stage testing.
It is important to note that the modelling results are associated with the optimized parsimonious model for each region. Specifically, in an optimized parsimonious model the number of time-dependent variables as well as the number of adaptive basis functions for each time-dependent variable is optimized, in the sense that increasing the number would not noticeably improve the calibration accuracy and decreasing the number would significantly degrade the accuracy. Finally, for an illustration on the degree of accuracy the model has achieved, the model calibration results of Hubei are shown in Fig. 6. The calibration for the other countries and their limitations are reported in B.

Robustness on the transmission trend
A natural concern regarding the discovered transmission trend is that if the trend is a genuine underlying structure of the epidemic, or it is merely some artificial/superficial structures from the specific timedependent model. It is challenging to (perfectly) resolve this concern because a compartmental model (or any mathematical model) is inevitably an approxi-  21) is fixed to 100, assuming the error in the increment of infected is of the order of a few hundreds. The figure suggests a highly accurate calibration on data using at most one adaptive basis function for each parameter mation on the real epidemic. Moreover, even an exact model exists, it is still challenging (if not impossible) to accurately identify the model due to the presence of endogenous variables. However, at least we could show that the proposed framework is self-consistent. In C, we simulate artificial epidemics from analytical SEIR laws and investigate whether the proposed modelling approach could identify correct transmission trends.

Incorporating the undetected cases
In Sect. 5, the reported/observed population in each compartment is used to calibrate the model, and the kernel function in Eq. (21) only flattens the likelihood function instead of altering its intrinsic shape. Consequently, the model describes an epidemic scenario consistent with but also confined by the reported cases. An important missing issue to address is to incorporate the undetected cases to fully uncover the magnitude of the epidemic. A practical modelling strategy is to introduce a probability distribution assumption on the (possibly time-dependent) ratio between reported and undetected cases and rewrite the likelihood function similar to Eq. (19). Clearly, the critical ingredient is the model assumption on the undetected. The ongoing studies on blood test for antibodies of SARS-CoV-2 [47] can be useful for this future research direction.

Application to more complex compartmental models
Depending on the modelling purposes, one could introduce additional compartments, e.g. the tested/suspected, the ICU case, the female and male, the old and young, etc., to study the interactions between different groups. It is also straightforward to include spatial distributed information by including adjacency and incidence matrices. However, one should be aware that the model variance and the possibility of converging to local insignificant likelihood modes in general would increase with model complexity. Therefore, it would be crucial to collect robust prior knowledge regarding the modelling parameters.

Conclusions
In this study, we have proposed a stochastic compartmental modelling framework of epidemics equipped with entropy-based metrics to measure both the impact and the evolution of a pandemic event. The model belongs to the nonlinear Markov processes class, which allows a robust formulation and a natural setting for developing entropy-based metrics. In addition, we have provided a complete Bayesian inversion scheme to calibrate the model parameters with related uncertainties. Subsequently, we specialized the proposed structure to a modified SEIR model and the COVID-19 pandemic.
In particular, we used the framework to investigate six regions: Hubei, South Korea, Italy, Spain, Germany, and France. We showed that the change in entropy in the selected areas (which is associated with the impact of an epidemic) is significantly different. However, it is surprising to note that the dynamic evolution of pandemic waves shows very regular trends and very similar concentration measures.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.   Figure 7 illustrates entropy-based measures for the two systems. Next, we fix α 2 = 0.15 and α 3 = 0.05, and let α 1 vary within [0. 2,2], so that the basic reproductive ratio varies within [1,10]. The initial condition is set to P(t 0 ) = [0.99, 0.01, 0, 0]. Figure 8 illustrates the entropy rate, total entropy, and concentration factor. Both Figs. 7 and 8 imply that the basic reproductive ratio and the entropy-based measures describe different aspects of the epidemic dynamics, and they do not have a one-on-one mapping. Crudely speaking, the entropy rate characterizes the instantaneous intensity and flatness of the epidemics, while the total entropy measures the overall impact. Note that in Figs. 7 and Fig. 9 Illustration of the entropy-based measures using the reshaped P and S. The figure suggests that when the infected and recovered states are merged into a compound state, the contribution (to the entropy rate and entropy) from the death cases is increased. The concentration factor is not sensitive to this technique. Note that without using a reshaped P and S, the entropybased measures for the two systems will be exactly the same 8, the entropy-based measures are computed using the original distribution vector P and stochastic matrix S. In Sects. 3 and 4.3, it is mentioned that P and S can be reshaped to highlight the contribution from the death cases. Figure 9 illustrates the effect of this technique. Specifically, two SIR models with the same initial conditions P(t 0 ) = [0.99, 0.01, 0, 0] are considered. The first model has α = [0.6, 0.19, 0.01], and the second model is obtained by swapping the recovery and death rates of the first model. The basic reproductive ratios for the two models are both 3.0.

B Results of model calibration
The model calibration results for Italy, South Korea, Spain, France, and Germany are shown in Figs. 10, 11, 12, 13, and 14. In the figures, the solid lines correspond to posterior mean estimations, and the shaded areas correspond to {10%, 20%, . . . , 99%} credible intervals around the posterior mean. Fig. 10 Modelling the overall epidemic dynamics of Italy with the modified SEIR model. The α 1 (t), α 4 (t) and α 5 (t) are modelled with a single sigmoid basis, and α 2 (t) is modelled as a constant variable. The n in Eq. (21) is fixed to 1000, assuming the error in the increment of infected is in the order of a few thousands Fig. 11 Modelling the overall epidemic dynamics of South Korea with the modified SEIR model. The α 1 (t) is modelled with a single sigmoid basis, and α 2 (t), α 4 (t), and α 5 (t) are modelled as constant variables. The n in Eq. (21) is fixed to 100, assuming the error in the increment of infected is in the order of a few hundreds Fig. 12 Modelling the overall epidemic dynamics of Spain with the modified SEIR model. The α 1 (t) and α 5 (t) are modelled with a single sigmoid basis, and α 2 and α 4 are modelled as constant variables. The n in Eq. (21) is fixed to 1000, assuming the error in the increment of infected is in the order of a few thousands Fig. 13 Modelling the overall epidemic dynamics of France with the modified SEIR model. The α 1 (t) and α 5 (t) are modelled with a single sigmoid basis, and α 2 and α 4 are modelled as constant variables. The n in Eq. (21) is fixed to 1000, assuming the error in the increment of infected is in the order of a few thousands

C Results on robustness/self-consistent test on the transmission trend
We simulate a random epidemic scenario from a modified stochastic SEIR model with parameters α = [α 1 (t), 0.10, 0.20, α 4 (t), α 5 (t)]. The duration of the simulation is set to 40 days. The time-dependent parameters are specified as α 1 (t) = 0.6 − 0.5t/40, α 4 (t) = 0.05 + 0.30t/40, and α 5 (t) = 0.15 − 0.10t/40, and consequently, the instantaneous repro- Fig. 15 The calibrated time-dependent modified SEIR model. The α 1 (t), α 2 (t), α 4 (t), and α 5 (t) are modelled with a single sigmoid basis, and α 3 is modelled as a constant variable. The n in Eq. (21) is fixed to 1, assuming no error in the dataset. The figure suggests an accurate calibration. Even the unobserved exposed population and its initial condition are accurately identified. However, this is because the dataset is generated from an analytical model. For real dataset where a mathematical model is only an approximation, it is preferable to use real clinical data to specify the mean incubation period and its epistemic uncertainty (as it is performed in the main text)  16, it is seen that the bias is noticeably decreased ductive ratio decreases from 3.50 to 0.75 as t grows from 0 to 40 (days). The population size is fixed to N = 10 3 , and the initial condition is set to P(t 0 ) = [1 − 101/N , 100/N , 1/N , 0, 0]. The specification of a relatively large initial condition for the unobservable exposed state is to "challenge" the proposed approach and investigate whether the exposed population can be accurately identified even without observing it. We calibrate a time-dependent modified SEIR model, and the results are illustrated in Figs. 15 and 16. Finally, to verify that the bias in Fig. 16 is mainly caused by the inherent stochastic variation of the model, we calibrate the model on the expectation of the analytical model, and the results of transmission properties are illustrated in Fig. 17.