Characterizing stages of COVID-19 epidemics: a nonlinear physics perspective based on amplitude equations.

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.

Introduction

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 instability-induced phenomena [4-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-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-14]. In particular, it has been demonstrated that 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-19]. Such stage schemes in turn have been a key tool to understand COVID-19 waves and possible intervention impacts [20-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 model [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.

Infected subsystems: an amplitude space perspective

Model formulation and decomposition of states

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 that form the state vector . The state space is given by . The state is assumed to satisfy a compartmental model [29] of the formFor analysis purposes, it has been suggested to split the compartments into those of infected and non-infected individuals [4, 5]. Without loss of generality, let the first m compartments with and denote various types of infected individuals and the remaining compartments denote non-infected individuals. Accordingly, and denote vectors of the respective subspaces and (with ) and . The disease-free state is defined by with . As indicated, it is assumed that the disease-free state is a fixed point (or stationary state) of Eq. (1). Let us split the right-hand side function into . Then Eq. (1) can equivalently be expressed in terms of two coupled dynamical systems likeIt 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 involves a smaller set of variables as compared to the entire system. An amplitude equation approach that focuses only on 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 .

Amplitude space description of the infected subsystem

The m-dimensional subspace with state vector can be mapped to a m-dimensional amplitude space described by the amplitude vector . 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 denote a perturbation in anddescribe the linearized dynamics of Eq. (1) at with the linearization matrix L composed of coefficients . Then the aforementioned assumption that the linearized infected system at is independent of the non-infected system implies that L can be decomposed into three matrices , A, and B likewhere is the upper, left-corner matrix with coefficients and . Importantly, as indicated in Eq. (4), holds for all and . Consequently, the linearization of at is given bywhere is the linearization matrix in . In what follows, it is assumed that the matrix exhibits m linearly independent (right) eigenvectors associated with the eigenvalues . If the vectors are taken as column vectors, then they constitute the matrix M defined bywhose inverse matrix exists. Let describe the rows of the inverse matrix likeThen, holds, where is the Kronecker symbol. Accordingly, and form a biorthogonal basis. are referred to as biorthogonal vectors or left eigenvectors. The matrix M is also called diagonalization matrix because holds, where D is the diagonal matrix with eigenvalues . The amplitude space spanned by the amplitude variables can be defined with the help of the mappingwhich impliesEquation (8) describes the superposition of a state in in terms of amplitude variables. As such, Eq. (8) defines a mapping of points from the amplitude space to the state space . Vice versa, Eq. (9) describes a mapping of points from the state space to the amplitude space. To derive the evolution equations for the amplitudes , the dynamical system is decomposed into the linear part (5) and a remainder term such thatMultiplying Eq. (10) by and using , we obtainFinally, occurring in Eq. (11) is expressed like and the identity is used, which leads toIn addition, from Eq. (2) the evolution equation for can be cast into the formEquations (12) and (13) read in componentswith and , respectively. Equations (15) for provide the m-dimensional 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 the amplitude dynamics depends in general on the non-infected variables , as indicated. Let denote the eigenvalue with the largest real part. Then, we can distinguish between two cases.

Epidemic outbreak

If holds the disease-free fixed point in is an unstable fixed point. Let define the index of the eigenvalue . 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 corresponds to a saddle point with one unstable direction defined by the unstable eigenvector . All other directions in are characterized by stable eigenvectors. Consequently, the dominant dynamics takes place along . The corresponding amplitude 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.

Subsiding epidemic driven by the linear system

Amplitudes may decay in magnitude for two reasons. First, if holds then due to the impact of the nonlinear terms the amplitudes 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 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 discussed above that is qualitatively different from all other (negative eigenvalues) with , let us assume that there exist a gap in the eigenvalue spectrum between and the remaining eigenvalues with . To describe this gap, we consider the time constants and . If then holds. Consequently, the dynamics of along is slow relative to the dynamics of the amplitudes in other directions. In other words, the amplitudes with quickly decay to zero, while decreases slowly. If so, the subsiding of an epidemic is determined by dynamics of and direction .

SIR and SEIR models

The derivation of the amplitude equations (14) holds for arbitrary dimension m. To begin with, let us briefly consider the trivial case of . To this end, let us consider the susceptible-infectious-recovered (SIR) model [29]where denotes the effective contact rate and is the recovery rate of infectious individuals related to the recovery period T like . We have , , . The space is one dimensional with . The space exhibits the state vector . Let us consider the fixed point with and perturbations with , , . Linearization of the evolution equation of I at the fixed point yields the eigenvalue [31]. Formally, the subspace can be mapped to a one-dimensional amplitude space describe by the variable . In fact, both spaces are identical, which implies . Substituting and into the evolution equation of I (see Eq. (16)), we obtain the amplitude equationwith p denoting the linear function . From the evolution equations of S and R, we obtainEquations (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 that was worked out previously [12, 15] for the susceptible–exposed–infectious–recovered (SEIR) model defined by [29]where and denote the effective contact rate and the recovery rate again. The parameter is the rate of progression from being exposed but non-infectious to being infectious. We have , , , , and . In what follows, we are interesting 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 with . In this case, the linearization matrix is given byand exhibits the eigenvalueswhere the plus sign holds for and the minus sign for . From Eq. (21) it follows for any parameter set that holds and . For () we have (). The eigenvectors of readThe perturbation involves . The amplitude equations can be found in previous studies [12, 13, 15] and readwithwhere |M| denotes the determinant of the eigenvector matrix: and . The dynamics of the non-infected system is given byFurthermore, the mappings and read explicitlyandThe 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.

Stages of COVID-19 epidemics: characterization in amplitude space

Epidemic stages

Frequently the impact of intervention measures to stop the spread of COVID-19 has been studied with the help of epidemic models [16–22, 32–41]. In this context, it is often assumed that a given period of interest can be partitioned into several intervals or stages [20-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 . Each stage j is given by an interval such that and 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 . 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 , the eigenvalues , eigenvectors 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.

Analytical approach

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 and the corresponding state in the subspace of the infected compartments, the initial amplitude vector can be computed from , 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 , 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 of stage 1. Using this final state as initial state of stage 2, the initial amplitudes of stage 2 can be computed from . In general, this implies that the amplitudes described by the vector change in a discontinuous manner at the stage boundary . That is, computed from stage 1, in general, does not correspond to the initial stage 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 and ) in stage 2 in general differs from the basis of the amplitude space used in stage 1. In short, while the state as function of time is a continuous function at , the amplitudes and of the SEIR model (and in general the amplitudes ) are functions of time that are discontinuous at the stage boundary .

Given the initial states and 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 , then the analysis is completed after the second stage. If , 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, 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 decreases but is still positive in stage 2. The remaining parameters and remain constant across S1 and S2. For the simulation we used /d, /d, , /d for the first 20 days and /d after day 20. From Eq. (21) it then follows that /d and /d. The entire simulation period was days and . We used as initial conditions exposed individuals, infectious individuals, and susceptibles.

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 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 to . Panel A of Fig. 2 shows the amplitudes and computed from Eqs. (23) and (25). In particular, at days the amplitudes have been set to their respective initial values and of stage 2. On the relatively large scale of 2000 individuals for , a small discontinuity of can be seen at days in the top subpanel of panel A. The expected discontinuity of is clearly visible in the bottom subpanel of panel A. Panel B shows the trajectory as phase curve in the E–I subspace. The phase curve describes a loop that starts and ends close to the disease-free state . Importantly, the directions of the unstable eigenvectors of stages 1 and 2 are depicted in panel B as well. Panel B illustrates that in stage 1 the trajectory follows the eigenvector of stage 1. When changes at the trajectory quickly adjusts to follow the new eigenvector of stage 2. Subsequently, the trajectory branches off from the direction specified by 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 determines the dynamics in stage 1 and part of the dynamics in stage 2. The stable eigenvector plays a negligible role during these periods. Likewise, the amplitude during the first 20–30 days is large as compared to 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.

Fig. 1

State space description of the first example of a simulated two-stage epidemic. The solid lines show plots of E (A) and I (B) versus time for the SEIR model (19) with a sudden change of at days. The dashed lines indicate the dynamics if would remain unchanged. Model parameters and initial conditions see text

Fig. 2

Amplitude and state space description of the first example of a simulated two-stage epidemic with a sudden change in at days. A The amplitudes and as functions of time. B The phase curve of the epidemic in the E-I subspace. The directions of are shown as a dotted lines for stages 1 and 2, respectively. C A detail of B

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 and can be used to compute the time course of the state variables E(t) and I(t) in a given stage j likeEquation (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 from the state vector without solving the amplitude equations (23). That is, we may compute in a given stage j fromAgain, Eq. (29) holds in general, that is, for any n-dimensional epidemic model (1) that involves m infected compartments and satisfies the requirements discussed in Sects. 2.1 and 2.2. We applied Eq. (29) to the state functions E(t) and I(t) shown in panels A and B of Fig. 1. The amplitudes thus obtained are shown in panel A of Fig. 2 as full circles. As expected, the amplitudes and based on Eq. (29) are identical with the amplitudes and computed directly from Eqs. (23) and (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 with to a stable one with . For the simulation we used again /d, /d, , and /d for the first 20 days. After day 20 we put /d. From Eq. (21) it then follows that /d and /d. We used the same initial conditions as in the previous simulation. Figure 3 presents the simulation results. Panels 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 and as functions over time as computed from Eqs. (23) and (25). As expected, dominates the amplitude dynamics and increases dramatically in stage 1, while 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 . In stage 2 the fixed point exhibits two negative eigenvalues. The simulation reveals that the linearized model (see Eq. (5)), which implies , 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 , whereas and exhibit jumps at that boundary. Nevertheless, 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 in stage 1 and in stage 2. Since holds (see Eq. (26)) and , we obtain . Panel D shows the trajectory in the E-I subspace. The directions of the eigenvector for stages 1 and 2 are shown as well. Panel D demonstrates that during stage 1, follows the unstable eigenvector . When becomes negative at days, the trajectory quickly converges to the direction defined by the new eigenvector of stage 2, which corresponds in S2 to a stable eigenvector. The fast approach towards the direction specified by indicates that there is gap in the eigenvalue spectrum as discussed in Sect. 2.2. In fact, for stage 2 we obtain /d and /d and the corresponding time constants given by days and days, respectively. That is, exhibits a time constants that is by a factor 24 smaller than the time constant of . In this sense, describes a fast dynamics, while describes a slow dynamics. The slow dynamics determines the vanishing of the simulated epidemic as can be seen in panel D.

Fig. 3

Second example of a simulated two-stage epidemic for the case in which the disease-free fixed point becomes stable at days due to a sudden change of . A and B Show state and amplitude space trajectories as computed from Eq. (19), on the one hand, and Eqs. (23) and (25), on the other hand. C Shows a detail of the function to highlight the discontinuity of at days. D Shows the phase curve E(t)-I(t) (solid line) and the directions of for stages 1 and 2 (dotted lines)

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.

Computational approach

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 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 (), exposed (), and symptomatic infectious () individuals. On the side of the quarantined individuals, the model describes susceptible (), exposed (), hospitalized symptomatic infectious (), and hospitalized symptomatic intensive care unit (ICU) cases (). Moreover, asymptomatic infectious COVID-19 cases () 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 denotes the force of infection [29] or rate constantThe model involves the parameters p, q, , , , , , , , , r, , , , , , , , , and , which are semi-positive. The parameters denote the probability of infection per contact (p), proportion of being quarantined (q), proportion of exposed who transition to (), proportion of exposed who transition to (), efficacy of quarantine (), general efficacy of quarantine (), quarantine rate of exposed (), incubation period of exposed (), incubation period of exposed (), isolation rate of (), proportion of exposed who transition to (r), rate of progression to ICU case (), death rate of non-isolated infectious (), death rate of isolated infectious (), death rate of asymptomatic cases (), death rate of ICU patients (), reduced infectiousness for asymptomatic cases (), reduced infectiousness for isolated cases (), removal rate of (), removal rate of (), removal rate of (), removal rate of (), removal rate of (), removal rate of (), and the removal rate of (). 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 . While the individuals in intensive care units were also infected, they were neglected in the analysis because individuals 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 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 of the five-dimensional subsystems as obtained from at the disease-free fixed point with . 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.

Table 1

Eigenvalues of the three stages of the first COVID-19 wave in the state of New York, 2020, as obtained for the infected five-dimensional subspace

	Rank
Stage	1	2	3	4	5
1	0.22/d	-0.28/d	-0.32/d	-0.42/d	-0.64/d
2	0/d	-0.28/d	-0.33/d	-0.42+i0.04/d	-0.42-i0.04/d
3	-0.06/d	-0.27/d	-0.33/d	-0.39+i0.04/d	-0.39-i0.04/d

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.

Fig. 4

Model-based analysis of the death data of the first-wave epidemic in the state of New York, 2020. A Cumulative deaths during the period from March 1 to June 30 (circles) and model fit (solid line) obtained from Eq. (30). Vertical dotted lines show stage boundaries. B Phase curves (solid lines) in two-dimensional subspaces and the directions specified by the dominant eigenvectors and in stages 1 and 3 (dotted lines). In all subpanels, the lower (upper) direction refers to ()

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) . Likewise, in stage 3 there was one real-valued eigenvalue with a maximal real part. The corresponding stable eigenvector is considered as the remnant of the unstable eigenvector 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 of Eq. (30) in the two-dimensional subspaces –, –, –, and . The eigenvectors and were computed from the matrix , see Eq. (5) for stages 1 and 3, respectively. itself for the model (30) and readsThe projections of and are plotted in panel B as well. As can be seen, the epidemic followed closely the unstable eigenvector and its remnant in stages 1 and 3, respectively.

The trajectories were substituted into Eq. (29) to obtain the amplitudes as functions of time over the three stages. As such the eigenvectors were defined as the eigenvectors related to the eigenvalues to 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: to 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 . Panel A shows the real-valued amplitudes 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 and related to the two remaining eigenvalues. These eigenvalues and their corresponding amplitudes were real-valued 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 played the dominant role among the three amplitudes , , and . Comparing panel A with panels B and C it also becomes clear that with respect to and 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 in stages 1 and 3. During stage 2, remained almost constant (see panel A of Fig. 5) due to the fact that . During stage 2, the amplitudes and (panel A) and the imaginary parts of and (panel C) decayed in magnitude. Interestingly, the real parts of and (panel B) increased in magnitude during the very short bifurcation period of stage 2. Overall, the dynamics in stage 2 was determined by variations in , , , and . In other words, the stage 2 dynamics shown in the phase curves in panel B of Fig. 4 was primarily due to changes of the remaining four amplitudes to related to eigenvalues with non-zero real parts.

Fig. 5

Amplitude description of the first-wave epidemic in the state of New York, 2020. A Amplitudes (thick solid line), (thin solid line), and (dotted line) as functions of time in all three stages S1, S2, and S3. B and C Real (B) and imaginary (C) parts of (solid line) and (full circles) in stages S1, S2, and S3. The real parts shown in B are identical in stages 2 and 3. The imaginary parts shown in C are mirror images of each other in all stages

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 first-wave simulation, the number of susceptibles decreases by about 260,000 individuals. As in Ref. [19] the stage boundaries were varied to obtained obtain an optimal fit between data (circles) and model (solid line).

Fig. 6

Model-based analysis of the death data of the second-wave epidemic in the state of New York, during the period from November 1, 2020, to May 31, 2021. Cumulative deaths (circles) and model fit (solid line) obtained from Eq. (30) are shown. Vertical dotted lines show stage boundaries

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 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 ). 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.

Table 2

Eigenvalues of the three stages of the second COVID-19 wave in the state of New York during the winter period 2020/2021

	Rank
Stage	1	2	3	4	5
1	0.10/d	-0.28/d	-0.32/d	-0.43/d	-0.51/d
2	0/d	-0.28/d	-0.33/d	-0.42+i0.04/d	-0.42-i0.04/d
3	-0.01/d	-0.28/d	-0.33/d	-0.41+i0.04/d	-0.41-i0.04/d

Figure 7 compares the dynamics of the amplitude during stage 1 for the two waves. The amplitude as plotted in panel A of Fig. 5 is shown. In addition, 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 as exponentially increasing functions (gray dotted lines) with exponential coefficients as given in Tables 1 and 2 for stage 1. The exponential approximations fit the numerically obtained (exact) solutions very well.

Fig. 7

Dynamics of the dominant amplitudes during the outbreak stages (stages 1) of the first and second-waves hitting the state of New York. Solid black lines corresponds solutions computed from Eq. (30), whereas dotted gray lines stand for exponential approximations using the maximal eigenvalues listed in Tables 1 and 2

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 components of susceptible individuals (S), exposed individuals (E), symptomatic infectious individuals (), asymptomatic infectious individuals (), hospitalized infectious individuals (), and infectious individuals who are treated in intensive care units (). 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 given byThe model parameters correspond to the proportion of asymptomatic infections (), incubation period (), 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 of individuals out of certain compartments. The rate constant 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 . Again, the variable can be neglected because it does not play a role for the stability of the disease-free fixed point . 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 disease-free 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 of cumulative COVID-19 cases was introduced. From Eq. (33), it follows that satisfies [17]With the help of , the model parameters and in particular were fitted to confirmed COVID-19 cases of Pakistan [17, 40]. The best-fit parameters for the three-stage 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 as obtained by solving Eqs. (33) and (35) for the best-fit parameters. Moreover, the linearization matrix readsThe eigenvalues of 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).

Fig. 8

Model-based analysis of confirmed COVID-19 cases of the first-wave epidemic in Pakistan, 2020. A Cumulative COVID-19 cases during the period from March 1 to September 30 (circles) and model fit (solid line) computed from Eqs. (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 of stage 1 and its stage-3 remnant . In all four subpanels, the lower and upper directions refer to and , respectively

Table 3

Eigenvalues of the three stages of the first COVID-19 wave in Pakistan, 2020, as obtained for the infected five-dimensional subspace

	Rank
Stage	1	2	3	4	5
1	0.05/d	-0.49/d	-0.81/d	-1.31+i0.18/d	-1.31-i0.18/d
2	0/d	-0.51/d	-0.78/d	-1.10+i0.14/d	-1.10-i0.14/d
3	-0.03/d	-0.53/d	-0.75/d	-1.10+i0.10/d	-1.10-i0.10/d

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 (order parameter) . Likewise, in stage 3 there was one real-valued eigenvalue with a maximal real part. The corresponding stable eigenvector can be interpreted as the remnant of the unstable eigenvector of stage 1. Panel B of Fig.  8 shows the phase curves of the COVID-19 epidemic as obtained from the solutions of Eq. (33) in the two-dimensional subspaces , and . The eigenvector for stages 1 and 3 is shown there as well. As can be seen, the epidemic followed closely in both stages 1 and 3.

The functions were substituted into Eq. (29) to obtain the amplitudes 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 . Panels B and C show the real- and imaginary-parts of the amplitudes and related to the two remaining eigenvalues that were complex-valued in all three stages. Comparing the amplitudes and with , it is clear that 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 . In contrast, while 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 to (which again is consistent with the fact that was the amplitude in stage 2 that exhibited a zero eigenvalue).

Fig. 9

Amplitude description of the first-wave epidemic in the Pakistan, 2020. A Amplitudes (thick solid line), (dotted line), and (thin solid line) as functions of time in all three stages S1, S2, and S3. B and C Real (B) and imaginary (C) parts of (solid line) and (full circles) in stages S1, S2, and S3. The real parts shown in B are identical in all stages. Likewise, the imaginary parts shown in C are mirror images of each other in all stages

Discussion

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 is known to describe the emergence of convection rolls in fluid and gas systems heated from below [8, 44–46]. The emergence of Turing patterns in chemical systems [9, 47–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–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 epidemic, 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-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 , what is the corresponding increase in infectious individuals? If we distinguish between asymptomatic and symptomatic infectious individuals , we may ask, how are the changes and related to ? 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.

Concluding remarks

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, exponential-like 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 general, and for the time course of the function shown in Fig. 10, in particular. As far as the epidemic situation in December 2021 is concerned, the impact of the omicron variant of the SARS coronavirus 2 is difficult to predict. 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.

Fig. 10

Death data of the third epidemic outbreak in the state of New York, during the winter period from August 1 to November 31, 2021. Cumulative deaths (circles) are shown. The three months period from September to November is fitted by a regression line (solid line) to highlight the apparent straightness of the function of observed cumulative deaths during that period