Integrating psychosocial variables and societal diversity in epidemic models for predicting COVID-19 transmission dynamics

During the current COVID-19 pandemic, governments must make decisions based on a variety of information including estimations of infection spread, health care capacity, economic and psychosocial considerations. The disparate validity of current short-term forecasts of these factors is a major challenge to governments. By causally linking an established epidemiological spread model with dynamically evolving psychosocial variables, using Bayesian inference we estimate the strength and direction of these interactions for German and Danish data of disease spread, human mobility, and psychosocial factors based on the serial cross-sectional COVID-19 Snapshot Monitoring (COSMO; N = 16,981). We demonstrate that the strength of cumulative influence of psychosocial variables on infection rates is of a similar magnitude as the influence of physical distancing. We further show that the efficacy of political interventions to contain the disease strongly depends on societal diversity, in particular group-specific sensitivity to affective risk perception. As a consequence, the model may assist in quantifying the effect and timing of interventions, forecasting future scenarios, and differentiating the impact on diverse groups as a function of their societal organization. Importantly, the careful handling of societal factors, including support to the more vulnerable groups, adds another direct instrument to the battery of political interventions fighting epidemic spread.


Introduction
The current Coronavirus disease 2019 (COVID- 19) pandemic is more than a health crisis; it is also an economic, social and, in many countries, political crisis. Governmental decisions to fight COVID-19 must consider a range of factors including infection spread, health care capacity, as well as economic and psychological variables, all embedded and linked in a worldwide dynamic context. Importantly, these variables influence each other and evolve over time. The specific trajectory of these interactions will determine the consequences for each nation (e.g., deaths, duration of lock-down, economic loss, etc.). Indeed, for the battle against COVID-19 to be successful, it has been suggested to "inform and qualify action with evidence from behavioral and cultural research" [1] to create high acceptance of and compliance with measures, even after restrictions are lifted. While science offers this evidence to inform policy making [2,3], the influence of psychosocial (PS) variables is not integrated into epidemiological models to predict further disease transmission. Feeling risk of contracting the disease but also other, more remote psychosocial constructs such as trust in the regulating bodies and acceptance of policies [1] may be important to take up protective behavior (e.g., physical distancing, wearing a face mask). In current approaches the time invariance of these behaviors is assumed and adapted when an intervention occurs [4]. Any violation of model assumptions reduces the predictive power of the model and limits its validity to short terms. A potential interaction of psychosocial variables with mobility is such a violation: high psychological strain felt during a lock-down may lead people to reject the preventative measures, which is then related to increased mobility [5]. Consequently, model parameters need to be continuously updated to counter for the accumulating model errors [6,7]. As governmental decisions and interventions need to be planned further ahead, the mechanistic correction of the forecast model by integrating the interactions of PS variables with mobility is crucial. Here we use the weekly cross-sectional COVID-19 Snapshot Monitoring (COSMO) data [5,8], which allow rapid and adaptive monitoring of psychosocial variables and gain a real-time insight into the knowledge, perceptions, and behaviors of the population-the "psychological COVID-19 situation" [9]. Our objective is to improve the accuracy of longer-term predictions of epidemic spread models by integrating the influence of PS variables into disease dynamics.

Integrating psychosocial variables in epidemic transmission model of COVID-19
infectious spread provides a framework that enables a causal analysis and testing of this hypothesis [10]. As a first step, we express the influence of PS variables upon behavior as a causal hypothesis, where its coupling strength and direction, i.e. how variables impact eachother, is inferred from the empirical data. Changes in behavior are expressed by modulation of the contact rate and quantified by the mixing matrix β, which describes the number of infection-producing contacts per unit time. R is the reproduction number of an infection and reflects the expected number of cases directly generated by one case in a homogeneously mixed population where all individuals are susceptible to infection [11]. Both quantities, β and R, cannot currently be modified through vaccination or other changes in population susceptibility, but can be modified by physical distancing and other political interventions. These measures, however, are only effective if large parts of the population comply with them. However, certain groups in the population may have less resources and capacity to adapt and thus perceive different levels of affect, e.g. due to age or variations in personal situations, such as a loss of income during a lock-down, exposed working conditions, or restricted living space. We hypothesize that these groups will show different degrees of acceptance of the measures, leading to variation in mobility and ultimately disease transmission, which then again causally loops onto itself, contributing to further changes in PS levels. However, the mixing between the populations is assumed to be homogeneous as a first implementation of this model, and in lack of any data pointing otherwise.
We formalize this emergent and dynamic process in Fig 1, using a compartmental SIR model for infectious disease transmission comprising three populations of susceptible (S) and infected (I) cases, as well as recovered or dead (R) cases. Note that including Exposed group that corresponds to the case with a significant incubation period during which individuals have been infected but are not yet infectious themselves wouldn't have significantly changed the dynamics from the model [12], and moreover this does not seem to be an important feature of the COVID-19 transmission [13]. The total number of people is N = S+I+R. As a government imposes political interventions J such as a lock-down, the mixing matrix β and the subsequent reproduction number R are reduced. Additional seasonal modulations of β add an Organization of the PS-SIR model. Psychosocial (PS) variables modulate human mobility through changes of the mixing matrix β i . As the perceived consequences of political interventions, such as lock-down effects, accrue, the subsequent accumulation changes the affective state ϕ i of the i-th sub-population. β i affects ϕ i through the influence f i of accumulated PS variables and ϕ i affects β i through the influence g i , e.g., affective risk perception. The index i signals group-specific variation of parameters and variables across the whole population, defining a group structure of vulnerable vs. resilient groups with varying resources and capacity to adapt. Changes in the mixing matrix β i directly affect the dynamics of the infectious disease in the SIR model. https://doi.org/10.1371/journal.pdig.0000098.g001 important temporal variation to the epidemic transmission mechanisms [14]. We expand β along two additional dimensions. First, β can slowly evolve over time as a consequence of modulations with the PS variables, which generally evolve on a time-scale of several weeks to months, as also observed in the rate of the changes observed by COSMO data. The time scale separation between infectious transmission and PS variables justifies the mathematical assumption of weak coupling and its expression as an additive modulation [15]. Second, groups within the population can be differentially affected depending on their resources and capacity to adapt to the situation [1], which we separate into vulnerable and resilient groups for the purposes of the model. The affective state ϕ i of the i-th group is a longer lasting state, which is not caused by a single stimulus but results from an accumulation of experiences. The accumulation of PS effects, u i , accrues over time due to prolonged restrictions, indicating the capacity of the i-th group to adapt to the imposed political interventions. It varies with the pandemic containment efforts such as physical distancing and travel restrictions, but also economic factors such as business closures. Behavioral changes are absorbed in the mixing matrix β i , which establishes the link to the SIR variables and only depends on the political interventions J and the affective state ϕ i of the i-th group.
To infer the strength and direction of the interactions, we use the COSMO data of Germany and Denmark. The measures derived from the COSMO data sets were collected weekly from early-March in Germany (overall N = 11,669) and started a few weeks later in Denmark (overall N = 5,312). The measures focused mainly on the perception of the COVID-19 pandemic, including: AFFECT variables that estimated how serious the threat was to the respondent (e.g., "Is the virus near to you?"; "Is the virus spreading fast?"; "Are you worried about the virus?"), POLICY adjustments and acceptance (e.g., "Should the government close schools?") and TRUST in institutions ("Do you trust the performance of the government?"; "Do you trust the healthcare system?"; see S1 Text in Supplementary Information for more details on questions). Germany and Denmark show both an early increase in the concern around the spread of the virus (AFFECT), and agreement with the governmental interventions (POLICY), especially around school closings. This peaked around March 23, when the country-wide lock-down was instituted in Germany (Fig 2). This trend declined across time, indicating less concern with the pandemic, to levels similar to where they stood when data were first collected. The data from Denmark were collected starting immediately after the country-wide lock-down but showed a comparable declining trend with more recent measurement waves. Fig 2 shows the time-aligned estimates [4] for COVID-19 infection rate and basic reproduction number over the same time period, which were used as a ground truth for comparison with the forward simulations and sampling performed with our model. Fig 2 also shows the mobility of the population based on cell phone tracking indicating sharply reduced mobility during the lockdown and then increasing mobility as AFFECT and POLICY declined [16].

Sampling of psychosocial and behavioral interactions during COVID-19
In order to assess the strength, direction, and significance of the coupling of PS variables on the time-varying reproduction number R, we submit the above extended SIR model with PS interactions to Markov Chain Monte Carlo (MCMC) Bayesian inference. Three parameters are inferred: γ and c are the strength and direction of the couplings g i and f i , respectively (Fig  1), and Γ is the time scale of PS accumulation. For simplicity we dropped the group-specific indices of γ and c, as we do not have access to group-specific mobility data from Google. Fig 3  shows the results of the inference and the data used in the process. The two top rows show the time series for the reproduction number R(t) estimates [4] for Germany (DE) and Denmark (DK), respectively, which we are using as ground truth in this workflow, and the corresponding fit of the model as lines with the same color scheme. The bottom row shows posterior distributions for model parameters of interest along with the priors (in grey) placed on those parameters. While in Bayesian terms, the coupling strength γ is well identified by the data of both Denmark and Germany, Γ, c are identified strongly only for data from Germany. For Denmark these posteriors follow the prior, which are likely due to the lack of COSMO data during the onset of confinement in Denmark. As our primary goal is to determine whether γ is non-zero, i.e., PS factors have an effect on reproduction rate, we implement our null hypothesis as a normal prior on γ, centered at 0.1, i.e., a small effect. MCMC sampling from the posterior distribution allow the data to then identify γ, Γ, c. The point estimates (medians) of γ are 0.21 and 0.23 for Germany and Denmark, respectively, and 95% confidence In the first two rows, estimates for daily infections [4] and reproduction number R are presented as time varying means and standard deviations. Third row shows three selected variables from the AFFECT group in the COSMO surveys (see S1 Text in Supplementary Information for more details; higher values indicate that COVID-19 was seen as something one does not think about all the time, rather not terrifying, and something not to worry about); their weighted mean (derived from PCA) labeled as COSMO(t) is used in the model inversion. The last row presents the scatter plots of the Google mobility report overlaid with lines marking the 7-day weighted moving averages, where the mean of the blue "retail and recreation" variable is used in the model inversion.
https://doi.org/10.1371/journal.pdig.0000098.g002 intervals are (0.17, 0.25) and (0.18, 0.30) for Germany and Denmark, respectively. This indicates that PS variables significantly impact the reproduction rate with a non-zero and positive coupling strength, i.e. the less stress people felt about COVID-19, the higher R. This defines a baseline, which allows individual parametrization of the models for Germany and Denmark.
Note that the surrogate model corresponding to a null hypothesis where γ = 0 is available to the inference process, because the prior on γ covers the value 0 with higher probability than the posterior values that we report. This indicates that the posterior has been well informed by the data. We could have evaluated an equivalent model with γ fixed to 0, and performed a model comparison, but this would be identical to allow the model to explore the parameter space around γ = 0 which is what happens with MCMC with the prior we have chosen.

Exploring the effects of PS variation on pandemic spread
We initialize each model with its respective past values for epidemic, mobility, and PS variables and provide simulations for the upcoming months. To estimate robustness of the forecast coupling strength of f i ) for Germany and Denmark (blue and orange histograms, respectively) and their respective priors (grey lines). On the left, the histograms of posterior samples for γ for both countries, which deviate significantly from the prior (grey line); the data thus identify a strong and robust effect of PS factors on SIR dynamics. Similar effects are seen for Γ and c for Germany, while Denmark follows the prior on these parameters likely because of missing COSMO data at the onset of confinement.
https://doi.org/10.1371/journal.pdig.0000098.g003 simulations and sensitivity upon the inferred parameters, we run simulations for all values of coupling parameters sampled from the posterior distributions. The panels of Fig 4 show a comparison of simulations for the scenarios with and without PS coupling, displaying the forecast trajectory evolution until April 2021 for the fraction of susceptible S and daily new infections derived from the infected individuals I, reproduction number R, affective state ϕ i , and estimated strength of political interventions J, which for March of 2020 is set so that the dynamics after the intervention would follow the observed infections and R. For the subsequent governments' relaxations and lockdowns, J is set to follow the relative strictness between the measures at different periods.
The left column of Fig 4 displays the situation for Germany, the right for Denmark. In addition to the good agreement with the estimates for the pandemic spread until the summer 2020 [4], the model is also capturing the subsequent estimated and reported peaks and declines in the infection rates of both countries all the way to the spring of 2021 [17,18]. Further differentiation comprises group structure as indexed by i, simulating the population into groups of vulnerable and resilient based on different levels of PS variables in response to political interventions. The group differentiation is realized by assigning different values to coupling parameters c i and Γ i .
The seasonality of transmission [14] is included in the simulations where more transmissions are assumed in autumn/winter than in spring/summer. Together, PS variables and seasonality cause the predicted epidemic wave in autumn/winter of 2020 to be larger than the previous one in spring 2020, despite the same stringency level of lock-down. This effect is also partly responsible (together with the government interventions) for the lower level of infections observed during May and June 2020 than later on in autumn. For every simulation of the SIR model without PS variables (magenta dashed lines in Fig 4A), the efficacy of political interventions J is overestimated, showing an order of magnitude smaller second wave, than the one that was observed and predicted with the full model. When comparing the forecast epidemic spread with and without PS coupling (γ = 0), the political intervention J would have proven sufficient to control the epidemic spread in the latter case, but not the former, because the PS variables cause a surge of infections towards the second part of 2020. This is even better illustrated when comparing the observed [17,18] curves for S and new daily infections in S5 Fig, which contain the same peaks and values within the same ranges as the model, up to the May 2021, i.e. S between 85-90%. The model without the PS subsystem misses the subsequent peaks, and has an order of magnitude smaller ratio of susceptible (S<2%). Even though the mean inferred impact of the PS subsystem is slightly lower for Germany than for Denmark, values for R and the second infection wave related to the seasonality is predicted to be larger in Germany.
The sensitivity of the simulations to the choice of parameters is explored in Fig 4B. Here the coupling parameters γ are drawn from the posterior distributions for a sample of 4,000 model realizations for each country, and then used in the simulations with all other parameter settings equal. The shaded areas indicate a confidence interval of 95%. All simulated trajectories show qualitatively the same behavior, that is a surge of infection towards the end of 2020. To provide a quantitative impression of the strength of PS coupling as compared to lock-down and seasonal variation, Fig 4C shows Supplementary Information, S2 Fig): Red is the impact of the government interventions (the largest on 19/03, and the subsequent relaxations on 13/05 and 30/06); blue is seasonality (the maximum impact, and the impact for 2 months after the relaxation on 10/05 to 10/07). Green is the impact of the PS subsystem, where light green is the cumulative impact before the relaxation on 14/04 and before the predicted intervention in autumn; and dark green is the predicted impact in autumn 2020, separated per group.
https://doi.org/10.1371/journal.pdig.0000098.g004 effects, the effect strength has been inferred for a particular time after the onset of the cause, in particular the cumulative impact of PS is shown before the relaxation on 14/04 and before the predicted intervention in autumn (light green); furthermore, the cumulative impact of group structure is also shown (dark green), clearly demonstrating its epidemic spread-reinforcing effect. Political interventions via mobility restrictions show the strongest effect, seasonality is factor 2 smaller than mobility constraint, and a factor 3 to 4 larger than PS factors.

COVID-PS Model
We express the influence of ϕ i and β i of the i-th group by the coupling functions f i (u i , J, ϕ i ) and g i (ϕ i , J, β i ), which shall be inferred from the data. As a first step in model building, we recognize similarities with many physical and biological systems, which are coupled indirectly through a slowly evolving dynamic medium. Such slow modulations of a coupling medium are effectively equivalent to a linear weak coupling of the dynamic systems, prescribing their interaction around a working point as additive in a first approximation, but leaving the mathematical form of the actual interaction functions themselves open. We harness this effect to extend the transmission dynamics of the SIR model to include psychosocial interactions and infer the precise form of the influence functions using Monte Carlo Markov Chain (MCMC) methods. The full set of equations reads: where the coupling functions g i (ϕ i , J, β i ) and f i (J, β i , t) affect the growth rates of the mixing matrix β i , effective state ϕ i , while the assumed infectious time is set at 14 days, i.e. α = 1/14 [4], instead of being modelled as a distribution [12]. The time is hence implicit and set in days. The latent period of the disease, even though important for the evolution on shorter time-scales [19], is ignored in our model due to the focus on the long term evolution. For similar reasons and to focus on the impact of the PS variables, the model also does not contain removed population, and S and I are relative to the total population, so that S+I = 1. Since we are interested in the dynamics of the pandemic variables, all the other units are assumed phenomenological and set accordingly to this. The death rate is set to be equal for the both groups, μ = 1/365/80 and the total birth rate of both groups the same, with each group having Λ i = μ/2. These rates are at much slower timescales than the relevant ones for the diseases, but we have nevertheless introduced them to allow for a consistent division on groups, which is performed through Λ i . For simplicity, and due to lack of other priors, the mixing matrix is assumed to be identical among the groups.
The back-coupling of β i on ϕ i is communicated via an accumulation of lock-down effects, aka deviations from normal mobility, through an auxiliary variable u i . For both types of interactions, the political intervention J is included in the coupling and subject to the inference.
The coupling function f i is assumed to increase linearly and then saturates at a maximal value of 1, whereas g i shows the same behavior but remains restricted to the linear signal range. The strength and direction of the coupling is equivalent to the scale and sign of γ i and c i of the couplings g i and f i , respectively, subject to inference from empirical data using a stochastic model (see below). As we will not discuss group-specific mobility changes due to changes in PS variables, we drop the index of γ i , leading to g i ð� i ; J; b i Þ ¼ Jb i þ g i � ð� i À � 0 Þ and f i ðJ; b i ; tÞ ¼ hJ þ c i tanhð2u i Þ, so that the following closed loop equations appear: Here u is an auxiliary variable for the accumulation of restrictions-related fatigue. As a such, u i = 0 for J = 0, i.e. when there are no restrictions, and during the relaxation of restric- During restrictions, it linearly increases with their severity, which is accounted through the impact on the mixing matrix, i.e. larger the deviation of β i from β 0 , larger is the fatigue. In the same time, each group can be differently affected by the interventions, and this is scaled with the parameter Γ i . However, the way this affects the affective state (or the happiness) ϕ, is set to be a sigmoid, saturating at levels that are driven by c i , which is inferred from the data together with Γ i and γ i . Here Γ i /σ <<1 reflects the time scale separation, that is both mobility via β i and affective state via ϕ i evolve on a time scale of the order of 1, whereas PS effects accumulate on a time scale of the order of Γ i , hence we set σ = 30. The impact of the government interventions J on ϕ is smaller than on β, and hence h = 0.2, meaning the mixing matrix is stronger influenced by the interventions, than is the affective state ϕ. The overall coupling of the SIR and PS sub-systems is driven by the parameter γ i , which describes how much PS variables impact the dynamics of the epidemic. For the asymmetry between the group parameters c i and Γ i , whose mean values were inferred, we set the first group to have 9 times larger values than the second, corresponding to 80% larger or smaller values than the mean for each respective group.
It should be noted that depending on the nature of the phenomena, the coupling between them is chosen to be sigmoid, if the impact is supposed to saturate, or linear otherwise. These are the simplest choices due to the phenomenological nature of the model, even though the coupling functions [20] and their proper recovery [21] are of crucial importance for biological oscillations [22].
To describe the logic behind the equations, absence of interventions or an impact of the PS variables, the first term in _ b leads to a relaxation of the transmission rate towards pre pandemic levels over time. The second term lowers this asymptotic value due to government interventions, and the last term increases it due to the decrease in the affective state ϕ i , as compared to its base value ϕ 0 = 1. The strength of this effect is driven by the coupling variable γ i that is being inferred. The affective state ϕ i is similarly relaxing to its base value ϕ 0 in absence of interventions J, which make it to decrease exponentially to a new level dependent on the strength of J. Similar impact on ϕ i is observed by the accumulated fatigue, which is a sigmoid with a time constant of months, and with a level set by the inferred parameters c i .
To include the impact of seasonality on the transmission ability of the virus [14], we adapt sinusoidal modulation of β 0 with a period of one year and the peak being in Mid-January for the northern hemisphere. The seasonally varying β 0 (t) hence reads: with J s being the magnitude of the variation that depends on the latitude and the climate. For both Germany and Denmark, it is fixed at 40%, which is the level predicted for New York [14].
To calculate the relative impact of the different parts of the model (PS part, seasonality, and the government interventions), for each country we run forward simulations that exclude certain parts of the model and we calculate the difference in the observed numbers for given periods of time.
The model is numerically integrated in Matlab using the Heun integration method with a time step of 0.05 days. A moving average of 1 day was then used to bring the results of the model to the same resolution as the empirical data. R(t) is calculated as R t ð Þ ¼ L i b i mðaþmÞ , and for better comparison between groups and models it is always normalized as R of different models/groups were referring to a same size of population (effectively the values in the case of groups are multiplied by 2).
The statistical model was implemented in v2.22.1 of the Stan probabilistic programming language [23], with parameters sampled using the No-U-Turn variant [24] of Hamiltonian Monte Carlo [25]. The posterior draws were verified with the Stan 'diagnose'utility, which checks for transitions that hit the maximum tree depth, divergent transitions, low estimated Bayesian fraction of missing information (E-BFMI), low effective sample sizes andR values below 1.1 (indicating convergence of the MCMC method). Posterior distributions for key parameters and posterior predictive distributions for data were constructed, verifying that the data lie in a 95% confidence interval, to ensure model fit is reasonable.

Grouping PS variables in the COSMO data
Given the importance of the group structure on the relation between PS variables and SIR suggested by the model, we further analyzed the COSMO survey data to determine if significant grouping variables could be extracted. We tested age (young, middle, and old age), sex, and education (three levels) and compared the aggregate PS variables across all available waves using multivariate Partial Least Squares (PLS) analysis [26,27]. Employment status and household income were less consistently collected so were not considered for this analysis. In both German and Danish data, differences based on sex ( Fig 5) and age (Supplementary Information S1 Fig) were observed. Differences related to education were less reliable but existent. We also observed that the aggregate PS measure was reliably related to subjective psychological stress and resilience in a subset of the COSMO data, suggesting the aggregate captured broader psychosocial aspects beyond risk perception (see Supplementary Information S4 Fig). For Germany, data collection started earlier and shows an initial overlapping transient with a faster rebound for males than for females (Fig 5 left), resulting in a group difference which remains preserved as time evolves. The same situation is observed for Denmark, even more pronounced (Fig 5 right), but without the initial transient as data collection started later. The interaction of age and sex was not significant, indicating these were additive factors in these data. Older adults and females showed larger overall values on the PS effect, suggesting these groups felt higher affective risk and that the COVID-19 pandemic was more of a crisis relative to younger groups or males.

Discussion
Although mechanisms of COVID-19 transmission are indiscriminate, this may not extend to a complex society fighting epidemic spread [28]. Societies have a rich inter-individual variability, different cultures, wealth distribution, and ways of functioning, but most epidemic models  Table A in S1 Text in Supplementary Information for a complete list of variables. Error bars: Bootstrap estimated 95% CI. Bootstrap ratios are the singular value weights divided by their standard error and are roughly equivalent to a z-score. The shape of the Psychosocial score curve relative to the PS Indicators suggests an increase in affect risk perception and great institutional and policy trust in March, followed by a general relative reduction.
https://doi.org/10.1371/journal.pdig.0000098.g005 describe the dynamics of a homogeneous population in different disease stages. The model parameters absorb all properties characterizing the diversity of society including psychosocial factors. With the introduction of psychosocial coupling in epidemic disease models, the consideration of effects linked to societal diversity becomes possible.
We estimated the strength and direction of psychosocial causes within a causal inference framework and demonstrated their effect upon the infectious spread to be smaller, but on a similar scale of magnitude as current political interventions (approx. factor 8 smaller than lock-down effects and social distancing) and seasonal variation (approx. factor 4-5 smaller). In particular, affect variables that estimated how serious the threat was to the respondent provide strong evidence for the causal effect of psychosocial factors on mobility and subsequently infectious spread. Our cause-effect analysis considered minor parametric changes due to psychosocial variation, but has the capacity to integrate a large range of multi-dimensional data including societal, national and cultural differences, political strategies, and trust in the government and media. For instance, multidimensional approaches could integrate additional psychological and sociocultural indicators, affecting the mixing matrix, thereby establishing a battery of mutually dependent behaviors. Given the here estimated size of psychosocial causes, the synergy of these interactions may likely reach scales that cannot be ignored in epidemic spread. Furthermore, heterogeneity within a society generates groups that are not equally affected by the same political intervention. The group structuring differences are plotted in Fig  5 and supplemental figures, demonstrating that a highly heterogeneous society would show acceleration and augmentation of infections across the entire population due to group structure (S1 Fig). Despite the fact that the resilient population experiences no change in its affective state, the uniform mixing of the populations causes increases of total numbers of infections in both groups, vulnerable and resilient. A detailed quantitative understanding of these effects provides an added means to guide political decisions to influence the infectious spread by considering parameters linked to demographic factors of societal organization. In the continuum spanned by group heterogeneity and timing of interventions, the choice of the best lock-down duration will depend on the heterogeneity within the social structure and allow the optimization of political interventions for a given society, and thus minimize economic damage while respecting the efficiency for epidemic transmission control. This is particularly important as interventions have optimal working domains (such as lock-down strength or duration), but outside of these windows they are inefficient. For example, in S1 Fig A for Denmark, we show that extending lock-down durations do not alter the trajectory of epidemic spread significantly. Alternatively, the same mechanism in our model may be used to increase the efficacy of political interventions by changing societal factors, in particular reducing societal heterogeneity and increasing support to the more vulnerable groups [1].
Partly because of a lack of specific data in the COSMO surveys, we have left explicit economic factors aside in this discussion, although such factors would be relevant to the grouping structure. It has been suggested that psych-socio-economic variables such as income, education, and ethnicity play an important role in the transmission of infectious diseases [29][30][31], not only because they reflect resilience, but also because they influence the adoption of protective behaviors during pandemics [32]. As a consequence, models taking into consideration PS-SIR interactions can provide stronger support for political interventions that support less resilient groups socially and economically as this is a factor in reducing heterogeneity [33,34]. The model framework also underscores that psychological knowledge from personality and health psychology can be put to very concrete use in a crisis [35]. Future efforts should quantitatively assess the impact of communication and behavioural interventions that address PS factors to enhance compliance.
An important caveat in our approach is the uncertainty of the ground truth, against which we are validating our results. Underreporting of the cases has varied greatly during the course of the pandemics [36], but still several reliable datasets appeared trying to infer and fit the pandemics evolution as reliably as possible [37]. Using this as ground truth, we demonstrate that the model with PS variables can qualitatively capture the evolutions of infections, which is not the case when these are omitted.
Similarly, besides the mobility which is the focus here, reactive behavioural changes, in particular, how we come into contact with one another, are also impacting the mixing matrix [38]. It would be important for these to be also accounted for in future models that account the PS factors. In the same direction, having a true mixing matrix β ij with group-specific transmission rate modulators, instead of the assumption about equal risk of infections per group, would also make the model closer to reality. Accounting for the latent period of the disease is also important for the future modifications, especially since it is strongly influencing the growth rate of the incidence [19]. In the same line, the impact of the seasonality in the spread of COVID-19 pandemic seems to be well established [39,40], and in further developments of the model, its impact would need to be also inferred from the multiseasonal data. Finally, COVID-19 pandemic is continuously evolving, with new variants appearing with different virulence, but also different vaccination strategies across countries [41], which both would need to be included in future models.
The COVID-19 pandemic is exacerbating existing global and national inequalities, which typically hit the most vulnerable groups hardest [42,43]. The focus of policymakers has been on containing the spread of COVID-19 and mitigating the socioeconomic effects of the pandemic. Our findings demonstrate that addressing these vulnerabilities nation-or even worldwide through policy action (for instance, temporary basic income strategy [28,44,45] may not only be a mean to fight the socioeconomic consequences of COVID-19, but also represent a mean to directly fight its spread across all groups, both vulnerable and resilient, through nonpharmaceutical actions containing the disease. Supporting information S1 Text. Table A Table A in S1 Text for a complete list of variables. Bootstrap ratios are singular value weights divided by their standard error and are roughly equivalent to a z-score. (TIFF)

S4 Fig. Age differences in PS indicators.
The dominant latent variable (both p<0.001 by permutation test) of a multivariate PLS analysis of mean change over time of PS indicators in German (A-B) and Danish (C-D) COSMO data by age group. Subjects were grouped into three equally sized bins of age categories. In the German data, age categories were 18-37, 38-55 and 56-87 years. In Danish data, age categories were 18-49, 50-63 and 64-92 years. AFFECT (AFF), POLICY (POL) and TRUST variables were included if data were available from all waves in each of the COSMO datasets. See Table A in S1 Text for a complete list of variables. Error bars: 95% CI. (TIFF) Fousek.