Social stress drives the multi-wave dynamics of COVID-19 outbreaks

The dynamics of epidemics depend on how people's behavior changes during an outbreak. At the beginning of the epidemic, people do not know about the virus, then, after the outbreak of epidemics and alarm, they begin to comply with the restrictions and the spreading of epidemics may decline. Over time, some people get tired/frustrated by the restrictions and stop following them (exhaustion), especially if the number of new cases drops down. After resting for a while, they can follow the restrictions again. But during this pause the second wave can come and become even stronger then the first one. Studies based on SIR models do not predict the observed quick exit from the first wave of epidemics. Social dynamics should be considered. The appearance of the second wave also depends on social factors. Many generalizations of the SIR model have been developed that take into account the weakening of immunity over time, the evolution of the virus, vaccination and other medical and biological details. However, these more sophisticated models do not explain the apparent differences in outbreak profiles between countries with different intrinsic socio-cultural features. In our work, a system of models of the COVID-19 pandemic is proposed, combining the dynamics of social stress with classical epidemic models. Social stress is described by the tools of sociophysics. The combination of a dynamic SIR-type model with the classical triad of stages of the general adaptation syndrome, alarm-resistance-exhaustion, makes it possible to describe with high accuracy the available statistical data for 13 countries. The sets of kinetic constants corresponding to optimal fit of model to data were found. These constants characterize the ability of society to mobilize efforts against epidemics and maintain this concentration over time and can further help in the development of management strategies specific to a particular society.

www.nature.com/scientificreports/ more rapidly than the SIR-type models prescribe, and then a second wave emerges that is impossible in SIR models. In a recent work 15 , an eight-dimensional SEIR-type model was proposed to analyze the ways to reduce the rate of SARS-CoV-2 virus transmission. Based on its analysis the double-peak behavior caused by removing restrictions was demonstrated 3 . It was also shown that the second peak could be much higher than the first one.
Other studies 16,17 focused directly on the dynamics of the second wave. They performed a comparative analysis of trajectories for different regions and identified the countries that most significantly reduced the mortality rate.
Various stress factors can significantly modify social dynamics 18 . In the modern world, TV, the Internet and social networks provide instant transmission of information. Information spreading in times of crisis affects human behavior and may create "infodemics" 19 . Changes in people's daily lives, either due to their own fears or due to forced administrative restrictions, lead to a significant change in the effective reproduction number of a viral disease. In a recent research 6 , the entire population was divided into the groups where the persons either obey or ignore social distancing rules. This assumption made it possible to determine the size of these groups, depending on social factors. Oscillatory patterns of repeated indulgences (when the situation seems to be improving) and tightening of restrictions (when the spread of infection becomes high) were found.
General adaptation syndrome (GAS) was discovered by Selye (1936) 20 . He found that "a typical syndrome appears, the symptoms of which are independent of the nature of the damaging agent". Later on, this concept was integrated with Cannon's fight-or-flight response 21 , and the general concept of stress was developed with a wide range of applications in physiology, psychology and beyond [22][23][24] . We will use these results as prototypes for our models.
GAS consists of three phases: a mobilization phase, alarm, a resistance phase, and an exhaustion phase 25 . To simulate the population dynamics within the conditions of aggressive COVID-19 spreading, the susceptible group of the persons can be divided into three subgroups of the individuals according to their behavior modes. The models of behavior in which people are seen as behaving in different ways (modes) at different times and in different contexts is widely used in the game theory and analysis of economic and social behavior 26 . According to the concept of GAS and stress, we consider three modes: (i) ignorant mode (persons living without any restrictions as it was before the pandemic), (ii) resistance mode (conscious persons practicing the social distancing rules to avoid the appeared danger) and (iii) exhaustion mode (the depletion of the person's resources leading to reduction of following social distancing rules). In the first approximation, we exclude the alarm mode that may be considered as a short delay before activation of resistant behavior. In more detailed models, alarm mode can be considered separately, but the more phases we consider, the more unknown constants we have to identify, and the simplest model should be tested first. We consider a loop of transitions between modes: This loop should be closed: after some time persons return from the exhaustion stage into initial mode and are susceptible to the alarm signals again. The probability of infection is much less in the resistant mode than in the ignorant or exhaustion ones.
The inclusion of social stress in the SIR type model gives rise to a new model, which we call SIR SS . It describes the multi-peak dynamics in the COVID-19 outbreaks. The model parameters have been successfully fitted to best match the statistical observations of epidemics in different countries of the world.
The parameters of social dynamics found for different countries demonstrate the differences between their social reactions to pandemic stress. These parameters are: a-morbidity rate, K 2 -exhaustion rate, q-stress response rate, I 0 -initial fraction of infected population. The qualitative picture meets the general expectations. Particularly instructive are the quantitative differences between countries. The fitted models can be used in the development of management strategies specific to a particular society. The modeling method can be directly transferred and applied to regions, cities and sufficiently large subpopulations.

Results
Comparative study concept and main results. We used simulations to fit the total cumulative cases (TCC) of COVID-19 observations in China and 12 other countries with a large number of coronavirus cases: Brazil, Colombia, France, Germany, India, Iran, Israel, Italy, Russia, Spain, UK, and USA. We took the duration of the simulation, which fully covers the first wave of the epidemic, but not too long so that it is not affected by seasonal fluctuations in immunity. An interval of 200 days was chosen as the priority period for simulations. We used the calculation of the R 2 coefficient as a goodness-of-fit measure to perform global optimization in the parameter space (a, K 2 , q, I 0 ). Finally, we have identified the precise values of these parameters, which provide the best fit between the theoretical curve of cumulative cases and the observed one. The results are shown in Table 1. Figure 1 presents the results of modeling. The notation used in the figures: TCC-the total number of cumulative cases (observed); DNC-daily new cases (observed); CC(t)-the predicted values of the cumulative cases; CC′(t)-the predicted daily new cases; the normalized values of CC and CC′ are the fractions of CC and CC′ in the population; S-the susceptible fraction of the population; S = S ign + S res + S exh for three modes of behavior: ignorant, resistance, and exhaustion; I-the infected fraction of the population; R-the removed fraction of the population.
In China, the resistance index q is unprecedentedly higher than in other densely populated countries. This immediately indicates that the highest discipline and state control that we observed in China allowed this country to avoid the second wave of the coronavirus outbreak (Fig. 1).
The first COVID-19 outbreaks in Europe. The general sample of 12 countries was divided into four groups of three countries each. First, we analyzed the large European countries, where the coronavirus out- www.nature.com/scientificreports/ break occurred at the beginning of the pandemic: Italy, Germany, and United Kingdom (Fig. 2). An interval of 200 days was also chosen for the simulations.
The results demonstrate an almost complete match of the cumulative cases, CC(t), for the SIR SS model data and the COVID-19 statistics, TCC, (R 2 > 0.997). The dynamics of the epidemic spread were accompanied by antiphase oscillations of S ign (t) and S res (t) with sustained leakage into the exhausted fraction S exh (t). In Fig. 2, the number of S exh reaches the saturation level after the end of the first epidemic wave. Curves for the infected state, I(t) and DNC(t) = TCC′(t), show the two-wave dynamics. The SIR SS model demonstrates the effect of reduction in morbidity at the peak of the first wave due to the rapid increase in the resistant fraction S res (t). Consequently, the I(t) dependence enters the second wave due to the accumulation of social stress and exhaustion caused by restrictions. This distinguishes our SIR SS model from the classical SIR, characterized by the depletion of potential carriers of the infection right after the first wave. For example, in the simplest SIR model with two parameters (a, b), only one epidemic wave I(t) can exist. The function CC(t) = I(t) + R(t) is a sigmoid with a saturation level of Italy Germany United Kingdom Countries with a high resistance index. Next, we analyzed the infection data of COVID-19 for countries where the coronavirus outbreak evolved dramatically like in China: it was relatively short-lived, and after it, the morbidity rate DNC(t) rapidly dropped to a low level. For this purpose data about France, Israel, and Spain were selected and analyzed (Fig. 3). Several trends were revealed from the results in Table 1. The parameter q, France Israel Spain Figure 3. Coronavirus outbreaks in countries with a high resistance index. Results for France (top panel), Israel (center), and Spain (bottom). The outbreaks dynamics are characterized by a rapid increase in the S res fraction during the peak of the first wave and the formation of protracted plateau on CC(t) with low DNC indices. To simulate the spread of COVID in Israel (population is less than 10 million), we take an interval of 100 days. Numerous data corrections in France and Spain for the period April'20-May'20 lead to the minor (practically negligible) discrepancies between the CC and TCC profiles. www.nature.com/scientificreports/ which is responsible for the mobilization response, for France and Israel is several times higher than the median value q* = 85.2. At the same time, in France and Spain, the infection rate, a, takes on abnormally high values: 0.3138 and 0.3770 relative to the median, a* = 0.2065. Thus, the tendency of society to trigger into the alarm state, combined with a high infection rate, leads to highly unified society's response 24 . Note, that at the peak the fraction of actively resistant people S res can rise up to 90% of total population.

Countries with a low resistance index.
On the contrary, when the parameter q takes small values, the behavioral reaction of the population to the growth of I(t) is weakly expressed. This case corresponds to the third group of countries selected for analysis: Brazil, India, and Russia ( Fig. 4). At q ≪ q*, the similarity of the SIR SS model with the classical SIR is high. Combined with the low infection rate, it leads to the "lengthening" of the first epidemic wave in the timeline. The fraction of resistant population, S res (t), over the entire interval of the outbreak for these countries is less than 40%. The K 2 parameter for India, which determines the predisposition to trigger into exhaustion mode, is close to this parameter for all five analyzed European countries. Nevertheless, in India this did not prevent S exh (t) from evolving at the lowest level among the countries of the general sample, not exceeding 20%, due to the constantly low level of the fraction S res (low resistance − low exhaustion). This parameter for Russia (13.45 × 10 -3 ) and Brazil (15.36 × 10 -3 ) was significantly higher than the median K 2 * = 7.27 × 10 -3 . However, the saturation level of the exhaustion fraction for Russia and Brazil was comparable to Western European countries with minimal K 2 : S exh * ~ 30%.
The distinctive feature of this group was the lower morbidity rate a, which was in the top-3 of the lowest values of a. It required much longer time intervals for simulation (up to 400 days). The persistence of the high fraction S ign , in contrast to the countries of the second group, can be a precondition for the imminent emergence of the second, more substantial epidemic wave without the rest period between waves. Compare the DNC(t) trajectory with a low plateau in Fig. 3 and the trajectories without such a plateau in Fig. 4. This fact corresponds to the official COVID-19 statistics.
Countries with a high exhaustion rate. Next, we analyzed COVID-19 data for the remaining three countries: Colombia, Iran, and USA (Fig. 5). This group was characterized by relatively high values of the exhausted population, S exh ~ 50%. As expected, the parameter K 2 ~ 20 × 10 -3 took extremely high values; this group is in the top-3 countries from the general sample in terms of the K 2 coefficient. Ordinary infection rates a ~ 0.15-0.2 indicated the average growth rate of the cumulative cases number. Large scatter of values of q illustrates the weak dependence of CC(t) trajectory on this parameter.
Similarly to the previous group, the formation of the second epidemic wave was also revealed. However, in this case, the forecast foreshadows a "huge" wave in which the DNC values may be many times greater than the first peak.

Discussion
We demonstrated how the dynamics of the COVID-19 epidemics are mainly driven by the dynamics of social stress. The large difference between epidemics in different countries is caused by social differences and not by real epidemiological (biological) factors, at least until new local viral mutations appear.
The simple SIR SS model works well for 200-300 days of pandemics. We cannot expect that such a simple model will work much longer for various reasons. It has limitations and does not take into account many factors that arise during epidemics, for example: • Improvement of medical and biological protection methods like new medicine and vaccination; • Changing government decisions affecting the behavior of millions of people; • Economic trends; • Viral mutations.
All these changes require modification of the coefficients of the model and, more importantly, entail the introduction of more processes, such as dynamics of attitudes towards vaccination, social conflicts, and changes in economics.
It may seem rather surprising that the simple SIR SS model, which hardly claims to account for all the epidemiological details of COVID-19, was flexible enough to describe the response of society to virus intervention and subsequent adaptation. At the scale of countries the model was sufficient to describe different kinds of statistics of available data with high accuracy. However, we should note that only available official COVID-19 statistics were taken into account in this research. The SIR SS model is not intended to predict how many times the true fraction of the infected population exceeds the number of total confirmed cases.
Our model allows us to assess the development of the epidemic and strategies to combat it for countries with different social structures and cultural traditions. This possibility could be more important than making predictions. The difference between countries is very clearly represented in kinetic constants.
The model parameters were fitted to various empirical data using a global optimization algorithm (a kind of machine learning). The parameters of social stress response in different countries were estimated. Two of them vary significantly between countries and characterize them: K 2 (exhaustion rate or the inverse time of exhaustion) and q (stress response or mobilization rate). The results ( www.nature.com/scientificreports/ • A very low K 2 value combined with an extremely high q allows countries to avoid the second wave and remain with the classic single peak profile in the dynamics of the infected fraction, I(t). • A very high K 2 will inevitably lead to a double-peaked wave, where the second can be particularly large.
• A low K 2 value, together with a high q and a high a (morbidity), provokes two split waves at I(t) with a protracted plateau at CC(t). The choice of the simplest model of epidemics driven by social stress was almost unambiguous. There are many avenues for future research. Many questions should be answered about the dynamics of immunity, the evolution of viruses, the dynamics of opinion, the logic of government decisions, feedback from social stress to governments, the dynamics of the economy during pandemics, and many others. Each complex model will have many parameters. They need to be determined from the data and the reliability of each value can be questioned. The trade-off between the model complexity and the reliability of parameter determination is a special problem Colombia Iran United States There are also several particular technical problems that are important for COVID-19 epidemics data analysis and modeling. For example, this is a question about hidden cases: a significant proportion of infected people are not counted in the TCC. The proportion of latent cases is not constant and depends on the workload of health services (and its possible overload), the organization of testing, the intensity of epidemics and other problems. To estimate the proportion of latent cases and to correct the TCC, special work on data analysis and modeling is required, in order to fit CC(t) not to the registered, but to the corrected TCC.

Russia
The main result of the research program should be a system of models that will provide us with tools for quantifying various situations and playing out various scenarios in silico to develop anti-epidemic strategies specific to a particular society: country, region or social group. Methods SIR SS model features description. SIR-type models mainly describe the following two types of transitions: S + I → 2I and I → R. The first one defines close contact infection of a susceptible person with an infected person. The second one describes the relaxation process of recovery or death from the disease. The reaction rate obeys the law of mass action, which sociophysics borrowed from chemical kinetics. The autocatalytic form S + I → 2I means that this is a transition S → I with the transition rate proportional to the fraction of infected people I. In state R, a person is no longer contagious and does not return back to state S, therefore, for relatively short periods of time (for example, no more than a year), one can take total population as constant: S(t) + I(t) + R(t) = 1.
However, the basic SIR models do not take into account a person's reactions to the information flow and events taking place in society. In particular, important administrative decisions, which require stricter hygiene standards and social distancing, affect human behavior and disease transmission.
We introduce an enhancement of the SIR model that takes into account the impact of social stress factor 25 . According to the general adaptation syndrome (GAS) theory, there are three chronological stages of the stress response: alarm, resistance, and exhaustion. After some rest the exhausted person can return to the initial state and become susceptible to alarm. Given the fraction of the population not under stress and neglecting the alarm stage, we assume that susceptible individuals (S) can be split in three subgroups by the types of behavior: ignorant or unaware of the epidemic (S ign ), rationally resistant (S res ), and exhausted (S exh ) that do not react on the external stimuli (this is a sort of refractory period). In other words: S(t) = S ign (t) + S res (t) + S exh (t).
In accordance with this add-on, three additional probabilistic transitions caused by social stress factor emerge in the detailed SIR SS model: (1) S ign + 2I → S res + 2I-mobilization reaction (fast). The autocatalytic form here means that the transition rate is proportional to the square of the infected fraction I. This assumption is needed to capture a super-linear increase in alarm reaction to an increase in the number of infected people; (2) S res → S exh -exhaustion process due to COVID-19 restrictions (slow); (3) S exh → S ign -slow relaxation to the initial state (end of the refractory period).
The first transition represents a stress response to a growing proportion of the infected. The second transition describes the exhaustion of resisting behavior. The third transition completes the loop within the S state and reflects the slow leak of the exhausted people to the state of ignoring the epidemic but susceptible to the alarm signals.
The infection rates of S exh and S ign are assumed to be the same as in the basic SIR model. For the fraction S res , the transition to I is much slower, and we take this category as not susceptible to disease.
Possible transitions in the model are illustrated in Fig. 6. www.nature.com/scientificreports/ Differential equations. According to the described transitions, the SIR SS model has five variables: S ign (t), S res (t), S exh (t), I(t), R(t), and can be formulated as follows: where K 1 = q is the rate of transition to the resistant state mediated by a stress response in society, K 2 is the exhaustion rate due to restrictions (K 2 ≪ 1), K 3 is the constant of transition to a state of ignorance (K 3 ≪ 1), K 4 = a is the morbidity rate, K 5 = b is the recovery constant. We kept the following parameters fixed: b = 0.1 (recovery rate I → R, expected period of being infected is b −1 = 10 days), K 3 = 0.01 (S exh → S ign , expected period of being exhausted is K 3 -1 = 100 days). Three other parameters: a (the morbidity rate, S + I → 2I), K 2 (the exhaustion rate, S res → S exh ), and q (the stress response rate, S ign + 2I → S res + 2I) were varied to fit experimental data. An additional parameter that affects the rate of the outbreak and, hence, the time of the first wave maximum, is the initial fraction of infected people in the population I(0) = I 0 ≪ 1. The rest of the population at the initial state is assumed to be ignorant by default: S ign (0) = 1 -I 0 ≈ 1.
The choice of parameters b and K 3 as constants is determined by the following considerations. The dynamics of the outbreak onset (when S res and S exh are small) depends on the a/b ratio. At large time scale, the trajectory of the epidemic's spread depends on the ratio between K 2 and K 3 . Thus, we choose b and K 3 as the references for a and K 2 , respectively. This means that dynamics is much more sensitive to the a/b and K 3 /K 2 ratios than to the absolute values of these constants. If we also vary b and K 3 , we get an unstable (strongly fluctuating in response to a small noise) solution of the optimization problem. Trial simulations showed that the model matches the infection data of COVID-19 for modified values of b and K 3 at the same top level preserving a/b and K 3 /K 2 . This matching confirms our choice.
The transition S exh → S ign closes the loop S ign → S res → S exh → S ign . It describes the slow transition to the initial ignorant state. Qualitatively, this means that humans cannot be exhausted indefinitely. They simply do not respond to external stimuli during a certain refractory period. The transition S exh → S ign significantly affects the dynamics of epidemics. For example, if we take K 3 = 0, then the fraction S res will decay faster, and the second wave of the epidemic will quickly absorb most of the population similarly to the first wave in the SIR model (see the right panel in Fig. 7). This is contrary to the observed COVID-19 data. Figure 7 illustrates the difference between the SIR SS model with and without such a transition.
For this example, we took the SIR SS model fitted to the German data (shown in Fig. 2) and compared it with the same model with K 3 = 0 and other unchanged reaction rate constants.
COVID-19 retrospective data are collected daily. Therefore, all calculations were performed in steps of one day or in steps of an integer number of times less (Δt = 1, 0.5, 0.1, 0.02 etc.). Equations (1)-(5) were transformed into an appropriate discrete-time system.
Comparison with COVID-19 observations. Cumulatively confirmed cases of COVID-19 and daily new cases were obtained from a publicly available source of statistical information: the COVID-19 Data Repository by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University (JHU). The calculation of the dynamics of phase variables according to Eqs. (1)-(5) was carried out using a specially developed MatLab software. The correspondence of the calculated the normalized number of cumulative cases as a function of time, CC(t) = I(t) + R(t), to the normalized observed TCC(t) was estimated using the coefficient of determination R 2 (normalization here means dividing by the total population). As the starting point for the algorithm to work (t = 0), we choose the day since at least 100 cases of COVID-19 were confirmed. The time interval for comparing the two samples was equal to the number of days of the modeled trace of CC(t). The search for the global maximum of the R 2 was implemented in the parameter space (a, K 2 , q, I 0 ), at which the model reproduces epidemics dynamics that are closest to the observed one. The level of R 2 is high: for all countries R 2 > 0.988, for almost all countries R 2 > 0.99, and for seven countries R 2 > 0.998. This means that the time dependence TCC(t) is very well described by the trajectories of the SIR SS model on the selected time intervals. The number of parameters to fit is small (four), therefore there is no sign of overfitting.