1 Introduction

The first patients with coronavirus disease 2019 (COVID-19) were reported from Wuhan city, China, in December 2019 [55]. Within one year, the disease spread over the whole globe and by January 2021 there were worldwide 93,000,000 confirmed COVID-19 cases and 2,000,000 COVID-19 associated deaths [56]. About half a year later, by July 2021, the death toll doubled and was at 4,000,000 people [57]. By the time of writing this paper, the pandemic is still unfolding.

Understanding the principles of nonlinear physics that determine how COVID-19 spreads in populations and how the virus spreads and replicates in the bodies of patients is a high-priority task that is far from being completed. In this study, the nonlinear physics perspective is taken that diseases emerge from instabilities as suggested by Mackey and Glass [17, 39]. Such instabilities traditionally are investigated by means of eigenvalue analysis that comes with fixed eigenvalues and eigenvectors [17, 32, 35]. Eigenvalue analysis that takes the time course of eigenvalues into account has also been turned out to be a fruitful tool in economy to investigate industrial processes [31, 38, 46, 49].

In the wake of the COVID-19 pandemic, various modeling approaches have been developed and applied to obtain insights into the nature and the possible time course of the pandemic [7]. The modeling methods and objectives are diverse. For example, various modeling studies have focused on modeling the impacts of intervention measures. Other studies have been concerned with the social and economic impacts of the pandemic [7]. One of the many modeling gaps in the field of COVID-19 research is the gap in understanding the nature of the virus dynamics within individual humans [7]. In order to close this gap and achieve such an understanding, the viral load dynamics of COVID-19 patients has been studied in previous works with the help of a standard, three-variable TIV model [24, 37, 52]. Higher-dimensional models beyond the three-variable TIV model have also been used, for example, to study SARS-CoV-2 infected cells in the eclipse phase [30, 37, 52] and the relevance of malfunctioning, non-infectious SARS-CoV-2 particles for the disease progression of COVID-19 [11, 43]. In particular, the aforementioned, standard TIV model was used to transfer a key finding from the epidemiology of COVID-19 to the virology of SARS-CoV-2 infections. Accordingly, previous work has suggested that COVID-19 spreads out in populations along unstable eigenvectors [20,21,22,23, 25]. A first attempt to show a similar result on the level of individual COVID-19 patients was made in Ref. [24]. However, Ref. [24] was based on a challenging three-dimensional vector analysis approach. Moreover, the study missed to take advantage of a time-resolved eigenvalue analysis as used in economic theory [31, 38, 46, 49]. The current study shows that a simplified, two-dimensional analysis is sufficient to demonstrate that the viral infections of COVID-19 patients are instability-induced phenomena, whose dynamics is determined by the eigenvalues and eigenvectors of the instabilities that give rise to COVID-19. Furthermore, the current study shows that the time-resolved eigenvalue analysis suggests that the switching of the sign of the leading eigenvalue indicates a transition from the stage in which virus load increases in patients to the stage in which the viral load decreases in patients. In fact, on the population level it has been shown that such sign-switching phenomena indicate that COVID-19 epidemics enter subsiding stages, in which the number of daily new COVID-19 infections decrease from day to day [20, 22, 23, 25].

In summary, the key contributions of this work are as follows:

  • In the current study, an analytical method is developed to identify the dynamics of coronavirus diseases 2019 in a suitably defined two-dimensional space during the initial stage of the disease. In this context, the novel aspect is to consider COVID-19 from the perspective of dynamical diseases and, accordingly, to focus on the unstable eigenvector that based on theoretical considerations should determine the initial disease stage. A simplified methodology is developed that works on the level of individual patients and holds the promise to produce robust results irrespective of patient-specific differences (Sects. 2.1 and 2.2).

  • The current work transfers the method of time-resolved eigenvalue analysis from the field of economics to the field of virology. In doing so, a novel mathematical framework is established that allows to observe and study the sign-switching phenomenon on the level of COVID-19 patients (Sect. 2.3). So far, this phenomenon has only be observed on the population level in the context of COVID-19 waves.

  • The proposed methods are applied to real data of a sample of eight COVID-19 patients. The study shows that despite individual differences between patients the proposed methods yield the theoretically predicted results. Accordingly, it is found that for all patients the initial stage COVID-19 dynamics followed patient-specific unstable eigenvectors and for all patients the disease progressions exhibited the sign-switching phenomenon (Sect. 3).

2 Model-based fixed and time-resolved eigenvalue analysis

2.1 Basics

The following standard three-variable virus dynamics model will be used [34, 45, 47]:

$$\begin{aligned} \frac{\text{ d }}{\text{ d }t}T = - \beta V T \ , \ \frac{\text{ d }}{\text{ d }t}I = \beta V T - k_1 I \ , \ \frac{\text{ d }}{\text{ d }t}V = p I - k_2 V \ . \end{aligned}$$
(1)

Here, T, I, and V denote non-infected target cells, infected target cells, and the free virus concentration, respectively. The model parameters are [34, 45] the infection rate \(\beta \), the production rate p of viruses by infected cells, the death rate \(k_1\) of infected cells, and the virus clearance rate \(k_2\). For sake of brevity, T and I will be simply referred to as target cells and infected cells. Let us identify the variables T, I, and V in the context of SARS-CoV-2 infections. Studies that have examined the SARS coronavirus outbreak in 2002–2004 have argued that the SARS coronavirus attacks epithelial cells in the lungs [8, 33]. Given that SARS coronavirus 2 is related to the earlier (2002–2004) SARS coronavirus, various authors have assumed that SARS-CoV-2 invades respiratory epithelial cells, in general [41, 50], and, in particular, in the lungs [44, 48, 59, 61]. In fact, some case studies on COVID-19 patients have produced evidence that supports this hypothesis [9, 40, 58]. Consequently, it is preliminarily assumed that the cells T of Eq. (1) corresponds to epithelial cells in the lungs. The virus load V reflects SARS-CoV-2 load measured in the sputum of COVID-19 patients [5, 53].

2.2 Fixed eigenvalue analysis using a simplified, two-variable approach

Let \((T_{st},0,0)\) denote the virus-free fixed point, where \(T_{st}>0\) is the target cell concentration of healthy adults. Linearizing Eq. (1) at \((T_{st},0,0)\), we obtain the \(3 \times 3\) linearization matrix

$$\begin{aligned} U= \left( \begin{array}{ccc} 0 &{} 0 &{} -\beta T_{st}\\ 0 &{}-k_1 &{} \beta T_{st}\\ 0&{} p &{} -k_2 \end{array} \right) \ . \end{aligned}$$
(2)

From U it follows that the fixed point \((T_{st},0,0)\) exhibits a zero eigenvalue \(\lambda _1=0\) and two eigenvalues \(\lambda _{2,3}\) that in general are finite. The latter eigenvalues can be computed from the dynamics in the I–V subspace. The linearized dynamics in this I–V subspace reads

$$\begin{aligned} \frac{\text{ d }}{\text{ d }t}I = \beta V T_{st} - k_1 I \ , \ \frac{\text{ d }}{\text{ d }t}V = p I - k_2 V \end{aligned}$$
(3)

and describes a closed dynamical system. In summary, the original three-variable model exhibits a neutrally stable direction which is given by the T axis. That is, \(T_{st}\) can assume any positive value. Regardless of \(T_{st}\) the eigenvalue \(\lambda _1\) equals zero. In contrast, the dynamics in the I–V subspace can exhibit a stable or unstable fixed point \(I=V=0\). Consequently, the subspace dynamics defined by I and V determines the overall stability of the fixed point \((T_{st},0,0)\) of the three-variable model. For this reason, a simplified, two-variable approach is motivated that focuses on that subspace dynamics.

Let us introduce the state vector \(\mathbf{X}=(I,V)\) with fixed point \(\mathbf{X}_{st}=(0,0)\). Then, Eq. (3) can be cast into the form

$$\begin{aligned} \frac{\text{ d }}{\text{ d }t}\mathbf{X}= L \, \mathbf{X} \ , \ L= \left( \begin{array}{cc} -k_1 &{} \beta T_{st}\\ p &{} -k_2 \end{array} \right) \end{aligned}$$
(4)

such that L is the linearization matrix. From Eq. (4) the eigenvalues \(\lambda _{2,3}\) can be obtained and read

$$\begin{aligned} \lambda _{2,3} = -\frac{k_1+k_2}{2} \pm \sqrt{ \frac{(k_1+k_2)^2}{4}+ p T_{st} \beta - k_1 k_2 } \end{aligned}$$
(5)

with and the plus (minus) sign holds for \(\lambda _2\) (\(\lambda _3\)). From Eq. (5) it follows that \(\lambda _3<0\) holds in any case. Equation (5) allows to distinguish between the cases \(\lambda _2>0\), \(\lambda _2=0\), and \(\lambda _2<0\) that hold for \(p T_{st} \beta >k_1k_2\), \(p T_{st} \beta =k_1k_2\), and \(p T_{st} \beta <k_1k_2\). Consequently, the fixed point \((T_{st},0,0)\) is unstable for \(p T_{st} \beta >k_1k_2\) and neutrally stable (in view of \(\lambda _1=0\)) for \(p T_{st} \beta <k_1k_2\).

Let \(\mathbf{v}_j=(v_{j,I},v_{j,V})\) denote the eigenvector j of L with \(j=2,3\). The eigenvector description of the viral load dynamics in the two-dimensional I–V subspace can be written like [17, 35, 49]

$$\begin{aligned} \mathbf{X}(t)=\sum _{j=2,3} A_j(t) \mathbf{v}_j \ , \mathbf{v}_j = \frac{1}{Z_j} \left( \begin{array}{c} \beta T_{st} \\ \lambda _j +k_1 \end{array} \right) \ , \end{aligned}$$
(6)

where \(Z_j\) are normalization constants and \(A_j\) denote time-dependent amplitudes. In the initial stage, when the virus multiplies in the human lung, we have \(\lambda _2>0\) versus \(\lambda _3<0\) and, consequently, \(\mathbf{v}_2\) is the unstable eigenvector, whereas \(\mathbf{v}_3\) is the stable eigenvector. Moreover, the infection dynamics follows primarily the direction given by \(\mathbf{v}_2\). Within the modeling framework given by Eq. (1), \(\mathbf{v}_2\) is the order parameter [17, 35] and characterizes together with the exponentially increasing amplitude

$$\begin{aligned} A_2(t)=A_2(0) \exp \{\lambda _2 (t-t_0)\} \end{aligned}$$
(7)

the SARS-CoV-2 infection in the human lung. In particular, after a transient period in which \(A_3\) decays exponentially to zero, we have

$$\begin{aligned} \Delta \mathbf{X} = \mathbf{v}_2 \Delta A_2 \end{aligned}$$
(8)

with \(\Delta \mathbf{X} = \mathbf{X}(t+\Delta t)-\mathbf{X}(t)\) and \(\Delta A_2 = A_2(t+\Delta t)-A_2(t)\). For more general discussions of the order parameter amplitude \(A_2\) and Eq. (8) see Refs. [17, 35].

2.3 Time-resolved eigenvalue analysis

The eigenvalue \(\lambda _2\) of the unstable direction at the instability fixed point \((T_{st},0,0)\) is given by Eq. (5) and (by definition of the instability) positive: \(\lambda _2>0\). In general, any state (T, 0, 0) with \(T \le T_{st}\) is a fixed point of Eq. (1) (whereas any state (T, I, V) with \(I>0\) or \(V>0\) is not a fixed point of Eq. (1)). Consequently, when during the course of the disease the number of target cells decays and is given at a particular time point t by T(t), then the state (T(t), 0, 0) is the special fixed point that is consistent with the number of target cells T(t) at the time t under consideration. We will refer to (T(t), 0, 0) as a sliding fixed point. A time-resolved eigenvalue analysis can be obtained by taking into account that the target cells T decay over time and considering the stability of the sliding fixed point (T(t), 0, 0). The stability changes when the number of target cells T is sufficiently low such that \(pT\beta < k_1 k_2\) (see Sect. 2.2). When the fixed point becomes stable, the infection subsides in the sense that the virus load V begins to decay. Therefore, a time-resolved eigenvalue \(\lambda _2\) can capture that important switch from an unstable to a stable virus-free fixed point. Accordingly, in order to conduct a time-resolved eigenvalue analysis, in Eq. (5) for \(\lambda _2\) the parameter \(T_{st}\) is replaced by T(t) such that

$$\begin{aligned} \lambda _2(t) = -\frac{k_1+k_2}{2} +\sqrt{ \frac{(k_1+k_2)^2}{4}+ p T(t) \beta - k_1 k_2 } \ . \end{aligned}$$
(9)

The model solution (T(t), I(t), V(t)) connects two fixed points defined by \((T(0),I(0),V(0))=(T_{st},0,0)\), on the one hand, and \((T(\infty ),I(\infty ),V(\infty ))=(T(\infty ),0,0)\) with \(p T(\infty )\beta < k_1 k_2\), on the other hand. While the first fixed point is unstable, the second fixed point is stable. They reflect the pre-infection and post-infection fixed points of the patients under consideration.

2.4 Application to patient data

The fixed and time-resolved eigenvalue analysis was applied to data from eight COVID-19 patients described in Refs. [5, 53]. Following the notation of Ref. [53], the patients are referred to as patients 1, 2, 3, 4, 7, 8, 10 and 14. The patients developed symptoms during the period of January and February 2020. During that time they were located in Germany. The disease was considered to be mild for all patients. SARS-CoV-2 viral load was measured from the sputum for all patients for periods of about 25 days [5, 53]. On the basis of the (lower respiratory tract) sputum data, the model parameters \(\beta \), p, \(k_1\), and \(k_2\) of Eq. (1) were estimated and reported in Ref. [52]. More precisely, in Ref. [52] several virus dynamics models were considered including the TIV model described by Eq. (1). Moreover, in Ref. [52] viral load data from the upper and lower respiratory tracts were analyzed. In the current study, the estimates \(\beta \), p, \(k_1\), and \(k_2\) presented in Ref. [52] for the TIV model (1) applied to the lower respiratory tract data were considered. Using those estimates, the fixed eigenvalue analysis was conducted for every patient by computing the fixed eigenvalue \(\lambda _2\) from Eq. (5) and the fixed unstable eigenvector \(\mathbf{v}_2\) from Eq. (6). Subsequently, the infection dynamics computed from Eq. (1) was plotted in the I–V subplane and compared with \(\mathbf{v}_2\) from Eq. (6). The time-resolved eigenvalue analysis was conducted for every patient by computing \(\lambda _2(t)\) from Eq. (9) using Eq. (1) to determine T(t). The time-course of \(\lambda _2\) was compared to the time-course of the viral load given in terms of V(t).

3 Results and discussion

Figure 1 shows the results of the eigenvalue analysis for patient 1. Panel A shows the measured viral load over time (gray circles) and the model fit (solid line) obtained by computing V(t) from Eq. (1) numerically. Time is given by days after onset of symptoms [53]. The data (gray circles) reveal that the viral load of patient 1 reached a peak at about 5 and 6 days after symptoms onset. Subsequently, the viral load decayed in a more or less monotonically fashion over a 20 days period. From the data it follows that the initial, increasing stage was relatively short as compared to the final, decaying state. The model solution reflects this property and exhibits a relatively quickly increasing viral load dynamics and a relatively slowly decreasing dynamics. Panel B shows the functions V(t) and I(t) as computed from the model (1) as a phase curve (solid thin line) in the I–V space. The circle indicates the initial state. The order parameter \(\mathbf{v}_2\) is plotted as well (dashed thick line). As can be seen in panel B, the virus dynamics closely followed the order parameter while the viral load was increasing. After the viral load reached its peaks the dynamics branched off in a different direction. Panel C shows the time-dependent eigenvalue \(\lambda _2\) (solid line) as computed from Eq. (9). The virus load curve V(t) as depicted in panel A (solid black line) is shown in panel C as well in as a rescaled function (dotted line). The eigenvalue \(\lambda _2\) dropped from a positive to a negative value at around day 5. Subsequent to this switch, the viral load reached its peak and began to decay. Consequently, for patient 1 there was a first, initial period during which the eigenvalue was positive and the virus-free fixed point unstable. During this period, viral load increased and the number of target cells decayed. There was a switching point at around 5 days when the number of target cells were sufficiently low and the eigenvalue \(\lambda _2\) turned from a positive to a negative value. There was a second, final period during which the eigenvalue \(\lambda _2\) was negative and the virus-free fixed point stable. During this stage the viral load decayed.

Fig. 1
figure 1

Eigenvalue analysis of the viral load dynamics of patient 1. Panel A: Observed viral load over time (gray circles) and model fit V(t) (solid line) computed from Eq. (1). Panel B: Phase space trajectory V(t) versus I(t) (solid thin lines) computed from model (1) and order parameter \(\mathbf{v}_2\) (thick dashed lines) obtained from the eigenvalue analysis based on Eq. (6). The order parameter is magnified in length for sake of visibility. Panel C: Time-course of the eigenvalue \(\lambda _2\) (solid line) as computed from Eq. (9). The viral load trajectory V(t) shown in panel A is depicted in panel C as well (dotted line). For visual inspection purposes, it was rescaled such that its peak corresponds to the maximum value of \(\lambda _2\)

Full size image
Fig. 2
figure 2

As in Fig. 1 but for patient 2

Full size image
Fig. 3
figure 3

Viral load dynamics observed (gray circles) and computed (solid black lines) for the six patients 3, 4, 5, 8, 10, and 14

Full size image

Figure 2 presents the same analysis for patient 2. Again, the virus dynamics was characterized by a short, initial stage during which viral load increased and a relatively long, final stage during which viral load decreased (panel A). The phase portrait (panel B) reveals that during the initial stage the dynamics followed the order parameter \(\mathbf{v}_2\) and subsequently branched off. The two stages can be explained using the time-resolved eigenvalue analysis. Accordingly, during the initial stage the eigenvalue \(\lambda _2\) was positive, while in the second, decaying stage the eigenvalue was negative (panel C). The switch of the eigenvalue was caused by the decrease in the number of target cells, see Eq. (9).

Figure 3 shows the viral load trajectories for the remaining six patients labeled 3, 4, 7, 8, 10 and 14 [53]. The measured data (gray circles) [53] and computed model solutions V(t) (solid lines) are presented. Panels A to F correspond to patients 3, 4, 7, 8, 10, and 14 as indicated. All patients showed the same characteristic temporal pattern composed of a relatively short increase and a slow decay of viral load.

Fig. 4
figure 4

Panels A, C, E: Phase trajectories computed from Eq. (1) for patients 3, 4, and 7, respectively. Order parameters \(\mathbf{v}_2\) (thick dashed lines) computed from Eq. (6). Panels B, D, F: Graphs \(\lambda _2(t)\) (solid lines) obtained from Eq. (9). For comparison purposes, the corresponding viral load trajectories V(t) are shown as well (dotted lines) as functions rescaled to the maximum values of \(\lambda _2\)

Full size image
Fig. 5
figure 5

As in Fig. 4 but for patients 8, 10, and 14

Full size image

Figures 4 and 5 present the eigenvalue analyses for the aforementioned remaining six patients in a shortened form. Figure 4 presents the analyses for patients 3, 4, and 7. Panels A, C, E and B, D, F show the results of the fixed eigenvalue analysis and time-resolved eigenvalue analysis, respectively. Likewise, Fig. 5 presents the analyses for patients 8, 10, and 14. Again, panels A, C, E and B, D, F show the results of the fixed eigenvalue analysis and time-resolved eigenvalue analysis, respectively.

As can be seen in panels A, C, E of Figs. 4 and 5, for all remaining patients the virus dynamics (solid thin line) in the initial stage of increasing viral load followed the respective SARS-CoV-2 order parameters \(\mathbf{v}_2\) (dashed thick lines). After reaching the respective peak loads, the trajectories branched off from their order parameters. Consequently, the results of those patients are consistent with those of patients 1 and 2. As can be seen in panels B, D, F of Figs. 4 and 5, the time-resolved eigenvalue analysis shows that for all but two patients the switch of the leading eigenvalue from a positive to a negative value took place after at least 1 day after symptoms onset. Accordingly, the initial stage of increasing viral load was at least 1 day long. The two exceptions were patients 4 and 8. The respective eigenvalues switched at about 0.10 days (i.e., about two-and-a-half hours) and 0.15 days (i.e., about three-and-a-half hours). In order to make these very early switching dynamics visible, panel D of Fig. 4 (patient 4) and panel B of Fig. 5 (patient 8) only show 1 day on the horizontal axis. For all patients (i.e., including patients 4 and 8) the initial stages of increasing viral load and the final stages of decreasing viral load were characterized by a positive and negative eigenvalue \(\lambda _2\), respectively. This is tantamount to say that these two stages were characterized by unstable and stable virus-free fixed points, respectively.

4 General discussion

Eigenvalue analysis is an indispensable tool for studying instabilities in physics and other disciplines. For example, instabilities in fluid and gas dynamical systems and optical systems leading to convection rolls and laser light, respectively, can be understood by determining the conditions under which relevant eigenvalues become positive (or exhibit positive real parts) [10, 19, 35]. Stripe and dot patterns emerging on animal skins [42], brain activity patterns [6, 18, 26, 27], and Turing patterns in chemical and biological systems [1, 12, 15, 16, 28, 29] can be explained with the help of eigenvalue analysis that determines the circumstances under which positive eigenvalues or eigenvalues with positive real parts occur. Infectious diseases emerging in populations are long known as phenomena that are due to instabilities featuring positive eigenvalues [14]. For the COVID-19 outbreak in Wuhan city, China, and for first-wave COVID-19 epidemics in European countries in 2020, in Pakistan, Thailand, and USA this has been shown explicitly in a series of recent studies [20,21,22,23, 25].

However, while eigenvalues may be considered as constants in various settings, in general, they may vary over time or as a result of a change in the state of the system under consideration. This kind of time-dependency and elasticity has been studied for various industrial processes [31, 38, 46, 49]. Likewise, it has been argued that COVID-19 waves are the product of such changing eigenvalues. More precisely, it has been argued that intervention measures during the year 2020 such as mask-wearing policies, physical distancing requirements, border closures, and the lockdown of businesses have affected eigenvalues and resulted in a switching phenomenon that turned positive eigenvalues into negative ones [20,21,22,23, 25]. Instabilities on the level of populations leading to COVID-19 outbreaks turned into stable fixed points under which the number of daily new COVID-19 infections decreased from day to day. In fact, the notion of eigenvalues that depend on states or time has been formalized in the field of human performance and perception as grand state dynamics [17]. A plenitude of phenomena have been explained in this context such as: scene decomposition, lonely speaker dynamics, Freud’s slip of the tongue phenomenon, child play, retrieval-induced forgetting, priming, overcoming of functional fixedness, oscillatory brain activity induced by visual forces with certain symmetry features, learning, negative hysteretic transitions in human performance, motion-induced blindness, rituals of obsessive-compulsive-disorder patients, social talking, and life trajectories [17]. The current study adds to these examples of systems exhibiting time-varying eigenvalues two more systems: the virus dynamics in human individuals, in general, and the SARS-CoV-2 dynamics in COVID-19 patients, in particular. While there are peculiarities of viral infections that need to be studied in detail, it is also important to see the fundamental underlying principles of nonlinear physics that determine not also SARS-CoV-2 infections but instability-induced phenomena in general. The current study makes a step into this direction and provides a general framework to study SARS-CoV-2 infections in the context of such generally applicable principles.

In order to demonstrate the applicability and robustness of the proposed eigenvector and eigenvalue analysis methods, the current study looked at a sample of eight COVID-19 patients. COVID-19 data from patients has been published in various studies (see, e.g. Ref. [51] for a review). Variability across patients as shown when comparing patients 1, 2, 3, 7, 10, and 14 versus patients 4 and 8 (see Figs. 1, 2 and 4) has been observed in other studies as well (e.g., see Ref. [36]). Despite the variability across patients, SARS-CoV-2 infections seem to exhibit a generic pattern consistent across patients [51]: after the onset of symptoms, viral load increases quickly (i.e., within a few days) towards a peak value. Subsequently, viral load decays slowly (i.e., within ten days or more) until virus measurements fall below the detection threshold. Data as shown in Figs. 1, 2 and 3 are consistent with this general temporal pattern. In other words, the sample of eight viral load trajectories can be considered as a representative data sample that exhibits, on the one hand, a common pattern and, on the other hand, variations of that pattern. Despite these variations, the proposed methods produced for all patients the theoretically predicted results. Accordingly, it was found that for all patients the initial stage COVID-19 dynamics followed the patient-specific unstable eigenvectors (see Figs. 1B, 2B, 4A,C,E, and 5A,C,E) and for all patients the disease progressions exhibited the sign-switching phenomenon (see Figs. 1C, 2C, 4B, D, F, and 5B, D, F). Having said that, future studies may be devoted to apply the methods illustrated in Figs. 1, 2, 4, and 5 to larger samples. In particular, the question raised in Ref. [7] about the nature of the difference between mild and severe infections may be addressed by studying such larger samples.

The model-based analysis presented in the current study involved a standard three-variable virus dynamics model used in the literature [34, 45, 47, 52]. In this context, it should be mentioned that there are various generalizations of the three-variable model leading to models in higher-dimensional spaces [41, 45, 47, 52]. Generalized TIV models involving more than three variables have been used to describe the viral load dynamics of COVID-19 patients. The motivation to use such generalized models typically is to address specific aspects of the progression of COVID-19. For example, as pointed out in the introduction, the relevance of infected cells in the eclipse phase [52] and malfunctioning, non-infectious virus particles [43] has been addressed in the context of COVID-19. Importantly, the primary goal of studies using generalized TIV models typically is to address such specific aspects of virus dynamics. In other words, the primary motivation to use models beyond the TIV model is not to achieve better fits between models and data [2, 3]. Consequently, although some details in the time course of the viral load are not addressed by the three-variable model (1), the model can address very well the general temporal pattern of SARS-CoV-2 infections. In view of the plenitude of models that could be used to address SARS-CoV-2 infections, the results presented in this study should be regarded as a step in the direction of understanding SARS-CoV-2 infection instabilities by means of fixed and time-resolved eigenvalue analysis approaches. The eigenvalue analysis presented in this study may be used in future studies as a tool to address aspects of the virus dynamics in COVID-19 patients that go beyond those aspects captured by the TIV model (1).

Fig. 6
figure 6

Impact of a SARS-CoV-2 mutation that exhibits an increased infectivity rate \(\beta \) on the critical target cell concentration \(T_{crit}\). The graph was computed from \(T_{crit}=k_1 k_2/(p \beta )\) with parameter values taken from patient 1. \(T_{crit}\) is shown on the vertical axis as \(T_{crit,rel}\) as percentage value of the pre-infection level \(T_{st}\)

Full size image

Since the beginning of the pandemic, several variants of the original SARS coronavirus 2 have occurred. In particular, variants such as Alpha, Delta, and Omicron exhibit a mutation on the spike protein of the virus. Accordingly, at a particular location (which is labeled as 614) the Adenine nucleotide is replaced by Guanine [4]. It has been suggested that due to this so-called D614G mutation, variants such as Alpha, Delta, and Omicron can spread out faster in populations [4]. On the level of individuals, the D614G mutation seems to increase the infectivity rate [4, 60], which corresponds to the parameter \(\beta \) in the TIV model (1). At the same time, when comparing the Delta and Omicron variants, preliminary results suggest that the latter variant is associated with less severe diseases progressions [54]. Let us briefly indicate how such issues could be addressed in the context of the eigenvalue analysis discussed in Sect. 2.3. From Eq. (9) it follows that the eigenvalue \(\lambda _2\) becomes zero at a critical target cell concentration \(T_{crit}=k_1 k_2/(p \beta )\). Consequently, in mutations with increased infectivity parameters \(\beta \), the critical value is lower which implies that there is larger damage in the affected sites before the disease enters the decline stage. Figure 6 illustrates this issue. All model parameters are taken from patient 1 except for \(\beta \). The parameter \(\beta \) is varied from the estimated value to a value that is two times as high. Accordingly, \(\beta _{rel}=100\%\) reflects the estimated value (as used to produce Fig. 1), while \(\beta _{rel}=200\%\) in Fig. 6 reflects an infectivity rate that is twice as high as this estimated value of patient 1. In Fig. 6, the critical value of target cells at which the disease enters the decline stage is expressed as \(T_{crit,rel}\) in percentage of the pre-infection level \(T_{st}\). For the estimated value \(\beta \) of patient 1 (i.e., \(\beta _{rel}=100\%\)), the disease decline begins when T reaches about 2.4% of \(T_{st}\). When \(\beta \) is increased, this percentage value decays as shown in Fig. 6. In doing so, Fig. 6 illustrates that the severity of the disease increases when the infectivity rate due to mutations increases. In fact, in the context of the Alpha variant, it has been suggested that the disease severity of the Alpha variant is increased as compared to the original SARS coronavirus 2 [13].

Fig. 7
figure 7

Comparing mutations that vary with respect to the clearance parameter \(k_2\) on the basis of \(T_{crit,rel}\). Mutations X and Y indicate mutations that exhibit the lowest and highest values of \(k_2\), respectively, considered in the comparison. The graph was computed from \(T_{crit}=k_1 k_2/(p \beta )\) with parameter values taken from patient 1. See text for details

Full size image

Figure 7 demonstrates a mutation-related mechanism leading to a lower chance of the occurrence of severe disease cases. As reviewed above, preliminary data suggest that the Omicron mutation may induce primarily mild cases of COVID-19 as compared to the Delta mutation [54]. This would be consistent with assuming that the clearance parameter \(k_2\) is larger for the Omicron variant. In order to demonstrate the effect of a mutation-related increased \(k_2\) value, the parameters of patient 1 were used again. \(k_2\) was varied from the estimated value of patient 1 to a value 5 times larger than the estimated one. Accordingly, \(k_{2,rel}\) ranged from 100% to 500%. The relative critical target cell concentration as defined above, was computed as a function of \(k_{2,rel}\).

Figure 7 shows that the relative critical value \(T_{crit,rel}\) is a linearly increasing function of \(k_{2,rel}\), which also can be deduced from the aforementioned equation \(T_{crit}=k_1 k_2/(p \beta )\). High levels of \(T_{crit,rel}\) imply that the disease causes less damage in the human body before it enters the decline stage. Consequently, Fig. 7 illustrates how a mutation-related increase of \(k_2\) may result in a less severe progression of COVID-19. Having said that, more detailed discussions based on patient data are needed to obtain a clear picture about the virus dynamics of variants of the SARS coronavirus 2. In this context, Figs. 6 and 7 illustrate that the framework discussed in the current study may be applied.

5 Conclusion

Although eigenvalues analysis is a frequently used tool in various disciplines, its full scale application to examine viral load data of COVID-19 patients has not yet been discussed in the existing COVID-19 research literature. This paper proposes to take advantage of the concepts of unstable eigenvectors, on the one hand, and time-dependent eigenvalues, on the other hand. Disease-related unstable eigenvectors and time-dependent eigenvalues can be identified in COVID-19 patients using a model-based data analysis approach. That is, in order to conduct the proposed approach a combination of modeling and data analysis is required. The paper demonstrates that it is the unstable eigenvector that determines the initial stage of the disease in a patient and it is the sign-switching of the time-dependent associated eigenvalue that indicates the beginning of the decline of the disease in the patient. In short, the conclusion can be drawn that eigenvalue analysis in its comprehensive way as it is presented in the current study helps to understand both the initial and the final stages of the disease. As a promising result, the analysis of data from a sample of eight patients suggests that differences in disease progression between patients are picked up by the aforementioned time-resolved measure, that is, the time-dependent eigenvalue. Therefore, the suggested approach may open a new avenue for personalized medicine, in particular, in the field of COVID-19.