Complexity of Coupled Behaviour-Disease Models and Their Relative Performance Against Empirical Data

Anonymous

Abstract

Non-pharmaceutical interventions in the early stages of the COVID-19 pandemic reduced case incidence, but cases later resurged once restrictions were relaxed, contributing to additional pandemic waves. The temporal unfolding of behavioural changes in populations presents a challenge to mathematical models for disease dynamics. Coupled behaviour-disease models with varying complexity accounting for diverse factors have been used to capture behavioural dynamics and S.A.R.S-CoV-2 transmission, yielding different results. We study the impact of model complexity on their predictive power by formulating five coupled behaviour-disease models with varying structure and number of parameters. We fit these models to case notifications and stringency of control interventions from five European countries in the first Spring 2020 wave, and study how well these fitted models predict the second wave. Models with more parameters do not necessarily have a greater ability to explain and predict key features of a pandemic wave. Thus, a relatively simple coupled behaviour-disease model with important parameters can provide adequate information about the pandemic wave. Additionally, complex models can be country-specific, working better for some countries and poorly for others. We conclude that modellers should not always use the most complicated available models if the data do not support their use.

Plain Language Summary

The spread of any infectious disease largely depends on the population's reaction. The population either facilitates a rapid spread or curtails it through control strategies. This phenomenon reveals the relationship between an infectious disease and the population and how they influence one another.

The population's goal is to eradicate the disease or at least mitigate it; adopting control measures will imply a change in population behaviour. Mathematical modellers have formulated various models to study how the disease behaves and/or how the control measures change the disease behaviour. The evolving nature of the behaviour of individuals creates a modelling challenge, particularly one that can be tested against data.

We consider five different ways to build a behavioural model, and then we add it to a standard disease model for the COVID-19 virus. We estimate the parameters of the models using daily reported case data from five European countries during the Spring period of 2020 and predict the disease's behaviour for Fall 2021. Our findings show that models with more parameters did not necessarily help to understand the disease better or predict significant features of the disease behaviour. But a relatively simple model that accounts for key parameters is sufficient to explain the disease behaviour and produce accurate predictions. Thus, mathematical models are useful and can play a significant role in decision-making when formulation is based on data.

Introduction

Mathematical models have become useful tools for understanding the dynamics of infectious disease spread within a population. These models have benefited from several contributions since the foundational compartmental model was formulated by Kermack and McKendrick from 1927 to 1933. These contributions have resulted in a better understanding of disease dynamics in the context of data availability and improved screening processes. Likewise, social dynamics, environmental changes, climate conditions, migration, urbanisation, and many other factors have now been identified as elements that"manipulate" the dynamics of disease outbreaks. Although these factors have added to the complexity of epidemiological models, they have greatly contributed to new discoveries, better decision-making processes, and the implementation of mitigation strategies, resulting in reduced infection and death rates, as well as the formulation of public health policies. Aside from the benefits to epidemiological studies, mathematical models have facilitated the design and implementation of experimental scientific studies involving either living organisms or chemicals, as well as curated data collection involving human subjects.

One form of adding complexity to models has been through the channel of social dynamics, whose impact on disease dynamics cannot be overlooked. It has become evident that the dynamics of a disease influence social dynamics, which, in turn, also affect the dynamics of the disease – illustrating a feedback mechanism. Such models can be used to study the influence of non-pharmaceutical interventions such as social distancing, mask-wearing, and international and local border controls, as well as pharmaceutical interventions such as vaccination and antiviral drug administration, factors like fear, opinions, rumours, or information spread, and the effects of social norms. These are characterised by a feedback loop between the disease system and the intervention strategy (which affects social dynamics). The COVID-19 pandemic illustrated a strong feedback loop between social behaviour and the disease. Many mathematical models have attempted to capture this interplay in various forms, largely by considering different approaches to modelling the behavioural responses of the population.

Models of the evolution of social dynamics formulation have been based on at least two main factors in a broad context: belief-based and prevalence-based, while others include incidence-based (new cases instead of prevalence), and the payoff (utility) individuals derive from their decisions. Prevalence-based behaviour changes take different forms based on current prevalence, the current and past history of disease prevalence, and media reported prevalence instead of actual. Utility-based behaviour can be estimated as a function of prevalence, as a fixed rate independent of current prevalence, or through social learning (imitation) influenced by the structure of society. The formulation of these models, coupled with disease models, results in coupled behaviour-disease models that take different structural forms of the interaction between the social system and the disease system. One such structure is interacting networks, which consist of at least two layer networks where one network layer represents disease transmission, and the second network or other network layers could represent social communication for the influence of information or behaviour in a scale-free or random nature.

A simple yet useful way to formulate a coupled behaviour-disease model is to use the concept of compartmentalisation of the population. Most network models make use of the foundational concept of this approach. In this modelling technique, the complete system is made up of two subsystems: a disease subsystem and a behaviour dynamics subsystem. The disease (epidemiological) subsystem divides the population into groups based on their status: susceptible, exposed, infectious, and recovered, where transitions from one state to another are described by ordinary differential equations. The behaviour dynamics subsystem can take two forms: 1) a further subdivision of disease groups based on their beliefs, risk perception, or adherence to interventions such as mitigating susceptible and non mitigating susceptible groups; 2) a game theoretical model of behaviour response when individuals engage in a game of decision-making. The two subsystems are coupled mostly through the transmission (contact) rate parameter (modelled as a function of behavioural factors such as the number of individuals engaging in mitigating strategies to decrease disease prevalence, increase the rate of individual awareness, or encourage the adoption of protective interventions with respect to prevalence level or information spread), as is the case with S.A.R.S-CoV-2 transmission. Behavioural models that are modelled using game theoretical frameworks, where individuals are primarily interested in maximising their utility, operate on the notion that individuals make decisions based on several factors and can change their decisions over time as the epidemic unfolds.

One challenging question most researchers face at this stage is how many factors to consider in the model formulation and which ones hold significant importance in arriving at the research objective. More complex models possess the ability to capture intricate patterns in the systems; however, the purpose, availability of data, and contextual structure guide the complexity of models. While some consider a simple behaviour-dependent transmission rate in the context of precautionary measures such as N.P.I's, others look at simple utility functions, such as the costs of mitigating and not mitigating; some extend this further to include social norms, while others incorporate the homophily of social groups, among many other factors. This mostly involves at least two strategies, where one strategy (mitigating) causes a change in disease prevalence and the other maintains a"status-quo" strategy (non-mitigating behaviour). The number of strategies (behaviour types/players) can exceed 2 to N finite numbers of behaviours, with utility functions assigned to each strategy based on perceived or expected payoffs (cost or benefit). These utility functions can take very simple to complex forms. All these considerations contribute to the complexity of the final model, but in the context of the model's ability to fit empirical data and perform predictions, how important are these?

This brings us to the focus of this study: how complex should models be for the purpose of fitting and prediction? Very few coupled behaviour-disease models are validated against data due to a lack of data on population behaviour. However, the S.A.R.S-CoV-2 virus has led to various data collection efforts on the behavioural responses of individuals in many countries.

We take advantage of this data to evaluate the predictive powers of five coupled behaviour-disease models with slightly varying behaviour models and different numbers of model parameters. Here, we investigate the potential of these different models to predict the second wave of case notifications during the COVID-19 pandemic, having fitted them to the first wave. We investigate which parameter considerations are significant for this purpose; that is, how the parameters influence the model prediction by working with data from five European nations (Austria, the Netherlands, Portugal, Denmark, and the United Kingdom).

Results

Modelling Framework

We consider five coupled behaviour-disease models for S.A.R.S-CoV-2 virus transmission in this study, where the structure of coupling the disease model and the behaviour model is the same for four of the models. The disease model is a seasonal Sir compartmental model with three groups: susceptible (S), infected (I), and recovered (R). The susceptible group consists of individuals who have a high probability of infection; the infected group consists of individuals who are hosts of the virus and can transmit it to others; and the recovered group consists of individuals whose health has been restored to its previous state or who are dead as a result of the infection. The disease model remains the same for all five coupled behaviour-disease models.

The differences in these models arise from the structure of the behaviour dynamics models, which are formulated using the concept of evolutionary game theory. We consider two behaviour types: mitigator (M) and non-mitigator (N) — with respect to N.P.I's within the population in the presence of S.A.R.S-CoV-2 virus transmission. Mitigators constitute a section of the population that supports the use of N.P.I's for controlling the spread of the virus, while non-mitigators are individuals who do not support N.P.I's. Each of these behaviour types is assigned an expression for perceived payoff, which forms the behaviour model, similar to the replicator equation.

Full details of each model can be found in the Methods section. The disease transmission and baseline behaviour-disease coupling framework adapts the model introduced in Ref., which we extend here to compare multiple behavioural model structures.

We fitted the models to empirical data from 5 European countries (Austria, the Netherlands, Portugal, Denmark, and the United Kingdom) to estimate the model parameters and initial conditions generating the posterior distributions for each model and country. The daily incidence from the model is fitted to the case notifications reported by each country, and the behaviour is fitted to the Oxford Stringency index. The total length of data considered for the study is identified as having two major waves: the first wave (March — May) and the second wave (September — December). The models are fitted to the first portion of the empirical data, which contains the first wave (see Table 1) using the Approximated Bayesian Computation Sequential Monte Carlo (see Methods for details).

: Table 1 summary: Data lengths for case studies in five European countries, all ending on January 28, 2021. The United Kingdom has the longest total data length at 364, followed by Austria at 339, and Denmark, the Netherlands, and Portugal with lengths ranging from 333 to 337. Cutoff dates for fitted data points vary by country, with the United Kingdom having the most at 250, followed by Denmark at 230, while Austria, the Netherlands, and Portugal each have 200.

First Wave Fit

The five coupled behaviour-disease models fit well to the first wave(March — May) of case notifications for Austria, the Netherlands, and Denmark (Figure 1 A, B, D). The models 1, 2, and 3 produce a better fit for Portugal compared to models 4 and 5 (Figure 1 C). We observe that models 1 and 2 fit the United Kingdom, whereas the other models fit poorly (Figure 1 E) based on a time series analysis.

Figure 1 summary: Five line charts showing disease incidence over time for Austria, Netherlands, Portugal, Denmark, and the United Kingdom. Each chart compares five predictive models (represented by colored lines and shaded 95% credible intervals) against empirical data (black line), with a vertical dashed gray line separating training data from testing data. Across the countries, models show significant variance in predicting the timing and magnitude of the second wave, with most failing to tightly align with the actual empirical data. The takeaway is that the models demonstrate varying and often inconsistent predictive performance for incidence trajectories following the initial training wave.

We consider an analysis of the models capability to fit the first wave peak such that the peak magnitude and day are closely matched by our model estimates. We observe from this analysis that the estimated peak magnitude for Austria was lower than the actual for all five models, while the model estimates closely match those for Denmark and the Netherlands. Models 2, 3, and 5 are identified as being below the actual value for Portugal, whereas Model 1 provides a better estimate. We observe in a similar manner that Models 1 and 2 provide good estimates for the United Kingdom compared to Models 3, 4, and 5, which have much higher peak values than the actual(Figure 2 A). A corresponding analysis of the peak day using a seven-day interval criterion from the actual peak period reveals that Models 1, 3, and 4 are accurate estimates for Austria. Models 2, 3, and 4 are accurate estimates for Denmark; Models 2 and 3 are accurate estimates for the Netherlands, whereas Model 2 is the only accurate estimate for Portugal. Our results reveal that no model accurately estimates the actual peak period for the United Kingdom with respect to our criterion. Nonetheless, when we consider the 95% credible interval, we observe that some of the models meet the criteria for all the countries (Figure 2 B).

Figure 2 summary: Two sets of plots comparing five models against observed data for five countries. The top set shows peak magnitude on a log scale, where most models align closely with observations across Austria, Denmark, the Netherlands, and Portugal, but show significant variation and overestimation in the United Kingdom. The bottom set shows the peak day, where model estimates are generally consistent with observations in most countries, except in the United Kingdom where Model 3, 4, and 5 significantly overestimate the peak day. The result is that while model fit is relatively consistent across most nations, there is substantial divergence when estimating the peak timing and magnitude for the United Kingdom.

Statistically, we employ the Adjusted Akaike Information Criterion (A.I.C.c) for evaluating model performance and comparison in two ways: incidence ( I i ) only and incidence plus behaviour ( I i + X ). We observe that Model 1 outperforms all other models for all five countries and has results similar to Model 3 for Denmark and the Netherlands using the incidence only A.I.C. Model 1 still outperforms all other models for all countries using incidence plus behaviour. In both approaches to calculating the A.I.C.c, our results show that Model 1 demonstrates the best goodness-of-fit in comparison to the other models for all five countries (Figure 3). Model 3 yields the second best A.I.C value on average, based on similar observations. Thus, considering all analyses of estimates from fitting, we can conclude that Model 1 fits the best compared to the other four models.

: Figure 3 summary: A grouped bar chart showing the AICc values across five countries for five different models, each compared in two versions: incidence only (I) and incidence plus behavior (I + x). Across all countries, the models using both incidence and behavior generally show lower, more negative AICc values than those using incidence alone. The point is that incorporating behavior into the model improved the fit for all models in every country tested.

Second Wave Prediction

We evaluate the predictive capabilities of the coupled behaviour-disease models using the time series trajectories, peak analysis, and statistical metrics. We observe from the time series trajectory that Model 1 predicts the second wave of case notifications in all five countries (Figure 1). Model 3 predicts that, for Austria and the Netherlands (Figure 1 A,B), Model 2 provides a good prediction for Austria and Denmark (Figure 1 A,D). All other results show an early and high amplitude second wave prediction. The worst performing model in terms of time series observation is Model 4 for all countries except Portugal (Figure 1 C). despite its ability to fit the first wave for those countries (Figure 1 A,B,D,E). Although the average simulation of the time series performs poorly, their corresponding 95% credible interval captures the actual case notifications for most of the models.

Prediction Peak Magnitude Analysis

The peak analysis provides useful information to assess the performance of the models in terms of prediction. We observe that the 95% credible interval of the predicted peak magnitude for Models 1 to 4 contains the actual peak magnitude for Austria, while Model 5 estimates a high range relative to the actual (Figure 4 A). Models 1, 3, and 4 show that their 95% credible interval contains the actual peak magnitude for Denmark, while Models 2 and 5 predict a high range (Figure 4 A). For the Netherlands, we identify that Model 3's prediction range was lower than the actual value, while the other models do contain the actual value.

Figure 4 summary: Two sets of multi-panel charts showing the estimated peak magnitude (top, log scale) and peak day (bottom) of a second wave for five countries. In both charts, five different models are represented by distinct colors. For peak magnitude, predictions vary significantly across countries and models, with some models showing much higher peaks and wider uncertainty intervals than others. For peak day, predictions generally cluster between day 200 and 350, though the degree of uncertainty varies by model and country. Overall, the figure demonstrates that model predictions for both the timing and scale of the peak are highly variable across the tested models and countries.

The similar outcome is observed for Portugal, but in this case, Model 4 predicts a much lower range relative to the actual peak magnitude. Lastly, for the United Kingdom (U.K), we identify that Models 1, 2, and 4 predict a 95% credible interval that contains the actual peak magnitude, while Models 3 and 5 predict a lower range. Interestingly, we observe that the average predicted peak magnitude values for Austria, Denmark, and Portugal lie in close proximity to the actual values for Models 1 and 3, although Model 1 appears to have a very good approximation for all five countries (Figure 4 A).

Prediction Peak Day Analysis

We assess the period when the peak magnitude occurred, first using the 95% credible interval and then a seven-day interval criterion. Our models estimate for the 95% credible interval show that, for Austria's case notification, Models 1,2,3, and 5 predict a range that contains the actual peak period, while Model 4 predicts a slightly lower range. For Denmark, we observe that Models 1 and 3 provide accurate predictions in terms of their 95% credible interval, while the other models predict much lower intervals compared to the actual values for Denmark. Models 1, 3, and 4 have predicted intervals that align with the actual values for the Netherlands, whereas Models 2 and 5 do not. Models 1 to 3 provide good interval predictions that contain the 95% credible interval for Portugal, while Models 4 and 5 predict lower intervals.

We observe that, for the United Kingdom, only Model 1 predicts the 95% credible interval that contains the actual peak period, while the others predict lower intervals (Figure 4 B). Now, using a seven-day criterion, we find that the actual peak day for Austria occurs on day 264 of our time series. We identify that only Model 1 meets this criterion with a 257 day prediction. None of the models is able to provide an accurate prediction for Denmark and the United Kingdom. Models 1 predict 267 days for the Netherlands, whose actual value is on day 264; Model 3 predicts 293 for Portugal, which has its peak period on day 296. Inductively, we can conclude that Model 1 gives the best prediction across all five countries on average compared to the other four models.

However, this conclusion is limited without an accompanying statistical evaluation of their performances.

Statistical Evaluation

We calculate numerical estimates (A.I.C.c and Area) for the purpose of model comparison with reference to the empirical data. We focus on the second portion of the data (second wave) for which our models attempt to predict. We again use the Adjusted Akaike Information Criterion (A.I.C.c), which penalises the number of parameters in each model and also serves as a measure of goodness-of-fit. We also calculate the area (Area) between the models predictions and the empirical data of case notifications. In a similar manner to the first wave, using A.I.C.c for evaluation and comparison, we consider A.I.C.c for the incidence ( I i ) only and A.I.C.c for the incidence plus behaviour ( I i + X ). Our results show that Model 1 outperforms the other models for both incidence ( I i ) and incidence plus behaviour ( I i + X ) for all five countries, followed by Model 3 for Austria and Denmark. Model 4 outperforms Model 3 for the Netherlands and the United Kingdom regarding A.I.C.c incidence, but we observe that Model 3 outperforms Model 4 when considering both incidence and behaviour for Portugal, with a slight difference (Figure 5).

: Figure 5 summary: A bar chart comparing AICc values across five countries and five model versions, with each model split between calculations for incidence only and incidence plus behavior. Across all countries, the bars for incidence plus behavior consistently reach lower AICc values than those for incidence alone. The takeaway is that incorporating behavior into the models significantly improves predictive performance across all tested countries.

The area is calculated using the residuals between the model output and case notifications. We observe that Model 1 gives a low value for the area for all countries, indicating better performance, except for the Netherlands compared to other models. Model 4 gives a lower value compared to Model 1 for the Netherlands with respect to the area, indicating better performance. Model 3 also predicts the second lowest values of area for Austria, Denmark, and Portugal (Figure 6). Following these results, we can observe that Model 1 predicts the features of the second wave of S.A.R.S-CoV-2 the best amongst the five models evaluated in this study.

Figure 6 summary: A bar chart showing the prediction area between model incidence predictions and empirical case notifications for five different models across five countries. Model 5 shows the highest area in Austria and Denmark, while the other models maintain relatively lower and more consistent area values across all countries. The takeaway is that Model 5 generally exhibits the largest discrepancy between predictions and observed data, particularly in Austria.

Mitigating Behaviour

The behavioural model is fitted to the Oxford Stringency Index from day 1 to the cutoff (Figure 7, Table 1) for all the models of the five selected European countries. We observe that Model 1 fits well with all the countries, whereas Model 4 fits very poorly. Model 2 also fits poorly for all the countries except Portugal; likewise, Model 3 fits poorly for the Netherlands, Denmark, and the United Kingdom. Model 5 fits quite well in terms of the trajectory of the data for all the countries; however, the estimates are higher than the empirical data. We observe similar results for the prediction of the second portion of the data.

Figure 7 summary: Four line charts showing the estimated trajectories of mitigators over time for Austria, Denmark, the Netherlands, and the United Kingdom, compared against empirical data. In all countries, Model 5 (purple) significantly overestimates the mitigators over time, especially after the cutoff mark, while Model 1 (blue) and the empirical data (black) track more closely. Models 2, 3, and 4 generally stay low and underestimate the stringency index. The point is that Model 1 provides the most accurate fit to the data across different countries, whereas Model 5 consistently over-predicts mitigation measures.

Model 1 predicts the behaviour of mitigators in the population with accuracy, while the others perform poorly, except for Model 5, which gives a good prediction for Austria and Portugal compared to the other models. The A.I.C.c for prediction affirms our results on the behaviour (Figure 5), indicating that Model 1 is both a good fit and a good predictor of COVID-19 case notifications and changes in population behaviour.

Time series trajectories of COVID-19 case notifications for coupled behaviour-disease models

Discussion

In this study, we evaluated the explanatory and predictive abilities of five different coupled behaviour-disease models, from which we established that complex models with more parameters do not necessarily translate to the model's ability to fit accurately to empirical data and also produce good predictions. Whereas simplified models are able to do so. However, these additional parameters of more complex models do provide useful information about the dynamics of S.A.R.S-CoV-2 virus transmission when the models fit well, as our findings reveal in the case of three models (Model 1, 2, and 3). Our findings also reveal that these extra parameters may not be important for prediction purposes; yet, they indicate that certain key parameters are relevant for all the models and demonstrate a significant level of sensitivity.

The usefulness of a coupled behaviour-disease model can be categorised into two major parts: 1) model analysis in the absence of empirical data and 2) model analysis in the presence of empirical data. The former mostly takes a theoretical approach and may involve evaluating model parameters to assess their influence on the model, presenting general ideas that are not verifiable by the study itself since it is conducted in isolation from empirical data. The second part considers models formulated based on the dynamics of the available data.

These models can provide very useful insights into the model; however, this is dependent on the accurate estimation of the model parameters. Our results show that, while some parameters are of great importance to the formulation of coupled behaviour-disease models, others make very little contribution to the model's ability to estimate the empirical data. The simplest of the five models demonstrates that fewer parameters can achieve an accurate fit by obtaining good parameter estimates.

These estimates are also easily attainable compared to models with more parameters due to the level of complexity. They also show the possibility of producing good predictions from the simple model, while the more complex models do not produce good predictions for all five countries, but only for a subset. The structure of the model significantly influences its predictability, as our findings reveal when comparing Models 1 to 4 (similarity in structure) and Model 5, which predicts early and high second waves.

Despite the challenges of the more complex models, we identified that useful information can be extracted from the extra parameters, particularly those that fit well with the empirical model. For instance, this study shows that Models 2 and 3 indicate that socio-economic losses are an important consideration based on empirical data for Austria, the Netherlands, and Denmark. This implies that while mitigation strategies lead to socio-economic losses; as such, the economic impact must be accounted for in models for such countries; although excluding these parameters does not defeat the purpose of fitting and prediction, they provide more information. This consequently reveals that while certain additional parameters may be useful for one set of empirical data, they may not be useful for others. The COVID-19 pandemic was a global health crisis; the dynamics within each country differed, which implies that one model structure may not be a good fit for all. A case-by-case model can prove to be more useful; complex models with more parameters can be case specific and useful in producing desirable results.

We clearly see this play out with the empirical data for Portugal and the United Kingdom, where only Model 1 fits and predicts the case notifications accurately. Additionally, this model supports mathematical models with"on-off switches" for certain parameters to suit different scenarios based on the dynamics of empirical data.

This study contributes to the field of coupled behaviour-disease models; nevertheless, there are a number of limitations that could be addressed for further improvement. The first challenge we identify for this study is the fitting algorithm used for all the complex models. Although Approximated Bayesian Computation Sequential Monte Carlo (A.B.C S.M.C) is a good parameter estimation algorithm, it may not be suitable for different complex models.

Additionally, we assumed the same parameter range for the prior distribution to estimate the posterior distribution for all the models and countries. The different models may need to be sampled from different ranges. Another limitation is the number of fitting sequences that need to be carried out.

We considered 30 sequences (iterations), but due to the complexity of the models, more than 30 sequences of fitting may be needed to obtain a good fit and prediction. As we have already established, coupled behaviour-disease models in their simplest form are useful for fitting and prediction compared to the disease model. As such, for future work, we will consider all 13 European countries from our previous studies, perform stratification of the models with the empirical data to identify the best model-data outcome.

In doing so, we will explore fitting different models across different ranges, simulating more sequences using the A.B.C S.M.C algorithm and other fitting algorithms, such as Gradient Descent.

A further potential limitation of this study is the challenge of parameter identifiability of the models for various countries. We considered the Fisher Information Matrix (F.I.M), the Bayesian posterior analysis, and sensitivity analysis, but we could not concretely draw a conclusion on this problem. As a future study, an extensive investigation of the practical identifiability of the models for various countries can be undertaken by employing techniques like profile likelihood. This approach holds one parameter at a specific value while it optimises the other model parameters.

In summary, we have illustrated that increasing the complexity of coupled behaviour-disease models by the addition of more parameters does not automatically imply that such models can fit empirical data and produce good predictions. They can, however, provide useful information from the additional parameters included in the models. We also identified that complex models may be data specific and not a one-size-fits-all, revealing the usefulness of some parameters over others. A complex model, in their simplest form, can achieve a similar objective to what a complex model may not in terms of fitting and prediction, although the latter could perform better when the appropriate processes are followed.

Methods

We consider five coupled behaviour-disease models with a standard Sir model for the disease transmission of the S.A.R-CoV-2 virus. The disease system is compartmentalised such that we have susceptible, infected, and recovered (dead) groups. The Susceptible group consists of individuals who, when exposed to the virus, can become infected; the infected group is the set of individuals who have come into contact with the virus, have been infected by it, and may be able to transmit the virus to healthy susceptible individuals; and the recovered (dead) group includes individuals who have either recovered from the virus or died as a result of complications from the infection. The disease system remains the same for all five models, but they have different behavioural models and different coupling approaches.

We consider that the virus is affected by seasonality as a result of environmental or climatic changes. The seasonality of the virus is viewed as occurring in phases, as the case notifications reported by countries reveal. Thus, the seasonal Sir model equations are formulated as:

Equations 1, 2, 3 summary: This seasonal SIR model calculates the rates of change for the proportion of susceptible, infectious, and removed individuals over time. The process adjusts the baseline transmission rate using a cosine function to account for seasonal amplitude and phase, then tracks the flow of individuals from susceptible to infectious and finally to recovered based on recovery rates.

where S is the proportion of susceptible individuals (“susceptible”), I is the proportion of individuals who are both infected and infectious (“infectious”), and R is the proportion of individuals who are no longer infectious (“removed”). beta is the baseline transmission rate in the absence of non-pharmaceutical interventions, and gamma is the rate at which an infectious person recovers. b represents the amplitude of the seasonal force, and phi is the phase of seasonality.

Behaviour Models

The differences and complexity of the coupled behaviour-disease models considered in this study are accounted for by the differences in the behavioural models. Although all five models are built on the concept of evolutionary game theory by considering two behaviour types, their respective utility functions (expressions) differ. There is no doubt that the introduction and implementation of N.P.I's significantly influence the dynamics of the S.A.R.S-CoV-2 virus positively by decreasing infections and S.A.R.S-CoV-2 virus related deaths. In spite of these successes, due to the longevity of the survival of the virus, it has become evident that the use of N.P.I's may not be sufficient to mitigate the spread to low levels to enable the population to return to their pre-covid-19 state of living. N.P.I's such as social distancing, isolation, and quarantine, as well as practicing personal hygiene, have become costly for individuals, societies, and nations. This implies that the initial decision to adhere to the use of N.P.I's had to be reconsidered. At this point, the decision to mitigate or not to mitigate became a game where individuals had to weigh their options and make decisions. The average individual may want to make decisions that are most beneficial to themselves first, then to their household and community, by assessing the costs and benefits of available options (strategies).

We formulate the behavioural model by assuming that only two options (strategies) are available for the individual to choose from. The mitigating (M) strategy, where individuals support the use of N.P.I's, and the non-mitigating (N) strategy, where individuals do not support the use of N.P.I's. We suppose that an individual's decision on a strategy is not made in isolation from others in the community; that is to say, individuals imitate others in their decision-making since they aim to maximise their utility. We estimate a perceived payoff or utility for each strategy; mitigators payoff ( E M ) and non-mitigators payoff ( E N ). If the payoff for mitigators is greater than the payoff for non-mitigators, we expect a change of behaviour from non-mitigators to mitigators and vice versa. We represent the proportion of mitigators by the variable x; then the proportion of non-mitigators becomes 1 — x. Using the proportional imitation rule (replicator equation), we are able to model the changes in the proportions of both behavioural types, forming our behaviour model, which is then coupled with our seasonal Sir disease model to constitute the coupled behaviour-disease models.

Base Model

Models 1 to 4 have the same structure regarding how the disease system is coupled with the behaviour system, although the behaviour has different model structures (forms). Since the behaviour consists of contact precautionary measures, changes in the population influence the transmission potential of the S.A.R.S-CoV-2 virus by either enabling or reducing transmission.

Equations 4, 5, 6 summary: This system of differential equations calculates the rate of change for susceptible, infected, and recovered populations over time. It determines transmission by multiplying a time-varying seasonal factor and a mitigation reduction factor by the product of susceptible and infected individuals, while the recovery rate is based on a constant scaling factor applied to the infected population.

where x represents the proportion of mitigators in the population, and epsilon represents the strength of the use of N.P.I's by mitigators. We discuss the behaviour model's structure for models 1 to 4 in the following sections.

Model 1

This model is the same as in our study of coupled behaviour-disease models, and its ability to produce pandemic waves and demonstrate improved predictability compared to that of a disease model. Full details for such a model can be found in the paper, and the behaviour model equation is:

: Equation 8 summary: This expression calculates the rate of change of the behavior for a population over time. The computation multiplies the social learning rate, a quadratic term representing social learning dynamics, and the difference between a decaying sensitivity factor and the cost of maintaining measures.

where kappa is the social learning rate, lambda is the decay in sensitivity over time as the population learns more about the virus and becomes less afraid of it, and c is the cost of maintaining measures. The quadratic term x(1-x) represents social learning dynamics, where individuals sample others at some rate and change their opinions based on the utility differences in the brackets.

Model 2

This model is an extension of model 1, where we include socio-economic losses to influence the cost c of mitigation, and the accumulation of socio-economic loss structure changes to account for the strength of the use of N.P.I's (epsilon). This implies that the perceived payoff for mitigators E M becomes E M equals minus c times L. Following the same concept of evolutionary game theory with the imitation dynamics of the proportional imitation rule (E M minus E N equals I times e to the power of minus lambda t minus c times L), the model equation for the behaviour dynamics becomes:

Equation 9 summary: This expression calculates the rate of change in behavior dynamics over time. It multiplies a growth constant, the current proportion of the population, and a term representing the difference between a decaying payoff and the product of mitigation cost and socio-economic loss.

where L is a phenomenological representation of accumulated socio-economic losses that obeys

Equation 10 summary: This expression calculates the rate of change for accumulated socio-economic losses over time. The process scales the product of a proportional economic impact of mitigation, an error term, and input values, then subtracts the product of the social norm strength and the current losses.

where alpha represents a proportional economic impact of mitigation, and delta 2 represents the strength of the social norm.

Model 3

This model is a simplification of model 2, where we include socio-economic losses to influence the cost c, but the strength of mitigation represented in model 2 is defined as delta 2 = alpha epsilon as a proportional constant of the rate of mitigators. The model equation for the behaviour dynamics becomes:

Equations 11, 12 summary: This system of equations computes the rate of change for a population of mitigators and their associated socio-economic losses over time. The change in mitigators is determined by a growth constant multiplied by the current population, the remaining available population, and the difference between a decaying influence and the scaled economic loss, while the change in losses is the difference between the economic impact of the mitigators and the current loss value.

where delta 3 represents the economic impact of the mitigation.

Model 4

This model introduces social norms in the estimation of the respective perceived payoffs for each behaviour type: mitigators (E M equals c plus delta 4 times x) and non-mitigators (E N equals minus I times e to the power of minus lambda t plus delta 4 times open parenthesis 1 minus x close parenthesis); thus (E M minus E N equals minus I times e to the power of minus lambda t minus c plus delta 4 times open parenthesis 2 times x minus 1 close parenthesis). Following the proportional imitation rule, we obtain the model that accounts for the direct impact of social norms on mitigating strategies. Socio-economic losses are not included in this model.

Equation 13 summary: This expression calculates the rate of change of the proportion of mitigating strategies over time. The process involves multiplying a scaling constant and the current proportion by the difference between its complement and a composite value consisting of a decaying incentive, a constant cost, and frequency-dependent social processes.

where delta 4 of 2x minus 1 represents frequency-dependent social processes and delta 4 represents the strength of existing social norms of N.P.I's.

Model 5

This model takes a slightly different structure than the previous four models. The behaviour model still uses evolutionary game theory, but imitation does not follow the proportional imitation rule as in previous models; rather, players base their decisions on the difference in payoffs without necessarily sampling for imitation. The perceived payoff for mitigators is influenced by the proportion of mitigators at each time of decision making, that is, E M = -cx , and the ability of the population to learn about disease transmission also influences the perceived payoff for non-mitigators E N = kappa times -Ie superscript minus lambda t . The rate of the proportion of mitigators (x) per day is estimated by the difference between the two payoffs ( E M — E N = kappa Ie superscript minus lambda t — cx ). This is then coupled to the disease model by multiplying (1 — e superscript minus e x) with the transmission rate to account for the increase or decrease in transmissions as a result of behavioural changes. The coupled behaviour-disease model for this is:

: Equations 14, 15, 16, 17 summary: This system of differential equations computes the time rates of change for the proportions of susceptible individuals, infected individuals, recovered individuals, and mitigators. The transmission rate is adjusted by a scaling factor based on the mitigator proportion and a weight for non-pharmaceutical intervention strength, while the mitigator rate is determined by the difference between the weighted infection impact and a cost factor.

where epsilon 5 is the strength of the use of N.P.I's but no longer has to be constrained between 0 and 1 in this model.

Under-reporting

In a similar fashion to the study Ref. 38, we account for under-reporting in all five models since the chances of under-reporting data are very significant for the COVID-19 pandemic. Social distancing, isolation, and quarantine are effective mitigating strategies for the disease, but they present a disadvantage in terms of data collection, as many asymptomatic and symptomatic cases were not reported for verification. Additionally, the high stress level on hospital facilities and their personnel made it challenging to test all suspected cases for confirmation. To account for this case of under-reporting, we introduce a parameter eta , which is estimated for this purpose; the reciprocal is taken as the actual value, which is then multiplied by the estimated daily incidences.

Incidence and Parameterisation

The daily incidence is estimated by formulating the cumulative cases (C) for each day (t). The daily incidence is evaluated by taking the difference between successive cumulative cases and multiplying it by the under-reporting factor. Since we have two main model structures for coupling the disease system to the behaviour system, we have two expressions for the estimation of the cumulative cases. The first, which we represent as C 1 minus 4 , shows that for models 1 to 4, and C 5 represents that for model 5.

Equations 18, 19 summary: These expressions calculate the rate of change for cumulative cases over time for model structures one through four and model five. The process multiplies a transmission rate and a time-based cosine modulation by the product of susceptible and infected individuals, while further scaling the result by distinct factors involving an epsilon value and a behavioral variable.

Then the daily incidence for each day t is estimated as:

Equation 20 summary: This expression calculates the daily incidence for specific categories across different time periods. It computes the result by taking the difference between the current and previous cumulative case counts and multiplying that value by the inverse of a scaling factor.

This final expression is then fitted to the daily case notifications reported from each country under study for COVID-19 (I i reported of t) to estimate the model parameters for analysis and prediction.

We employ the same parameter estimation approach from Ref. of Approximated Bayesian Computation Sequential Monte Carlo (A.B.C S.M.C) for each model for the five selected European countries, chosen based on clear first and second waves.

We identify the number of parameters and initial conditions for each model that need to be estimated, as well as those needed for statistical analysis in model evaluation and performance comparison.

To estimate the parameter values for each model and for each country, we assume a range for each parameter, from which we create a uniform prior distribution for sample drawings.

Table 2 summary: Comparison of five models, their parameter counts and core assumptions. The models are similar in complexity, with parameter counts ranging from 12K for Model 1 to 15K for Model 2. All five models assume two behavior types and cost-induced utility functions. Model 2 is the most complex, incorporating socio-economic loss accumulation, social norms, and strategy imitation via sampling. Other models focus on specific subsets of these features, such as Model 3 and Model 4 which both have 14K parameters, or Model 5, which has 13K parameters and assumes non-mitigators can learn about the disease while the proportion of mitigators influences utility.

Table 3 summary: The parameter names, ranges, and sources for the study models. Key parameters include the transmission rate beta, ranging from 0.1 to 0.9, and the infectious period 1 over gamma, spanning 2 to 10 days. Model-specific parameters include intervention effectiveness epsilon for Models 1 through 4, ranging from 0 to 1, and epsilon 5 for Model 5, ranging from 0 to 1,000. Social norm strengths delta 3 and delta 4 for Models 3 and 4 both range up to 5. Other calibrated parameters include the proportionality constant alpha for economic impact of mitigation, ranging from 0 to 2, and the social learning rate kappa, ranging from 100 to 10,000 per day. The table also specifies initial conditions for population support for NPIs x zero, ranging from 0 to 1, and initial infection and incidence rates that are scaled by alpha and the country population size.

Parameter Identifiability and Uncertainty

The identifiability of a model is simply the uniqueness of model inputs producing unique outputs. It takes two forms: structural and practical. Structural identifiability can be considered to hold when the model parameters cannot be replaced by fewer parameters. Here, Model 2 is structurally unidentifiable since the parameters alpha, epsilon can be simplified into one parameter as we have in Model 3 through parameter substitution.

Since parameters are estimated from empirical data, we place more emphasis on practical identifiability for this study. Several methods, such as profile likelihood, the Fisher Information Matrix (F.I.M), Bayesian posterior analysis, and sensitivity analysis, can be employed to detect practical identifiability. We consider the last three for the analysis of our models.

The Fisher Information Matrix measures the amount of information about a set of parameters that can be extracted from a probability distribution given observable data.

Math summary: This expression calculates the Fisher Information Matrix. It computes a sum of the products of the transpose of the output derivatives with respect to parameters, the inverse of the measurement noise covariance matrix, and the output derivatives.

where (t i)^{-1} is the covariance matrix of the measurement noise. The covariance matrix is simply the inverse of the F.I.M, that is, Cov({p}) = F.I.M^{-1} . The Bayesian posterior analysis assesses the distribution of the posterior to make judgments on practical identifiability, while the sensitivity also provides additional support for the change in output with respect to changes in model inputs. A significant difference in parameter values corresponding to low changes in output can be classified as having weak identifiability or non-identifiability. Additionally, we consider the correlation matrix for the posterior distribution to assess the relationship between the model parameters. Details of these analyses are presented in the supplementary appendix.

References

1. Kretzschmar, M. & Wallinga, J. Mathematical models in infectious disease epidemiology. In Modern infectious disease epidemiology: Concepts, methods, mathematical models, and public health, 209–221 (Springer, 2009).