key: cord-0692010-ghdyp2bc
authors: Frank, T. D.; Smucker, J.
title: Characterizing stages of COVID-19 epidemics: a nonlinear physics perspective based on amplitude equations
date: 2022-03-16
journal: Eur Phys J Spec Top
DOI: 10.1140/epjs/s11734-022-00530-9
sha: ec50f80a0443aabe28e2c59ba37724ad4bb8f2ee
doc_id: 692010
cord_uid: ghdyp2bc

The relevant dynamics underlying COVID-19 waves is described from an amplitude space perspective. To this end, the amplitude dynamics of infected populations is considered in different stages of epidemic waves. Eigenvectors and their corresponding amplitudes are derived analytically for low-dimensional models and by means of computational methods for high-dimensional models. It is shown that the amplitudes of all eigenvectors as functions of time can be tracked through the diverse stages of COVID-19 waves featuring jumps at the stage boundaries. In particular, it is shown that under certain circumstances the initial, outbreak stage and the final, subsiding stage of an epidemic wave are primarily determined by the unstable eigenvector of the initial stage and its corresponding remnant vector of the final stage. The corresponding amplitude captures most of the dynamics of the emerging and subsiding epidemics such that the problem at hand effectively becomes one dimensional leading to a dramatic reduction of the complexity of the problem at hand. Explicitly demonstrated for the first-wave COVID-19 epidemics of the year 2020 in the state of New York and Pakistan are given.

The pandemic of the Coronavirus disease 2019 (COVID-19) claimed 2,000,000 lives worldwide in the year 2020 [1] . Despite the availability of vaccines since the beginning of the year 2021 in some parts of the world, in the first half of the year 2021 (i.e., from January 17 to July 4, 2021) 2,000,000 more deaths due to COVID-19 were registered worldwide [2] . At the end of September 2021, the death toll since the beginning of the pandemic was at 4,700,000 people and a total of 231,000,000 individuals had been diagnosed with the disease [3] .

Epidemic outbreaks and waves such as COVID-19 waves observed regionally can be understood within the framework of nonlinear physics as instabilityinduced phenomena [4] [5] [6] . While in general for such instability-induced phenomena the underlying relevant, low-dimensional dynamics can be identified with the help of amplitude equations [7] [8] [9] [10] [11] , the identification of the relevant dominant sub-dynamics of COVID-19 waves from an amplitude equation perspective is a research field in its infancy.

In previous studies it has been pointed out that the relevant unstable eigenvector [9] or order parameter [7, 11] determines the initial stage of a COVID-19 wave [12] [13] [14] . In particular, it has been demonstrated that a e-mail: till.frank@uconn.edu (corresponding author) the dynamics of an epidemic converges towards the dominant sub-dynamics along the unstable eigenvector [15] and it has been recognized that the amplitude equation perspective based on eigenvalues and eigenvectors naturally leads to a three-stage scheme of COVID-19 waves [16] [17] [18] [19] . Such stage schemes in turn have been a key tool to understand COVID-19 waves and possible intervention impacts [20] [21] [22] [23] [24] [25] [26] . In particular, they can be regarded as approximation to more fine-grained modeling approaches that assume that model parameters vary continuously (e.g., on a daily basis) [27, 28] . However, what has not been discussed so far is an approach to determine the amplitude dynamics for the outbreak and subsiding stages of COVID-19 waves in general for epidemiological models that involve an arbitrary number of infected compartments.

In Sect. 2, a fundamental decomposition of epidemic models into infected and non-infected compartments that has been frequently used in the literature [4, 5] will be discussed and it will be shown how to derive amplitude equations of infected subsystems in arbitrary dimensions. The stage concept of epidemics will be introduced in Sect. 3.1 and stages of COVID-19 waves will be discussed from the amplitude equation perspective both on an analytical (Sect. 3.2) and computational level (Sect. 3.3). The analytical discussion will be restricted to SIR and SEIR models [29] . Both types of models have been extensively used in the literature to examine COVID-19 waves [20, 21, 24, 25, 27, 28, 30] . The computational approach will be exemplified by means of high-dimensional models.

Let us consider epidemics that can be described in terms of individuals in n compartments or subpopulations. The compartment sizes are described by the variables X 1 , . . . , X n ≥ 0 that form the state vector X = (X 1 , . . . , X n ). The state space D 0 is given by lR n + . The state is assumed to satisfy a compartmental model [29] of the form

For analysis purposes, it has been suggested to split the compartments into those of infected and noninfected individuals [4, 5] . Without loss of generality, let the first m compartments with m ≥ 1 and m < n denote various types of infected individuals and the remaining compartments denote non-infected individuals. Accordingly, X + = (X 1 , . . . , X m ) and X − = (X m+1 , . . . , X n ) denote vectors of the respective subspaces D + and D − (with D 0 = D + ∪ D − ) and X = (X + , X − ). The disease-free state X st is defined by X st = (X + st , X − st ) with X + st = (0, . . . , 0). As indicated, it is assumed that the disease-free state is a fixed point (or stationary state) of Eq. (1). Let us split the righthand side function N into N = (N + , N − ). Then Eq. (1) can equivalently be expressed in terms of two coupled dynamical systems like d dt X + = N + (X = (X + , X − )) ,

It is known that under certain circumstances the stability of the infected subsystem determines the stability of the entire system [4, 5] . Therefore, in what follows, an amplitude equation description for the infected subsystem will be developed. The benefit of this approach is that the subsystem D + involves a smaller set of variables as compared to the entire system. An amplitude equation approach that focuses only on D + makes the analysis feasible. That is, when focusing on the infected subsystem dynamics, it might be possible to derive analytical expressions that could not be derived for the entire system. Likewise, when conducting a computational approach, the primary focus can be directed to the dynamics in the subspace D + .

The m-dimensional subspace D + with state vector X can be mapped to a m-dimensional amplitude space described by the amplitude vector A = (A 1 , . . . , A m ). In this context, it is frequently assumed that the linearized dynamics of the infected subsystem at the disease-free fixed point is independent of the variables of the non-infected subsystem [4, 5] . Let u = X − X st denote a perturbation in D 0 and

describe the linearized dynamics of Eq. (1) at X st with the linearization matrix L composed of coefficients L ik . Then the aforementioned assumption that the linearized infected system at X st is independent of the non-infected system implies that L can be decomposed into three matrices L + , A, and B like

where L + is the upper, left-corner m × m matrix with coefficients L ik and i, k ≤ m. Importantly, as indicated in Eq. (4), L ik = 0 holds for all i ≤ m and k ≥ m + 1. Consequently, the linearization of dX + /dt = N + (X) at X st is given by

where L + is the linearization matrix in D + . In what follows, it is assumed that the matrix L + exhibits m linearly independent (right) eigenvectors v k associated with the eigenvalues λ k . If the vectors v k are taken as column vectors, then they constitute the matrix M defined by

whose inverse matrix M −1 exists. Let w k describe the rows of the inverse matrix like

Then, w i v k = δ ik holds, where δ ik is the Kronecker symbol. Accordingly, w i and v k form a biorthogonal basis. w i are referred to as biorthogonal vectors or left eigenvectors. The matrix M is also called diagonalization matrix because M −1 L + M = D holds, where D is the diagonal matrix with eigenvalues λ k . The amplitude space spanned by the amplitude variables A 1 , . . . , A m can be defined with the help of the mapping

which implies

Equation (8) describes the superposition of a state in D + in terms of amplitude variables. As such, Eq. (8) defines a mapping of points from the amplitude space to the state space D + . Vice versa, Eq. (9) describes a mapping of points from the state space D + to the amplitude space. To derive the evolution equations for the amplitudes A 1 , . . . , A m , the dynamical system dX + /dt = N + (X) is decomposed into the linear part (5) and a remainder term R such that

Multiplying Eq. (10) by M −1 and using M −1 X + = A, we obtain

Finally, X + occurring in Eq. (11) is expressed like X + = M A and the identity M −1 L + M = D is used, which leads to

In addition, from Eq. (2) the evolution equation for X − can be cast into the form

Equations (12) and (13) read in components

with k = 1, . . . , m and j = m + 1, . . . , n, respectively. Equations (15) for k = 1, . . . , m provide the mdimensional amplitude space description of the infected subsystem of the epidemic under consideration. The amplitude space description in general does not correspond to an autonomous system. Due to the nonlinear terms w k R the amplitude dynamics depends in general on the non-infected variables X − , as indicated. Let λ max denote the eigenvalue with the largest real part. Then, we can distinguish between two cases.

If lR{λ max } > 0 holds the disease-free fixed point X + st = (0, . . . , 0) in D + is an unstable fixed point. Let k(max) define the index of the eigenvalue λ k(max) = λ max . For illustration purposes, let us assume that there is only a single positive eigenvalue (for an example in this regard see the SEIR model in Sect. 2.5). In this case, the fixed point X + st ∈ D + corresponds to a saddle point with one unstable direction defined by the unstable eigenvector v k(max) . All other directions in D + are characterized by stable eigenvectors. Consequently, the dominant dynamics takes place along v k(max) . The corresponding amplitude A k(max) increases in magnitude during the initial stage of the epidemic [12, 15] . Using the terminology of synergetics [7, 11] , the unstable eigenvector may be considered as order parameter of the infectious disease outbreak and the corresponding amplitude as order parameter amplitude.

Amplitudes may decay in magnitude for two reasons. First, if lR{λ max } > 0 holds then due to the impact of the nonlinear terms the amplitudes A 1 , . . . , A m may converges to zero despite the fact that the maximal eigenvalue exhibits a positive real part and the epidemic subsides. Second, in contrast to this kind of subsiding of an epidemic driven by nonlinear terms, let us consider the case lR{λ max } < 0 and let us assume that the linear terms in Eq. (14) dominate over the nonlinear terms. Furthermore, just as in the special case of a single positive eigenvalue λ max discussed above that is qualitatively different from all other (negative eigenvalues) λ j with j = k(max), let us assume that there exist a gap in the eigenvalue spectrum between λ max and the remaining eigenvalues λ j with j = k(max). To describe this gap, we consider the time constants τ max = 1/λ max and τ j = 1/lR(λ j ). If |λ max | |lR(λ j )| then τ max τ j holds. Consequently, the dynamics of A k(max) along v k(max) is slow relative to the dynamics of the amplitudes in other directions. In other words, the amplitudes A j with j = k(max) quickly decay to zero, while A k(max) decreases slowly. If so, the subsiding of an epidemic is determined by dynamics of A k(max) and direction v k(max) .

The derivation of the amplitude equations (14) holds for arbitrary dimension m. To begin with, let us briefly consider the trivial case of m = 1. To this end, let us consider the susceptible-infectious-recovered (SIR) model [29] 

where β > 0 denotes the effective contact rate and γ is the recovery rate of infectious individuals related to the recovery period T like γ = 1/T . We have m = 1, n = 3, X = (I, S, R). The space D + is one dimensional with X + = I. The space D − exhibits the state vector X − = (S, R). Let us consider the fixed point X st = (0, N, 0) with S st = N and perturbations u = (u 1 , u 2 , u 3 ) with u 1 = I, u 2 = S − N , u 3 = R. Linearization of the evolution equation of I at the fixed point yields the eigenvalue λ = β − γ [31] . Formally, the subspace D + can be mapped to a one-dimensional amplitude space describe by the variable A 1 . In fact, both spaces are identical, which implies I = A 1 . Substituting I = A 1 and S = N + u 2 into the evolution equation of I (see Eq. (16)), we obtain the amplitude equation

with p denoting the linear function p(A 1 ) = βA 1 /N . From the evolution equations of S and R, we obtain

Equations (17) and (18) correspond to the explicit forms of Eqs. (14) and (15) , respectively, in the case of the SIR model.

Next, let us illustrate the case m = 2 that was worked out previously [12, 15] for the susceptible-exposedinfectious-recovered (SEIR) model defined by [29] 

where β > 0 and γ > 0 denote the effective contact rate and the recovery rate again. The parameter α is the rate of progression from being exposed but noninfectious to being infectious. We have m = 2, n = 4, X = (E, I, S, R), X + = (E, I), and X − = (S, R). In what follows, we are interested in studying outbreaks of novel infectious diseases such as COVID-19. In this case, it is generally assumed that the whole population of interest is initially susceptible to the virus. Accordingly, the fixed point of interest reads X st = (0, 0, N, 0) with S st = N . In this case, the linearization matrix L + is given by

and exhibits the eigenvalues

where the plus sign holds for λ 1 and the minus sign for λ 2 . From Eq. (21) it follows for any parameter set α, β, γ > 0 that λ 2 < 0 holds and λ 1 = λ max > λ 2 .

The perturbation u = X − X st = (E, I, u 3 , R) involves u 3 = S − N ≤ 0. The amplitude equations can be found in previous studies [12, 13, 15] and read d dt

with

where |M | denotes the determinant of the eigenvector matrix:

The dynamics of the noninfected system is given by

Furthermore, the mappings ( (26) and

The amplitude space description of the SEIR model will be used below to show how to describe stages of COVID-19 epidemics in amplitude space. and how to track characteristic changes across such stages.

Frequently the impact of intervention measures to stop the spread of COVID-19 has been studied with the help of epidemic models [16] [17] [18] [19] [20] [21] [22] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] . In this context, it is often assumed that a given period of interest can be partitioned into several intervals or stages [20] [21] [22] [23] [24] [25] [26] . During each stage, the parameters of the virus-human system under consideration are assumed to be approximately constant. In contrast, across stages these parameters vary. Let us assume a number of s stages j = 1, . . . , s. Each stage j is given by an interval [t j−1 , t j ] such that t 0 and t s denote the initial and end time points, respectively, of the total period of interest. Model parameters can then be considered as step functions over time. They change discontinuously at the stage boundaries t j . This implies for the amplitude space description of the system under consideration that the key characteristics such as eigenvalues and eigenvectors vary across stages. To track an epidemic through a sequence of stages in amplitude space, for each stage the matrix L + , the eigenvalues λ k , eigenvectors v k and other coefficients must be computed that occur in the amplitude equation description defined by Eqs. (14) and (15) . The benefit of the amplitude space description is that under appropriate conditions (see Sect. 2.2) it allows us to determine for each stage the dominant dynamics.These ideas will be illustrated in what follows for analytical and computational approaches.

Let us consider an epidemic model for which the amplitude space description in terms of Eqs. (14) and (15) is explicitly given. For illustration purposes, let us consider the SEIR model discussed in Sect. 2.2 for which the amplitude space description is given in form of Eqs. (23) and (25) . Given the model parameters α, β, γ in stage 1 and the initial state X(t 0 ) and the corresponding state X + (t 0 ) = (E(t 0 ), I(t 0 )) in the subspace of the infected compartments, the initial amplitude vector A(t 0 ) can be computed from A(t 0 ) = M −1 (Stage 1)X + (t 0 ), which is explicitly given by Eq. (27) . Subsequently, the state space equations (19) and the evolution equations (23) and (25) of the amplitude space description can be solved numerically for stage 1. At the end of stage 1, that is, at the time point t 1 , it is assumed that at least one of the model parameters α, β, γ changes. The new initial state of stage 2 is given by the final state X + (t 1 ) of stage 1. Using this final state as initial state of stage 2, the initial amplitudes A(t 1 ) of stage 2 can be computed from A(t 1 ) = M −1 (Stage 2)X + (t 1 ). In general, this implies that the amplitudes described by the vector A change in a discontinuous manner at the stage boundary t 1 . That is, A(t 1 ) computed from stage 1, in general, does not correspond to the initial stage A(t 1 ) computed for stage 2. The reason for this is that the basis of the amplitude space (given for the SEIR model in terms of the eigenvectors v 1 and v 2 ) in stage 2 in general differs from the basis of the amplitude space used in stage 1. In short, while the state X + (t) as function of time is a continuous function at t 1 , the amplitudes A 1 and A 2 of the SEIR model (and in general the amplitudes A 1 , . . . A m ) are functions of time that are discontinuous at the stage boundary t 1 .

Given the initial states X + (t 1 ) and A(t 1 ) for stage 2 and the model parameters for stage 2, the state dynamics and amplitude dynamics for stage 2 can be computed numerically again from Eq. (19) and Eqs. (23) and (25) , respectively. If s = 2, then the analysis is completed after the second stage. If s > 2, then the second stage is followed by another stage and the aforementioned procedure is repeated until all stages have been processed.

Let us illustrate these general considerations for a simulated epidemic that is computed from the SEIR model (19) and involves two stages S1 and S2 that both are characterized by positive eigenvalues. It is assumed that in stage 1, λ 1 is positive and relatively large. At the end of stage 1, the effective contact rate β is reduced (which mimics the impact of intervention measures) such that λ 1 decreases but is still positive in stage 2. The remaining parameters α and γ remain constant across S1 and S2. For the simulation we used α = 0.1/d, γ = 0.4/d, N = 10000, β(S1) = 2.0/d for the first 20 days and β(S2) = 1.0/d after day 20. From Eq. (21) it then follows that λ 1 (S1) = 0.22/d and λ 1 (S2) = 0.10/d. The entire simulation period was t s = 100 days and t 0 = 0. We used as initial conditions E(0) = 10 exposed individuals, Figures 1 and 2 show different aspects of the simulated epidemic. The solid lines in panels A and B of Fig. 1 show the trajectories E(t) and I(t), respectively, as computed from Eq. (19) . The dashed lines show how E(t) and I(t) evolve if β would not change at t 1 = 20 days. Comparing the solid lines with the dashed lines, the impact of the decrease of the effective contact rate is visible. Accordingly, the increase of the exposed and infectious individuals becomes less dramatic when β is changed from β(S1) to β(S2) < β(S1). Panel A of Fig. 2 shows the amplitudes A 1 and A 2 computed from Eqs. (23) and (25) . In particular, at t 1 = 20 days the amplitudes have been set to their respective initial values A 1 (t 1 ) and A 2 (t 1 ) of stage 2. On the relatively large scale of 2000 individuals for A 1 , a small discontinuity of A 1 can be seen at t = 20 days in the top subpanel of panel A. The expected discontinuity of A 2 is clearly visible in the bottom subpanel of panel A. Panel B shows the trajectory X + (t) = (E(t), I(t)) as phase curve in the E-I subspace. The phase curve describes a loop that starts and ends close to the disease-free state E = I = 0. Importantly, the directions of the unstable eigenvectors v 1 of stages 1 and 2 are depicted in panel B as well. Panel B illustrates that in stage 1 the trajectory X + (t) follows the eigenvector v 1 of stage 1. When β changes at t 1 the trajectory quickly adjusts to follow the new eigenvector v 1 of stage 2. Subsequently, the trajectory branches off from the direction specified by v 1 of stage 2, completes the loop, and converges towards the disease-free state. Panel C highlights the adjustment dynamics at the beginning of stage 2. The unstable eigenvector v 1 determines the dynamics in stage 1 and part of the dynamics in stage 2. The sta- ble eigenvector v 2 plays a negligible role during these periods. Likewise, the amplitude A 2 during the first 20-30 days is large as compared to |A 2 | and in this sense dominates the dynamics during that initial period. In summary, the unstable eigenvector or order parameter and its amplitude play the key role in the simulated epidemic outbreak, as discussed in Sect. 2.2.

The equivalence of the state space description (19) of the SEIR model and the amplitude space description of the SEIR model via Eqs. (23) and (25) can be illustrated in two ways. First, the amplitude dynamics A 1 (t) and A 2 (t) can be used to compute the time course of the state variables E(t) and I(t) in a given stage j like

Equation (28) has been written to hold for the general case of an epidemic model with n state variables and m infected compartments. We substituted the amplitude functions shown in panel A of Fig. 2 into Eq. (28) . In doing so, we obtained the graphs E and I, which are plotted as full circles in panel A and B of Fig. 1 . As expected, the curves E and I derived from Eq. (28) are identical with the functions E and I obtained by solving Eq. (19) directly. A second way to demonstrate the equivalence of the state space and amplitude space descriptions is to derive the amplitude vector A from the state vector X + without solving the amplitude equations (23) . That is, we may compute A in a given stage j from

Again, Eq. (29) holds in general, that is, for any n-dimensional epidemic model (1) (25). Next, we modeled a two-stage epidemic for which the effective contact rate drops in the second stage low enough to change the fundamental form of the epidemic. That is, we modeled a bifurcation induced by the impact of highly effective intervention measures: a change from an unstable disease-free fixed point E = I = 0 with λ max > 0 to a stable one with λ max < 0. For the simulation we used again α = 0.1/d, γ = 0.4/d, N = 10,000, and β(S1) = 2.0/d for the first 20 days. After day 20 we put β(S2) = 0.3/d. From Eq. (21) it then follows that λ 1 (S1) = 0.22/d and λ 1 (S2) = −0.02/d. We used the same initial conditions as in the previous simulation. Figure 3 presents the simulation results. Panel A shows the trajectories E(t) and I(t) as computed from Eq. (19) . At day 20 there is a change from the increasing trend of the number of exposed individuals to a decreasing one. The number of infectious individuals follows this trend change with a short delay. Panel B present the amplitudes A 1 and A 2 as functions over time as computed from Eqs. (23) and (25) . As expected, A 1 dominates the amplitude dynamics and increases dramatically in stage 1, while A 2 remains almost constant. At day 20, the basis of the amplitude space description changes and the amplitudes exhibit discontinuities. Panel C highlights the discontinuity for A 1 . In stage 2 the fixed point E = I = 0 exhibits two negative eigenvalues. The simulation reveals that the linearized model (see Eq. (5)), which implies dA j /dt = λ j A j , captures the essential dynamics in stage 2: both amplitudes rapidly approach zero in a more or less exponential manner indicating that the epidemic is vanishing. Note that E and I are continuous functions at the stage boundary t 1 , whereas A 1 and A 2 exhibit jumps at that boundary. Nevertheless, A 1 captures not only qualitatively but also quantitatively most of the dynamics of E. The reason for this is that for the selected parameters we have v 1,E = 0.99 in stage 1 and v 1, A 1 (t). Panel D shows the trajectory X + (t) in the E-I subspace. The directions of the eigenvector v 1 for stages 1 and 2 are shown as well. Panel D demonstrates that during stage 1, X + (t) follows the unstable eigenvector v 1 . When λ max becomes negative at t 1 = 20 days, the trajectory quickly converges to the direction defined by the new eigenvector v 1 of stage 2, which corresponds in S2 to a stable eigenvector. The fast approach towards the direction specified by v 1 indicates that there is gap in the eigenvalue spectrum as discussed in Sect. 2.2. In fact, for stage 2 we obtain λ 1 = λ max = −0.02/d and λ 2 = −0.48/d and the corresponding time constants τ j = 1/|λ j | given by τ 1 = 47.9 days and τ 2 = 2.1 days, respectively. That is, A 2 exhibits a time constants that is by a factor 24 smaller than the time constant of A 1 .

In this sense, A 2 describes a fast dynamics, while A 1 describes a slow dynamics. The slow dynamics determines the vanishing of the simulated epidemic as can be seen in panel D.

Finally, the equivalence of the state space and amplitude space descriptions was examined. Just as for the previous example, we computed states from Eq. (28) and amplitudes from Eq. (29) . That are plotted as filled circles in panels A and B of Fig. 3 . We found that the solutions thus obtained were identical to the solutions obtained directly from the respective evolution equations.

The computational approach is based on Eqs. (1), (5), (6) , and (29) . Accordingly, an epidemiologic model that is given in the form of Eq. (1) is solved numerically. The linearization matrix L + and the eigenvector matrix M (see Eqs. (5) and (6)) are computed for each stage of the epidemic of interest. Subsequently, the amplitudes as functions of time are computed from Eq. (29) . Let us illustrate the computational approach for first-wave COVID-19 epidemics that have been fully worked out as three-stage epidemics in the literature but for which the relevant eigenvectors and amplitudes have not been derived.

First, let us consider the first-wave COVID-19 epidemic of 2020 in the state of New York. In an initial study, a ten-variable generalized SEIR model was developed to examine the COVID-19 outbreak in the state of New York [39] . In a subsequent study, the model was used to examine the time course of the entire first-wave COVID-19 epidemic in the state of New York during the first half of 2020 [19] . Detailed descriptions of the model can be found in Refs. [19, 39] . In what follows the model will be reviewed only briefly. The compartmental model addresses both non-quarantined and quarantined individuals. On the side of the non-quarantined individuals, it describes susceptible (S u ), exposed (E u ), and symptomatic infectious (I u ) individuals. On the side of the quarantined individuals, the model describes susceptible (S q ), exposed (E q ), hospitalized symptomatic infectious (I h ), and hospitalized symptomatic intensive care unit (ICU) cases (I icu ). Moreover, asymptomatic infectious COVID-19 cases (I a ) constitute a compartment of their own. Finally, the model also addresses the compartment of recovered COVID-19 cases (R) and the compartment of COVID-19 associated deaths (D). The model equations read [19, 39] 

where k 0 denotes the force of infection [29] or rate constant

.

The model involves the parameters p, q, f 1 , f 2 , θ j , θ q , α, σ u , σ q , φ, r, ν, δ u , δ h , δ a , δ icu , η a , η h , k 1 , ..., k 7 , and β, which are semi-positive. The parameters denote the probability of infection per contact (p), proportion of being quarantined (q), proportion of exposed E u who transition to I u (f 1 ), proportion of exposed E u who transition to I h (f 2 ), efficacy of quarantine (θ j ), general efficacy of quarantine (θ q ), quarantine rate of exposed E u (α), incubation period of exposed E u (1/σ u ), incubation period of exposed E q (1/σ q ), isolation rate of I u (φ), proportion of exposed E q who transition to I h (r), rate of progression to ICU case (ν), death rate of nonisolated infectious I u (δ u ), death rate of isolated infectious I h (δ h ), death rate of asymptomatic cases I a (δ a ), death rate of ICU patients (δ icu ), reduced infectiousness for asymptomatic cases (η a ), reduced infectiousness for isolated cases I h (η h ), removal rate of S q (k 1 ), removal rate of E u (k 2 ), removal rate of E q (k 3 ), removal rate of I u (k 4 ), removal rate of I h (k 5 ), removal rate of I a (k 6 ), and the removal rate of I icu (k 7 ). Note that the model equation for the recovered is not presented above because it does not play a role for the subsequent discussion.

In Ref. [39] , the model was fit to death data observed during the COVID-19 epidemic in the state of New York. Moreover, the parameter β was used to account for the impacts of intervention measures such as physical distancing between people and wearing face masks.

To study the impact of such intervention measures the parameter β was varied while all remaining parameter values were fixed. Consistent with this approach and in line with previous three-stage epidemic models [16, 18] , in Ref. [19] a three-stage application of the model (30) was discussed in which β was varied throughout three stages. The analysis focused on the five infected compartments E u , E q , I u , I u , I h . While the individuals in intensive care units were also infected, they were neglected in the analysis because individuals I icu could not affect the stability of the disease-free fixed point (they were assumed to be perfectly isolated, see Eqs. (30) and (31)). Focusing on the five-dimensional subspace X + = (E u , E q , I u , I u , I h ) the three stages were defined as follows. Stage 1 was defined as the stage of the unstable disease-free fixed point with β such that there was at least one positive eigenvalue. Stage 2 was defined as the bifurcation point at which intervention measured lowered β such that all eigenvalues with positive real parts became zero. Stage 3 was defined as the subsiding stage in which intervention measures pushed the parameter β to an even lower level such that all eigenvalues exhibited negative real parts. The fixed model parameters and the parameters β for the three stages can be found in Ref. [19] . Here we only report the eigenvalues λ 1 , . . . , λ 5 of the five-dimensional subsystems as obtained from L + at the disease-free fixed point with X + = (0, 0, 0, 0, 0). They are reported in Table 1 . By definition, stage 1 exhibited an eigenvalue with a positive real part. Stage 2 exhibited as a maximal eigenvalue with zero real parts. Finally, stage 3 exhibited eigenvalues for which all real parts were negative.

Panel A of Fig. 4 shows the death data [42] during the observation period considered in Ref. [19] , which was March 1 to June 30. In addition, it shows the model fit D(t) to the data for the best-fit parameters. That is, Eq. (30) was solved numerically for the three stages to obtained D(t). The two vertical lines indicate the boundaries between stages 1 and 2 and stages 2 and 3, respectively.

In Ref. [19] , the relevant eigenvectors and amplitudes were not determined. Therefore, let us apply the computational approach outlined in the beginning of the section to examine the evolution of the epidemic. As can be seen in Table 1 , in stage 1 there was only one eigenvalue with positive real part. Consequently, there was one unstable eigenvector (order parameter) v 1 (S1). Likewise, in stage 3 there was one real-valued eigenvalue with a maximal real part. The corresponding sta-ble eigenvector v 1 (S3) is considered as the remnant of the unstable eigenvector v 1 (S1) of stage 1. According to our discussion in Sect. 2.2, the order parameter and the corresponding remnant order parameter should determine the dynamics of stages 1 and 3. Panel B of Fig. 4 shows the phase curves of the COVID-19 epidemic as obtained from the solutions E u , E q , I u , I h , I a of Eq. (30) in the two-dimensional subspaces E u -E q , E u -I u , E u -I h , and E u − I a . The eigenvectors v 1 (S1) and v 1 (S3) were computed from the matrix L + , see Eq. (5) for stages 1 and 3, respectively. L + itself for the model (30) and X + = (E u , E q , I u , I u , I h ) reads

The projections of v 1 (S1) and v 1 (S3) are plotted in panel B as well. As can be seen, the epidemic followed closely the unstable eigenvector v 1 (S1) and its remnant v 1 (S3) in stages 1 and 3, respectively.

The trajectories E u , E q , I u , I h , I a were substituted into Eq. (29) to obtain the amplitudes A 1 , . . . , A 5 as functions of time over the three stages. As such the eigenvectors v 1 , . . . , v 5 were defined as the eigenvectors related to the eigenvalues λ 1 to λ 5 listed in Table 1 . That is, the eigenvalues were sorted according to the magnitudes of the real parts as shown in Table 1 . Using this assignment of the indices 1-5 to the eigenvalues, the eigenvector indices 1-5 were assigned in the same way. As a result, the indices 1-5 of the amplitudes were assigned in this way as well: A 1 to A 5 for all three stages were defined as the amplitudes corresponding to the eigenvalues listed in Table 1 . Figure 5 illustrates the time course of the amplitudes A 1 , . . . , A 5 . Panel A shows the real-valued amplitudes A 1 , A 2 , A 3 related to the eigenvalues that exhibited the largest real parts and were real-valued in all three stages (see Table 1 , rank 1 to 3). Panels B and C show the real-and imaginary-parts of the amplitudes A 4 and A 5 related to the two remaining eigenvalues. These eigenvalues and their corresponding amplitudes were realvalued in stage 1 but assumed complex numbers in stages 2 and 3 (see Table 1 again). As can be seen in panel A, in stages 1 and 3, the amplitude A 1 played the dominant role among the three amplitudes A 1 , A 2 , and A 3 . Comparing panel A with panels B and C it also becomes clear that A 1 with respect to A 4 and A 5 played the dominant role in stages 1 and 3. This is consistent with the phase curves shown in panel B of Fig. 4 that suggest that the epidemic evolved primarily along the eigenvector v 1 in stages 1 and 3. During stage 2, A 1 remained almost constant (see panel A of Fig. 5 ) due to the fact that λ 1 = 0. During stage 2, the amplitudes A 2 and A 3 (panel A) and the imaginary parts of A 4 and Fig. 4 was primarily due to changes of the remaining four amplitudes A 2 to A 5 related to eigenvalues with non-zero real parts. Subsequent to the first COVID-19 wave, the state of New York was hit by a second wave during the winter period 2020/2021. Figure 6 shows the reported cumulative deaths (circles) during the period from November 1, 2020 to May 31, 2021 [42] . The death toll increased from about 33,000 to 52,000 people. The data were fitted again to the three-stage scheme of COVID-19 waves based Eq. (30) . To this end, the number of susceptibles on November 1 was taken as the number of susceptibles as obtained from the previous simulation that remained after the first COVID-19 wave. According to our firstwave simulation, the number of susceptibles decreases by about 260,000 individuals. As in Ref. [19] the stage boundaries were varied to obtain an optimal fit between data (circles) and model (solid line).

Let us compare the first and second-waves on the basis of their eigenvalues. Table 2 shows the eigenvalues obtained from the three-stage analysis of the second wave. Comparing Tables 1 and 2, it can be seen that the eigenvalues for j = 2, 3, 4, 5 do not vary much across the first and second waves. That is, the non-dominant first-wave stage 1 eigenvalues were comparable to the non-dominant second-wave stage 1 eigenvalues (except maybe for j = 5). Likewise, the second stage and the third stage eigenvalues were comparable across the two waves. However, the maximal eigenvalue was considerably smaller during the second wave both in stage 1 (by a factor 2) and stage 3 (by a factor 6). Consequently, both the increasing stage and the subsiding stage of the second wave was characterized by an overall slower dynamics as compared to the first wave. Figure 7 compares the dynamics of the amplitude A 1 during stage 1 for the two waves. The amplitude Tables 1 and 2 A 1 as plotted in panel A of Fig. 5 is shown. In addition, A 1 as computed for the second wave is shown. As excepted, the second wave dominant amplitude increased at a slower rate as compared to the first wave dominant amplitude. Figure 6 also shows approximations of A 1 as exponentially increasing functions (gray dotted lines) with exponential coefficients λ(max) as given in Tables 1 and 2 for stage 1. The exponential approximations fit the numerically obtained (exact) solutions very well. As a second example let us consider the first-wave COVID-19 epidemic in Pakistan during the period from March to September 2020. An epidemiologic model was developed in Ref. [40] to address the rising COVID-19 case numbers during the first three months of the epidemic in Pakistan. In a subsequent study, the model was further worked out to yield a three-stage model to explain both the stages of increasing and subsiding numbers of the first-wave epidemic in Pakistan from March to September 2020 and to address the stage in-between those two stages [17] . As such, the model for the COVID-19 epidemic in Pakistan involves eight compartments, which are the compartments of susceptible individuals (S), exposed individuals (E), symptomatic infectious individuals (I s ), asymptomatic infectious individuals (I a ), hospitalized infectious individuals (I h ), and infectious individuals who are treated in intensive care units (I icu ). In addition, the model considers infected quarantined individuals (Q) as well as recovered individuals (R). The model can again be regarded as a generalized SEIR model and reads [17, 40] 

with the rate constant k 0 given by

The model parameters correspond to the proportion of asymptomatic infections (ρ), incubation period (1/ω), hospitalization rate of symptomatic infectious individuals (η), hospitalization rate of quarantined individuals (δ), quarantine rate of exposed individuals (κ), rate at which hospitalized individuals require ICU care (φ), and removal rates k 1 , . . . , k 6 of individuals out of certain compartments. The rate constant k 0 involves as model parameters the effective contact rate (β) and the relative transmissibility parameters of the asymptomatic (ψ) and hospitalized (ν) individuals. Note that in Eq. (33) the evolution equation for the recovered individuals has been neglected since it will be of no concern in what follows. Furthermore, note that while in Ref. [40] demographic terms have been considered, in Ref. [17] and in Eq. (33) those terms have been neglected because on the relative short period of about half a year they do not affect the dynamics [43] . Just as in the previous example, it was assumed that the intervention measures such as physical distancing between individuals could lower the effective contact rate β [17, 40] . Accordingly, the three-stage model worked out in Ref. [17] defined the stages as in the previous example. Let us describe the stages with respect to the infected five-variable subsystem E, I s , I a , I h , Q. Again, the variable I icu can be neglected because it does not play a role for the stability of the disease-free fixed point (E, I s , I a , I h , Q) = (0, 0, 0, 0, 0). In stage 1 (the beginning of the epidemic in Pakistan), the disease-free fixed point of the population in Pakistan was unstable, which triggered an exponential increase in COVID-19 cases.

As a reaction to the COVID-19 outbreak, intervention measures were implemented that putatively decreased the effective contact rate β. Consequently, the diseasefree fixed point was stabilized. Stage 2 describes the bifurcation point at which the fixed point exhibited a zero eigenvalue. The disease-free fixed point was about to become stable. Subsequent to stage 2, that is, in stage 3, the intervention measures showed their full power and decreased the effective contact rate β to a sufficiently low value such that all eigenvalues of the disease-free fixed point exhibited negative real parts. The fixed point was stable and the epidemic subsided.

To fit model parameters, the variable P c of cumulative COVID-19 cases was introduced. From Eq. (33) , it follows that P c satisfies [17] 

With the help of P c , the model parameters and in particular β were fitted to confirmed COVID-19 cases of Pakistan [17, 40] . The best-fit parameters for the threestage model can be found in Ref. [17] . Panel A of Fig. 8 shows the cumulative COVID-19 cases [42] during the 7-months period (about 210 days) from March 1 to September 30. The data shows the typically sigmoid, three stage pattern [16, 18] of a first stage of accelerating increasing numbers characteristic for an unstable disease-free fixed point, a stage with a linear increase characteristic for the bifurcation point, and a third stage featuring a second de-accelerating bend characteristic for a subsiding epidemic under a stabilized disease-free fixed point. Panel A also shows the model solution P c (t) as obtained by solving Eqs. (33) and (35) for the best-fit parameters. Moreover, the linearization matrix L + reads

The eigenvalues of L + of the three stages are listed in Table 3 and demonstrate that in stage 1 (stage 3) the disease-free fixed point was unstable (stable). While in Ref. [17] the three stage model was discussed within the framework of the state space perspective, the amplitude space perspective was not worked out. In particular, the eigenvectors and amplitudes relevant for the rise and decay of the COVID-19 cases were not determined. Therefore, let us identify the relevant amplitude dynamics underlying the time course of the epidemic shown in panel A of Fig. 8 . Stage 1 exhibited only one eigenvalue with positive real part (see Table 3 ). Consequently, the epidemic outbreak in Pakistan was determined by a single unstable eigenvector 33) and (35) . Vertical dotted lines show stage boundaries. B Phase curves (solid lines) in two-dimensional subspaces and the directions (dotted lines) defined by the order parameter v1(S1) of stage 1 and its stage-3 remnant v1(S3). In all four subpanels, the lower and upper directions refer to v1(S1) and v1(S3), respectively The eigenvector v 1 for stages 1 and 3 is shown there as well. As can be seen, the epidemic followed closely v 1 in both stages 1 and 3.

The functions E(t), I s (t), I a (t), I h (t), Q(t) were substituted into Eq. (29) to obtain the amplitudes A 1 , . . . , A 5 as functions of time across the three stages. Figure 9 shows the amplitudes as functions of time. Panel A shows the three real-valued amplitudes A 1 , A 2 , A 3 . Panels B and C show the real-and imaginary-parts of the amplitudes A 4 and A 5 related to the two remaining eigenvalues that were complex-valued in all three stages. Comparing the amplitudes A 2 , A 3 , A 4 and A 5 with A 1 , it is clear that A 1 at any time point was at least by a factor 10 larger in magnitude than the other amplitudes. In particular, the dynamics in stages 1 and 3 was completely determined by A 1 . In contrast, while A 1 was larger than all other amplitudes in stage 2, it did not vary much over time during that intermediate bifurcation stage. Therefore, during the stage at which the epidemic in Pakistan was at its bifurcation point, the epidemic was determined primarily by the remaining amplitudes A 2 to A 5 (which again is consistent with the fact that A 1 was the amplitude in stage 2 that exhibited a zero eigenvalue).

Fundamental concepts of physics have been extensively used to better understand and prevent the spread of COVID-19. In this context, amplitude equations are a key tool to describe the dynamics of inanimate and animate systems at instability points. The benchmark cubic amplitude equation of the form dA/dt = λA − A 3 is known to describe the emergence of convection rolls in fluid and gas systems heated from below [8, [44] [45] [46] . The emergence of Turing patterns in chemical systems [9, [47] [48] [49] [50] and patterns on animal skins [10] has been extensively studied from the amplitude equation perspective. The emergence of brain activity patterns [51, 52] and, on the behavioral level, switches between gaits and other human reactions [11, [53] [54] [55] [56] [57] have been modeled with the help of amplitude equations.

In the current study, compartmental models as frequently used in epidemiology have been described in terms of this kind of amplitude equations. It has been demonstrated theoretically and in applications that the amplitude of the unstable eigenvector captures the initial time course of epidemics, in general, and COVID-19 outbreaks, in particular. In doing so, the present work supports a general interdisciplinary perspective according to which the spread of COVID-19 in populations should be seen in the broader context of physical, chemical, and biophysical/biochemical bifurcation phenomena that have been argued time and again to hold across various disciplines [7, 11, 58, 59] . Approaches alternative to the amplitude equation approach to describe COVID-19 waves such as fitting waves in terms of Gaussian functions [60] and piecewise-linear functions [25] have been suggested. However, while these approaches possess their own benefits, they do not relate the observed data to instabilities, that is, they do not relate data to the physical entities that lead to COVID-19 outbreaks and outbreaks of other infectious diseases [4] [5] [6] .

From a practical perspective, the unstable eigenvector or order parameter can be used to determine how the sub-populations involved in the epidemic of interest change in the amount relative to each other. For example, we may ask, when the number of exposed individuals E increases by a certain amount ΔE, what is the corresponding increase ΔI in infectious individuals? If we distinguish between asymptomatic I a and symptomatic infectious individuals I s , we may ask, how are the changes ΔI a and ΔI s related to ΔE? The unstable eigenvector in the initial stage (and its remnant in the subsiding stage) can be used to obtain answers to questions of this kind (see panel C of Fig. 2 , panel D of Fig. 3 , and panels B in Figs. 4 and 8) . Accordingly, the eigenvector provides approximate estimates for such relational changes. In summary, the amplitude space description of epidemics comes with a dramatic simplification of the problem at hand and provides predictive tools that might be exploited to forecast the course of emerging or subsiding epidemics.

By the time of writing this article, another winter wave seems to hit the state of New York. Figure 10 presents reported cumulative deaths (circles) during the period from August 1 to November 31, 2021 [42] . The function shows a first bend (or slightly nonlinear, exponentiallike increase) during August. During the subsequent 3 months from September to November the number of COVID-19 associated deaths increased in an almost linear fashion (see the regression line). Let us use Fig. 10 to present some concluding remarks of the current study in the context of future challenges. As such, the linear increase is consistent with a stage 2 dynamics characterized by a vanishing maximal eigenvalue. Alternatively, the increase may be interpreted in terms of an endemic fixed point, that is, a fixed point that describes a population that exhibits a finite number of infected individuals at baseline. In fact, the possibility of COVID-19 pandemic to become endemic is currently under debate [61] . A fixed point with a finite number of infected individuals typically implies that out of this pool individuals decease due to the disease at a certain rate. This mechanism yields a linear increase of cumulative deaths. Having said that, the effect of virus mutations plays a key role for the future of the COVID-19 pandemic, in gen- Therefore, in the context of Fig. 10 , it is difficult to predict how the function will continue quantitatively in the future. The benefit of the nonlinear physics perspective is that it allows to discuss various qualitatively different scenarios. Amplitude equations of endemic fixed points may be discussed in a similar way as the stage 3 amplitude equations discussed in Sect. 2. Bifurcations that are induced by the emergence of virus mutations, destabilize COVID-19 free fixed points, and trigger new COVID-19 waves may be discussed in analogy to intervention-induced bifurcation that stabilize COVID-19 free fixed point and lead to the subsiding of COVID-19 waves. In summary, the nonlinear physics approach outlined above provides a powerful tool to address challenges in the understanding and the modelling of phenomena that we may encounter in the near and far future of the COVID-19 pandemic.

World Health Organization, COVID-19 Weekly Epidemiological Update 23

World Health Organization, COVID-19 Weekly Epidemiological Update 47

World Health Organization, COVID-19 Weekly Epidemiological Update 59

Mathematical Epidemiology of Infectious Diseases

Synergetics (An Introduction

Introduction to Nonlinear Sciences

Mathematical Biology

Determinism and Self-organization of Human Perception and Performance

Proc. Natl. Acad. Sci. USA

ICCMS 20: Proceedings of the 12th international conference on computer modeling and simulation, (Association for Computer Machinery, 2020)

Timeline data from Johns Hopkins Center for Systems Science and Engineering

Dynamic Patterns-The Self-organization of Brain and Behavior