Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Adding a reaction-restoration type transmission rate dynamic-law to the basic SEIR COVID-19 model

  • Fernando Córdova-Lepe ,

    Contributed equally to this work with: Fernando Córdova-Lepe, Katia Vogt-Geisse

    Roles Conceptualization, Formal analysis, Methodology, Writing – original draft, Writing – review & editing

    Affiliation Facultad de Ciencias Básicas, Universidad Católica del Maule, Talca, Chile

  • Katia Vogt-Geisse

    Contributed equally to this work with: Fernando Córdova-Lepe, Katia Vogt-Geisse

    Roles Formal analysis, Funding acquisition, Methodology, Visualization, Writing – original draft, Writing – review & editing

    katia.vogt@uai.cl

    Affiliation Facultad de Ingeniería y Ciencias, Universidad Adolfo Ibáñez, Santiago, Chile

Abstract

The classical SEIR model, being an autonomous system of differential equations, has important limitations when representing a pandemic situation. Particularly, the geometric unimodal shape of the epidemic curve is not what is generally observed. This work introduces the βSEIR model, which adds to the classical SEIR model a differential law to model the variation in the transmission rate. It considers two opposite thrives generally found in a population: first, reaction to disease presence that may be linked to mitigation strategies, which tends to decrease transmission, and second, the urge to return to normal conditions that pulls to restore the initial value of the transmission rate. Our results open a wide spectrum of dynamic variabilities in the curve of new infected, which are justified by reaction and restoration thrives that affect disease transmission over time. Some of these dynamics have been observed in the existing COVID-19 disease data. In particular and to further exemplify the potential of the model proposed in this article, we show its capability of capturing the evolution of the number of new confirmed cases of Chile and Italy for several months after epidemic onset, while incorporating a reaction to disease presence with decreasing adherence to mitigation strategies, as well as a seasonal effect on the restoration of the initial transmissibility conditions.

Introduction

A novel Coronavirus (SARS-CoV-2) emerged from the city Wuhan in China in December 2019 and has caused a devastating public health impact across the world [1]. As of June 28, 2021, COVID-19 has caused over 180 million confirmed cases and over 3.5 million deaths worldwide [2]. The curves for daily confirmed new cases of COVID-19 in different countries present a high variability in their geometric forms. Every such curve shows a sequence of outbreaks and valleys when observed over time, while the sharpness of the outbreaks and the length of the valleys can vary [3].

In fact, in [3] we can find figures that show the daily number of confirmed new cases (7 day moving average), in particular from March, 2020, to January, 2021, for different countries. Most of them have experienced a second wave and others show even a third. Some European countries show a sharp first outbreak followed by a plateau of low height (that lasted several months) before the second wave, whose peak out-measures the peak of the first outbreak by at least three-fold. On the contrary, several South American countries presented an initial exponential phase of several months, soon after which a second wave of similar peak height as the first occurred. Finally, there are countries in which the curve of new cases of COVID-19 has shown extreme behaviors in some part of its evolution as compared to European and South American countries. For instance, Czech Republic has experienced a very weak first outbreak followed by a low plateau lasting for months; Iran did have a more pronounced first outbreak, but it was followed by a plateau of important height; and Indonesia experienced an almost unnoticeable first outbreak that is part of an exponential growth phase when zoomed out, which lasts several months.

Classical compartmental models based on the classical Kermack & McKendrick SIR model [4] with constant parameters often used to model epidemics do not reflect the behavior over several months described above. Neither the (β, γ)SEIR model for a population of size N– which is compartmentalized into susceptible (S), exposed (E), infectious (I) and removed (R)– given by (1) nor extensions of it have been efficient in adjusting the data well beyond the first epidemic outbreak when considering the transmission rate β and the removal rate γ constant. This is due to the unimodality of the active-infected-curve those models provide, i.e. one bell-shaped infection curve and an epidemic growth limited by the proportion of susceptible individuals [5, 6].

In general, the epidemiological data series do not reflect that the percentage variation of susceptibles per proportion of active cases, i.e. |S′/S|/(I/N), is approximately constant, as is assumed in the aforementioned classical epidemiological models. In fact, there exists literature that evidences the changing temporal behavior of disease transmission in epidemic or pandemic situations [715]. In particular, there are studies using mathematical models– some aiming to understand COVID-19 transmission– that include the decrease in the transmission rate [7, 9, 1620], and some incorporating human behavioral factors as part of the cause for a temporal change in transmission. For instance, in [7] the behavior of the transmission rate– providing exponential saturation for a large number of infectives– for three consecutive months is shown for four geographical settings: worldwide, United States, Russia and Canada. For each setting, the trend is a decrease in the transmission rate for then stabilizing at a value several times lower than initially. They as well include a parameter representing mask wearing within their transmission rate. Finally, this study also shows that the recovery rate R′/I is for each setting much more stable than the transmission rate. Supporting the idea that the lack of efficiency for the classical compartmental models to adjust well to data is due to that the transmission rate β is assumed constant over time.

In this article we attempt to break the unimodality of the active-infected-curve of the classical epidemiological models. We introduce a novel way to model the behavior of the transmission rate β, considering a balance equation between a reaction rate and a restoration rate; and including the resulting dynamic law for the transmission rate into the classical SEIR model. The paper is structured in the following way: In the next section The Transmission Rate we provide some understanding about the transmission rate of infectious diseases. In the section The βSEIR Mathematical Model we introduce and describe a new basic model, which we call βSEIR model, by adding to the classical SEIR the aforementioned dynamic law for the transmission rate, and show some mathematical analysis. In the section Numerial Results we provide simulation results; and finally, the last section contains the discussion and conclusions.

The transmission rate

There exist at least two groups of epidemic control measures. The first, aims to reduce the population that is being hit by the disease, i.e., the susceptible population. Such measures are, for example, vaccination or limiting the mobility of individuals. The second, intents to reduce the force of infection that is defined as the product of three quantities: number of close contacts per unit of time of a susceptible individual (pC), probability that a close contact is with an infectious individual (pI) and probability of transmission given a close contact with an infectious (pT). Reducing pI results in that per unit of time, there exist fewer active cases in the population, which is accomplished by eradication, i.e., removal from the system (e.g., slaughtering of sick animals which is widely used in animal epidemics, or banishing infectious individuals as was done aforetime), or by applying actions for a rapid recovery. Notice that the product of pC and pT is called the transmission rate and is usually denoted by β (see e.g. Eq (1)). Hence, the objective of most mitigation strategies that aim to reduce the force of infection, aim to reduce the transmission rate (lower β) by either increasing physical distance and hence reducing the number of close contacts (lower pC) or blocking the transfer of pathogens to a new host (lower pT). There are secondary measures such as for example reducing population movement (which is not reducing physical distancing nor blocking transmission), which make close encounters less likely.

When a highly transmissible disease with high mortality or morbidity invades a population of mostly susceptible individuals, and a vaccine is not in sight in the short term (as was initially the case for COVID-19), health authorities’ only way for reducing morbidity and mortality is mitigation, while the general population’s duty is to comply to the new norms and desired behavior. In other words, the efforts are put into reducing the transmission rate β. In this sense, β is a time-spacial dependent variable, i.e. it changes according to time and location. Also less controllable physical-environmental aspects in relation with the bio-chemical characteristics of the pathogen may influence it (but we consider those factors constants in this study). It is also worth mentioning that in general populations, individuals may live and participate in several cultural regions, which may also determine the variability of the transmission rate. We will assume in this study that the population stays within its territory during the time horizon studied, behaving homogeneously in this respect. We suppose that the disease studied will have a base line transmission rate β0, that we will call natural transmission rate, measured for a population that is initially free from the disease and does not consider any mitigation strategies or personal protective measures.

One of the characteristics of COVID-19 was that it has had a large media coverage since the first confirmed cases appeared in December 2019. This provoked sentiments of fear in the general population and played an educational role for pandemic preparedness (e.g. global media emphasizing on washing hands and physical distancing) before the imminent arrival of COVID-19 in many countries. It shaped how countries would confront COVID-19 right from the appearance of their first confirmed case and even before. The reproduction numbers corresponding to different geographical locations was most likely to be between 2 and 5 [21, 22], and was shown to rapidly decrease during the first weeks of the pandemic [1214], however reaching values above one that nevertheless allowed COVID-19 expansion.

To consider this decreasing effect, many authors assume an exponential decay of the transmission rate for a certain amount of time, for example, in [9, 18] they assumed in their continuous time model β(t) = β0exp(−b0 t), t ≥ 0, or in in [19, 20] they defined βk = β0 ak, 0 < a < 1, k ≥ 0, where k is a day-counting integer in their discrete time model. In order to extend the horizon of validity of their model some authors consider an exponential decrease from β0 to a minimum positive bound [15, 23]. To understand the rate of decrease from that baseline natural β0 and identify a time varying β(t) transmission rate, some researchers use mathematical expressions and the data on active cases I(t) and removed cases R(t) in a population of constant size N; for instance, one can use β(t) = −NS′/(SI) as in [7], which can be obtained from a SIR model and approximating the derivatives using the finite differences on one week running averages; or one could use β(t) = γ + I′/I at the beginning of an outbreak, assuming SN in the SIR model.

More time-varying transmission rates have been considered within mathematical models. For instance, the authors in [9] capture the early decreasing trend of COVID-19 in Malaysia using a time varying exponential decay log function β(t) = (1 − p)t for the transmission rate in an SIR model, that uses a fractional term z to measure the effectiveness of interventions and a proportion p to account for depletion. In the literature there are studies that, in order to capture realistic disease transmission, assume non-linear functions of S and I governing the force of transmission, as for instance in [16], where the force of transmission they use in an SIR-type model depends on the product of fractional powers of S and I. They use the model to fit COVID-19 data of Italy, Germany, France and Spain.

In [24] the authors include in an SIR model time-varying transmission rate, assuming that the probability of transmission of a susceptible is βλt(It/N), where λt(⋅) is a random variable, they refer to this model as a Spatial-SIR model. In [25] the authors assume a contagion rate as a sum of a base-line transmission rate and a component that satisfies a first order linear differential equation to represent the effect of non-pharmaceutical interventions (NPIs).

Behavioral factors represented in the transmission rate are also considered by more authors in order to represent the changing dynamics of the transmission rate. For instance in [26] the authors study a model that incorporates a non-constant transmission rate β(M) that depends not only on the current number of infectious individuals but on M, representing an information index that summarizes the current and past history of disease prevalence. Part of their results discuss that social behavioral change may trigger oscillations. The study in [27] extends an SIR type model defining a transmission rate that captures the impact of school and workplace closure through a function of time. The changing behavior of this function is based on Imitation Dynamics [28] and describes population-level support dynamics for closure. The article in [29] also uses Imitation Dynamics and studies a population in which individuals develop and learn a behavior of mutual protection.

The novelty of this article is that, at each moment in time we consider the variation of the transmission rate to be given according to a balance equation between two opposite thrives: a reaction rate and a restoration rate.

In what follows we are going to justify the functional forms of the reaction and restoration rates, as well as present the βSEIR model that incorporates the dynamic law for β. Further, we are going to analyze its effect on the shape of the main epidemiological curves.

The βSEIR mathematical model

In this section we present the dynamic law of the transmission rate β in order to introduce the βSEIR model and some mathematical analysis.

Transmission rate β

In this subsection we derive the form of the transmission rate at which susceptible individuals become exposed upon contacts with infectious. The only infectious class of the model is the I class. We model the case of an infectious disease transmitted directly from person to person, and assume that at the beginning of an outbreak, the appearance of first cases do not provide a reason for alarm and panic. Hence, initially the disease propagation is due to a high natural-transmission rate intrinsic to the population while no interventions to mitigate disease spread are in place. We call this natural-transmission rate β0. In general, β0 makes the disease expand rapidly, producing a fast initial increase in new cases.

We present a new form for the transmission rate, which is represented by a dynamic, time-dependent quantity that is governed by a balance equation between a reaction rate, g[t, I], and a restoration rate, f[t, β], i.e. the proportion by which the transmission rate decreases and increases per day, respectively, represented by (2)

Next we justify the introduction of a reaction rate and a restoration rate and propose a functional form for each.

Reaction rate.

During severe epidemic outbreaks that attract huge public attention and media coverage due to for instance a high morbidity and/or mortality in the population, a steady increase in the implementation of measures that aim to reduce the transmission rate can be observed. As long as there is no licensed vaccine or treatment, these measures are mainly based on non-pharmaceutical interventions. Here we are interested in those directed to diminish the factors pC and pT, whose product defines β, such as social distancing measures (that reduce the number of close contacts between people: large-scale or home quarantines, workplace non-attendance, travel restrictions, prohibition of social gatherings, school closures, etc.) that reduce pC, or blocking measures (that, given e contact, reduce the pass of the pathogen: hand-washing, respiratory etiquette, face-masks usage, etc.) that reduce pT, see [30, 31]. When fear governs the population, people react, complying with mandatory measures or adopting self-protecting measures to avoid infection [32]; we note that risk communication is also a factor to support the general public response [33]. It is to expect that the higher the severity of the disease is, the more effort is put into mitigation. In this sense, individuals’ reactions produce a decrease in the transmission rate from its initial natural-transmission rate value, β0. Hence, we define a reaction rate that we denote g[t, I], which is non-negative and positively correlated with the number of active cases I(t), in a way that it increases when I(t) does. It may also depend on other circumstantial conditions of the moment relative to the population. These we will discuss further later in the text.

We assume in this article, that the reaction rate, follows the Michaelis-Menten model [34] describing the reaction to the presence of the infectious (active cases) I(t) at any moment in time, i.e. we define (3) where we call μ(t) the reaction coefficient, which is a non-negative function that represents the daily maximum possible reduction at time t, and Im > 0 is the half-saturation constant, i.e. is the number of active cases, where the reduction is half-maximal. Notice that the parameter Im characterizes the population, i.e. it determines its sensibility to react to active cases.

Restoration rate.

It is important to observe that upon the appearance of a reaction rate there exist socio-environmental factors that tend to restore the transmission rate to the level observed at the beginning of the pandemic, i.e. to β0 [27, 35]. When in a certain location the health authorities cease to impose protective measures, e.g. the use of face masks, and individuals lost their initial fear, then the transmission rate that had been reduced due to these measures no longer stays low and returns to its natural level. Therefore, we introduce a restoration rate that at each instant t, t > 0, is responsible for an increase in the transmission rate. The restoration rate that we denote by f[t, β], is a non-negative function that correlates directly with the distance between β(t) and its natural value β0, as presented in the following equation: (4) where is a positive exponent. We call ν(t) the restoration coefficient, which is a non-negative function that regulates the daily form of the restoration rate. Also, f is an increasing function of |β(t) − β0|, and f[t, β0] = 0, which means that if the transmission rate β(t) reaches at a certain time point t its natural value, then there is no deviation to restore.

The βSEIR model: Formulation and analysis

We incorporate the differential Eq (2) to the classical SEIR model (1) obtaining the βSEIR model, given by the following system of equations (5) with some non-negative initial conditions β(0) = β0, S(0) = S0, E(0) = E0, I(0) = I0, R(0) = R0, and f[t, β], g[t, I] as in Eqs (4) and (3) respectively. Table 1 describes the variables and parameters of the model.

thumbnail
Table 1. Description of variables and parameters from the model in system (5).

https://doi.org/10.1371/journal.pone.0269843.t001

In the following we show that 0 ≤ β(t) ≤ β0, for all t ≥ 0 and β0 > 0. Observe that β′ ≥ − μ(t)β, and hence due to Grönwall’s inequality (see [36]) we can conclude that . Therefore, β(t) is non-negative for any non-negative initial condition β0. We also observe that β(t) ≡ β0 is an equilibrium solution of the first equation in system (5) as long as I(t) ≡ 0. In case I(0)>0, i.e. disease is present in the population, β ′ (0) < 0 and hence β(t) decreases initially. Since β′ ≤ ν(t){|ββ0|/β0}α β, using Grönwall’s inequality and the fact that the solutions to the differential equation β′ = ν(t){|ββ0|/β0}α β that pass through points with are increasing and bounded by β0, we obtain that β(t)≤β0 for all t ≥ 0 and β0 > 0, I0 > 0.

Just as for the classical SEIR model, we observe that the epidemiological state variables of the model in system (5) remain positive for positive initial conditions and are bounded by the total population size N. Also, adding the second, third and fourth equations in system (5) together we obtain (S + E + I)′ = −γI < 0 whenever I > 0. Hence, S + E + I is a non-negative smooth decreasing function, and therefore exists. On the other hand, the derivative of any smooth non-negative decreasing function must tend to zero, and hence , which implies . Similarly, by adding the second and third equation in system (5) together, we can prove that . Since the limit of S + E + I exists when time tends to infinity, we can then conclude that . The behavior of R can be obtained from N = S + E + I + R and we obtain . Hence, the long term behavior of the model we present holds (i.e. the limit of the epidemiological state variables exist for infinite time), just as for the classical compartmental epidemic models SIR or SEIR with constant transmission rate β [37]; in particular we have shown that the disease in the long term goes extinct.

Finally, as for classical models we can obtain a threshold condition that determines an initial epidemic outbreak. Given β0 > 0 and considering that β(t)≤β0 for all t ≥ 0 as well as that the state variables are bounded by N, we have that (I + E)′ = γI{(β/γ)(S/N) − 1}≤γI{(β0/γ) − 1}. Hence, we define the basic reproduction number [38, 39] as . This way, if , then the curve representing the infectious population is decreasing to zero and there is no epidemic outbreak. On the other hand, if , then (E + I) increases initially when we assume SN and β(0) = β0, and increases as long as (β(t)/γ)S(t)/N > 1 holds. Notice that, in this case, the curve of infectious individuals may be increasing at several time intervals depending on the behavior of the function β(t), and not only depending on the ratio of susceptible individuals in the population that is decreasing according to the second equation in system (5). In Section 1 we will call the quantity (β(t)/γ)S(t)/N the Effective Reproduction Number, , and its dynamics will determine disease dynamics.

In the context of our study, we consider short-medium term scenarios under an epidemic situation, i.e. when . We also consider no replacement of susceptibles, due to the time frames we describe and assuming that: first, infected individuals acquire some kind of protective immunity for several months after infection [40]; second, we do not consider multiple variants that may provoke new infection. In particular, we will point out the differences compared to the classical SEIR model with constant β, in which E and I are variables that describe the known unimodal behavior of one bell-shaped curve.

Our model generalizes an idea presented in [41], where the authors consider an SIR-type model with variable transmission rate of the form β(D(⋅)), in which the time dependent function D(⋅) represents social distancing that individuals in the population maintain to each other, governed by the dynamics of D′ = −λ1(DD*) + λ2 I/N, where λ1 and λ2 are positive constants, and D* is the culturally dependent natural social distance, to which D(⋅) converges in the absence of disease (I = 0).

Notice that if in their model with positive parameters, we observe that satisfies the equation , with and . Observe that regardless of the convergence of , we obtain when t → ∞.

Numerical results

In this section, we present simulations of disease dynamics assuming: first, constant reaction (μ(t)) and restoration (ν(t)) coefficients; second, we extend our results to consider time-varying reaction coefficients, representing a diminishing mitigation effect of government responses; third, an additional seasonal effect on the restoration coefficient and we show that our model is capable of capturing real COVID-19 disease data. For the simulations we use the Python programming language [42].

The autonomous βSEIR model: Constant reaction and restoration coefficients

We consider in this subsection the autonomous βSEIR model from system (5), assuming a constant reaction coefficient, μ(t) = μ; i.e. we suppose that the reaction to reduce the transmission rate just depends on point prevalence levels and not on an additional time factor (see Eq (3)); and also assuming a constant restoration coefficient, ν(t) = ν; i.e. the regulation on the restoration rate depends only on the deviation of the transmission rate from its natural value β0 and not on external temporal factors (see Eq (4)). We present qualitative results of our model through simulations considering an 18 months time span, and present in each of the Figs 13 five subplots that represent the dynamics for: restoration f(t, β) and reaction rates g(t, I); the transmission rate β(t); the effective reproduction number ; the new confirmed cases eE(t); and finally we illustrate the cumulative number of cases. The effective reproduction number is a time-varying threshold quantity– defined by – such that the number of cases increase while , reach a peak when and decrease when [43, 44], and in particular .

thumbnail
Fig 1. Simulations for a high restoration coefficient ν = 0.8.

The first subplot illustrates the restoration rate f(t, β) (dotted) and the reaction rate g(t, I) (solid). The remaining subplots show: The transmission rate β(t), the effective reproduction number , the new confirmed cases eE(t), and the cumulative cases E(t) + I(t) + R(t). The reaction coefficient in each subplot are chosen to be μ: 0.3 (blue); 0.4 (green); 0.5 (orange) and 0.6 (red), and the remaining parameter values and initial conditions are as in Table 2.

https://doi.org/10.1371/journal.pone.0269843.g001

thumbnail
Fig 2. Simulations for a medium restoration coefficient ν = 0.5.

The first subplot illustrates the restoration rate f(t, β) (dotted) and the reaction rate g(t, I) (solid). The remaining subplots show: The transmission rate β(t), the effective reproduction number , the new confirmed cases eE(t), and the cumulative cases E(t) + I(t) + R(t). The reaction coefficient in each subplot are chosen to be μ: 0.3 (blue); 0.4 (green); 0.5 (orange) and 0.6 (red), and the remaining parameter values and initial conditions are as in Table 2.

https://doi.org/10.1371/journal.pone.0269843.g002

thumbnail
Fig 3. Simulations for a low restoration coefficient ν = 0.2.

The first subplot illustrates the restoration rate f(t, β) (dotted) and the reaction rate g(t, I) (solid). The remaining subplots show: The transmission rate β(t), the effective reproduction number , the new confirmed cases eE(t), and the cumulative cases E(t) + I(t) + R(t). The reaction coefficient in each subplot are chosen to be μ: 0.3 (blue); 0.4 (green); 0.5 (orange) and 0.6 (red), and the remaining parameter values and initial conditions are as in Table 2.

https://doi.org/10.1371/journal.pone.0269843.g003

Additionally, the constant restoration coefficient ν takes in Figs 13 the values 0.8 (high), 0.5 (medium) and 0.2 (low), respectively, representing high, medium and low rates to restore transmission levels due to the urge to return to the natural transmission rate β0. For each constant restoration coefficient value, we choose within each figure the constant reaction coefficient μ to take the values 0.3 (low, in blue), 0.4 (medium-low, in green), 0.5 (medium-high, in orange) and 0.6 (high, in red), representing different constant levels of the daily maximum reaction to disease that reduces the transmission rate. The remaining parameter values and initial conditions used in the figures are described in Table 2. Note that the basic reproduction number obtained from the values β0 = 0.65 and γ = 1/14 from Table 2 is . For illustration purposes, we use a value significantly larger than one.

thumbnail
Table 2. Initial conditions and parameter values and their units, for the simulations in Figs 16.

https://doi.org/10.1371/journal.pone.0269843.t002

Fig 1 considers a high restoration coefficient, ν = 0.8. Observe that for high reaction coefficients, e.g. μ = 0.6 as well as for small reaction coefficients, e.g. μ = 0.3, the restoration rate f(t, β) and the reaction rate g(t, I) are very similar to each other. Initially the restoration rate is slightly smaller than the reaction rate, producing that β′(t)<0 (from the first equation in system (5)) and therefore being the reason for the initial decrease of in the transmission rate. After that drop, the reaction and restoration rates are almost equal, producing a plateau in the transmission rate due to β′(t)∼0, and subsequently, a slightly larger restoration rate than reaction rate produces an increase in the transmission rate, which converges to its original natural transmission rate value β0. Additionally, observe that for all values of the reaction coefficient μ, the duration of the plateau in the transmission rate is directly correlated with the value of μ: the higher the μ value, the longer its duration. Note that, while the transmission rate is at the plateau– i.e. β(⋅) behaves similar to constant– when equaling the first equation in our βSEIR model (system (5)) to zero, we conclude from f(t, β) = g(t, I) when α = 1 that

Hence, if ImI, then ββ0λ. Indeed, notice that the value of β0λ for μ = 0.3, 0.4, 0.5 and 0.6 in Fig 1 are respectively 0.406, 0.325, 0.244 and 0.163, which correspond very closely to the plateau levels of the respective transmission rates. For all μ cases, the effective reproduction number shows a decreasing shape, staying above one at least during the first months of a pandemic (the horizontal dotted line represents Re = 1 in the figure). While the effective reproduction number stays above one, one can observe that the number of new confirmed cases increases, reaching a peak when the effective reproduction number reaches the threshold value one. The lower the reaction coefficient value is, the larger is the transmission rate throughout the epidemic, which produces sooner and larger epidemic peaks of new confirmed cases, as well as a rapid increase in the number of cumulative cases, reaching quickly a number close to the final number of infected individuals throughout the whole epidemic. This happens shortly after the effective reproduction number reaches the value of one, and is due to a small number of susceptibles remaining in the population at that time.

Fig 2 considers a medium restoration coefficient, ν = 0.5. We observe that for a high reaction coefficient μ = 0.6, the restoration (f(t, β)) and reaction (g(t, I)) rates show an initial oscillation and differ from each other more clearly, being initially the reaction rate larger than the restoration rate and subsequently both intersecting several times. This produces an oscillatory behavior in the transmission rate due to the β equation in system (5), attaining the transmission rate a local maximum or minimum value each time the reaction and restoration rates intersect, i.e. f(t, β) = g(t, I). On the contrary, in the case of a low reaction coefficient value, e.g. μ = 0.3, we observe very similar restoration and reaction rates (blue curves) just as in Fig 1. We can also observe from Fig 2 that a high (red) or medium-high (orange) reaction coefficient μ, drives the effective reproduction number below one much faster than in was observed in Fig 1 (in the case of a higher restoration coefficient) and way before the cases for medium-low (green) and low (blue) reaction coefficient values. It is interesting to see that in the cases of higher reaction coefficients (red and orange), after the initial fast drop of the effective reproduction number, during the remaining time pictured it oscillates around one. When comparing the curves of the new confirmed cases and cumulative cases for these two reaction coefficient values, with the cases of low and medium-low reaction coefficients, we notice that the number of cases is way higher for the latter. Additionally, one bell-shaped curve of new confirmed cases is being observed for smaller reaction coefficient values (blue, green), since the effective reproduction number only manages to cross the threshold (see dotted line) once due to the small number of susceptibles remaining after that first large peak. On the contrary, an oscillatory behavior is seen for higher μ values, obtaining small epidemic peaks each time the effective reproduction number reaches one in a decreasing manner.

In other words, the unimodality of the behavior for new confirmed cases observed when μ is low (blue) or medium-low (green) is broken for high (red) and medium-high (orange) μ values, which was not seen in Fig 1. One can also observe that for higher μ values (red and orange), the increase in the cumulative cases is close to linear in the time-frame pictured, as opposed to the rapid increase in cumulative cases for smaller reaction coefficient values (blue, green).

Fig 3 shows the case of a low restoration coefficient ν = 0.2. We observe oscillatory behavior in the restoration and reaction rates, producing an oscillatory behavior in the transmission rate, and hence also in the effective reproduction number and in the number of new confirmed cases. We observe that for high reaction coefficient values, e.g. μ = 0.6, the absolute difference between the reaction and restoration rates are larger and their intersections occur sooner, and hence the oscillations in the transmission rate have a higher amplitude and their local maxima occur sooner, than in the case of lower reaction coefficient values, e.g. μ = 0.3. Also, transmission rates with higher amplitudes produce less pronounced peaks in the oscillatory behavior of the number of new confirmed cases. Additionally, sooner starting oscillations correspond to higher values of the reaction coefficient μ. The cumulative cases show a near to linear increase, where smaller slopes correspond to higher reaction coefficient values.

A non-autonomous βSEIR model: Time-varying reaction coefficient representing a diminishing mitigation effect

It is to expect that, when a disease enters a population, the reaction coefficient μ(t) right after the onset increases quickly, reaches a maximum value μ0, and then decreases due to many factors. This can be deduced (at least) from two sources: (a) The data and information provided by the Oxford COVID-19 Government Response Tracker (OxCGRT) and the time curves defined by the Stringency Index [45, 46]. This index records the intensity of several government responses combined, such as containment and closure policies by country. In general, we observe that, first, the time curves representing the stringency index rise, but then follow a decreasing behavior due to local socioeconomic reasons [4752]. (b) Studies in behavioral science explain the public fall of adherence to mitigation (distancing) measures, as a daily average compliance curve in times of COVID-19 shown for instance in [53] or discussed in the conclusions in [54]. Additionally, it is known that information-based interventions positively impact compliance with mitigation restrictions, such as keeping a certain social distance, decreasing the number of times individuals go out and their time spent outside. Nevertheless, if individuals have been restricted for a prolonged period of time, compliance with such mitigation strategies decreases [55]. In this manuscript we refer to a time period right after disease onset, where there is no vaccination yet available.

We can include such a behavior– representing a diminishing mitigation effect of restrictive measures– through a time-varying reaction coefficient μ(t), for instance, given by the following Eq (6), (6) with 0 < μ0 < 1.

Notice that h(⋅) is a non-negative, unimodal function such that h′(b) = 0 and h(b) = 1, i.e. it achieves its maximum value at t = b and then decreases at a rate that depends on the parameter a.

In each of the Figs 46 we show for high, medium and low constant restoration coefficients ν, respectively, the curves for the restoration (f(t, β)) and reaction (g(t, I)) rates, the transmission rate β(t), the effective reproduction number, the new confirmed cases and the cumulative cases, for a forgetting curve h(t) of the form given in Eq 6, with a = 40, b = 90, such that the maximum occurs at t = 90 days. Within each plot we present curves for different maximum values μ0 of the now variable reaction coefficient μ(t) = μ0 h(t). The other parameters used for these figures are given in Table 2.

thumbnail
Fig 4. Simulations for a high restoration coefficient ν = 0.8.

The first subplot depicts the shape of the compliance curve h(t) = (a + 2b)t/(t2 + at + b2), with b = 90, a = 40, and the second plot in the first row illustrates the restoration rate f(t, β) (dotted) and the reaction rate g(t, I) (solid). The remaining subplots show: The transmission rate β(t), the effective reproduction number , the new confirmed cases eE(t), and the cumulative cases E(t) + I(t) + R(t). The maximum value μ0 of the reaction coefficient (μ(t) = μ0 h(t)) used in each subplot is: μ0: 0.3 (blue); 0.4 (green); 0.5 (orange) and 0.6 (red), and the remaining parameter values and initial conditions are as in Table 2.

https://doi.org/10.1371/journal.pone.0269843.g004

thumbnail
Fig 5. Simulations for a medium restoration coefficient ν = 0.5.

The first subplot depicts the shape of the compliance curve h(t) = (a + 2b)t/(t2 + at + b2), with b = 90, a = 40, and the second plot in the first row illustrates the restoration rate f(t, β) (dotted) and the reaction rate g(t, I) (solid). The remaining subplots show: The transmission rate β(t), the effective reproduction number , the new confirmed cases eE(t), and the cumulative cases E(t) + I(t) + R(t). The maximum value μ0 of the reaction coefficient (μ(t) = μ0 h(t)) used in each subplot is: μ0: 0.3 (blue); 0.4 (green); 0.5 (orange) and 0.6 (red), and the remaining parameter values and initial conditions are as in Table 2.

https://doi.org/10.1371/journal.pone.0269843.g005

thumbnail
Fig 6. Simulations for a low restoration coefficient ν = 0.2.

The first subplot depicts the shape of the compliance curve h(t) = (a + 2b)t/(t2 + at + b2), with b = 90, a = 40, and the second plot in the first row illustrates the restoration rate f(t, β) (dotted) and the reaction rate g(t, I) (solid). The remaining subplots show: The transmission rate β(t), the effective reproduction number , the new confirmed cases eE(t), and the cumulative cases E(t) + I(t) + R(t). The maximum value μ0 of the reaction coefficient (μ(t) = μ0 h(t)) used in each subplot is: μ0: 0.3 (blue); 0.4 (green); 0.5 (orange) and 0.6 (red), and the remaining parameter values and initial conditions are as in Table 2.

https://doi.org/10.1371/journal.pone.0269843.g006

We present in Fig 4 the case of a high constant restoration coefficient ν = 0.8. We observe that the restoration and reaction rates are very similar regardless of the μ0 value. Despite their similarity, initially, up to day 90, a slightly higher reaction than restoration rate produces a sharp decrease in the transmission rate. The decreasing shape of the forgetting curve starting on day 90, immediately reduces the reaction rate to slightly below the restoration rate, producing an increase in the transmission rate, eliminating the plateau observed in Fig 1, in which the reaction coefficient was assumed constant. Due to the incorporation of a forgetting-curve in the reaction coefficient we can also observe higher and earlier occurring epidemic peaks in the curves of new confirmed cases as compared to the constant reaction coefficient case (see Fig 1), especially for large μ0 values.

Fig 5 depicts the dynamics for a medium restoration coefficient ν = 0.5. Here we can observe– especially for a large maximum value μ0 of the reaction coefficient μ(t) = μ0 h(t)– that the restoration rate f(t, β) and the reaction rate g(t, I) differ more from each other, as compared to the case depicted in Fig 4. If we observe the curves for μ0 = 0.6 we see that initially, up to day 90, the reaction rate is clearly higher than the restoration rate, which again produces a sharp decrease in the transmission rate. Every time the restoration and reaction rate curves intersect, we can observe a local maximum/minimum in the transmission rate, which eventually increases and converges back to its natural value β0. The changing human behavior reflected in the restoration and reaction rates, observed for μ0 = 0.6, produces a late occurring but large epidemic peak, preceded by one small peak. This dynamics can also be explained by observing the shape of the effective reproduction number, which crosses the threshold three times (red curve).

Fig 6 shows the dynamics for a low constant restoration coefficient ν = 0.2. Here we observe initial oscillatory dynamics for the restoration and reaction rates, which explain the oscillatory dynamics in the transmission rate and also in the effective reproduction number, which crosses the threshold (see dotted line) several times, producing small oscillations in the new confirmed cases, followed by a large epidemic peak.

A non-autonomous βSEIR model and the example of COVID-19 in Chile and Italy: Time-varying reaction and restoration coefficients

Although it has not been proven that the virus SARS-CoV-2 is seasonal in nature, seasonality of transmission may be an important factor to consider since social behavior is environmentally driven [56]. The main hypothesis regarding seasonality indicates that the higher the temperature the fewer infections occur [57, 58].

Therefore, additional to a time-varying reaction rate, the restoration rate could depend, among other, on environmental factors such as for instance temperature. For example, in regions with predominantly Mediterranean climate, an annual oscillation of the daily mean temperature can be observed [59], as is the case in a large part of Italy or central Chile. The seasonality could also represent peoples’ behavior during winter months vs summer months, or during vacation periods and special holidays. These factors affect the tendency of the population to return to its natural transmission rate β0, while making it seasonal in nature.

Hence, to capture seasonal factors we consider an annual periodic restoration coefficient of the form (7) where ν0 is the average restoration coefficient value, σ the amplitude of the oscillation and tmax the moment in time, where the restoration coefficient is maximum. In general, social studies would be needed to determine the nature of the seasonality depending on peoples’ behavior and habits in specific cultural and geographical contexts, independent of disease presence. This could help determine the parameters present in the seasonal part of the restoration coefficient. For instance, tmax could happen during winter (where individuals meet inside), with a larger amplitude σ if a holiday coincides with cold weather and thus people tend to gather and restore normal conditions. Also, σ ∼ 0 implies ν(t)∼ν0 constant, and hence there would be no seasonal variation in the restoration coefficient.

Considering this non-constant seasonal behavior for the restoration coefficient ν(t) as given in Eq (7), as well as a time-varying reaction coefficient as described in Eq (6), we illustrate through Figs 7 and 8 the applicability of our model during the first months of the pandemic, fitting the model to COVID-19 data of new confirmed cases of Chile and Italy [60, 61]– minimizing in the least square sense the residuals between the outcome of our model and the data set– respectively from March 16th, 2020, to February 16th, 2021, and from February 24th, 2020, to October 31st, 2020. Table 3 shows the respective fixed and fitted parameter values used in each figure. Chile and Italy served as examples of countries located at different hemispheres that experienced COVID-19 in distinct ways, due to cultural differences, distinct levels of initial knowledge of the virus, seasons, etc., and this way permitted us to show that our model can capture different COVID-19 dynamics and partially describe them through reaction and restoration thrives of the population.

thumbnail
Fig 7. Simulation fitting COVID-19 data from Chile.

The first subplot depicts the restoration rate f(t, β) (dotted) and reaction rate g(t, I) (solid); the second plot the transmission rate β(t); and in the third plot the blue dots represent data of daily confirmed new cases of COVID-19 in Chile, from March 16th, 2020, to February 16th, 2021 [60]. The red curve represents the least square fit of the model to the data with parameter values as in Table 3 for the population of Chile, with N = 18 million individuals and initial conditions of the model (5) taken to be E0 = 20, I0 = 81, R0 = 0, S0 = NE0I0R0. The fit produces a root-mean-square error (RMSE) of 792.89.

https://doi.org/10.1371/journal.pone.0269843.g007

thumbnail
Fig 8. Simulation fitting COVID-19 data from Italy.

The first subplot depicts the restoration rate f(t, β) (dotted) and reaction rate g(t, I) (solid); the second plot the transmission rate β(t); and in the third plot the blue dots represent data of daily confirmed new cases of COVID-19 in Italy, from February 24th, 2020, to October 31st, 2020 [61]. The red curve represents the least square fit of the model to the data with parameter values as in Table 3 for the population of Italy, with N = 60.5 million individuals and initial conditions of the model (5) taken to be E0 = 81, I0 = 566, R0 = 0, S0 = NE0I0R0. The fit produces a root-mean-square error (RMSE) of 1083.49.

https://doi.org/10.1371/journal.pone.0269843.g008

thumbnail
Table 3. Parameters for the simulations in Figs 7 and 8 for the case of Chile and Italy respectively.

https://doi.org/10.1371/journal.pone.0269843.t003

We can observe how the restoration rate and reaction rate differ among countries: In the case of Chile (see Fig 7), they are very similar, expecting a population whose reaction to disease presence and urge to return to their normal behavior (restoration) is similar in nature. Nevertheless, a small difference in those rates produce a large impact on the transmission rate, such as the steep decrease observed initially.

On the other hand, from Fig 8 one can see that the reaction rate (solid curve) and the restoration rate (dotted curve) differ more than in the case of Chile, having initially a population, which is reacting to disease presence faster than their urge to restore normal conditions, for later reversing their behavior three times, producing during that time period a long lasting plateau. Finally, one can see that a larger restoration rate than reaction rate produces the appearance of a large second peak.

Discussion and conclusions

In references from 1973 [69] and 1989 [70], strategic models are defined as those that, despite not having high resolution concerning a specific reality, have the advantage of containing all the minimum aspects of the referenced system. The main contribution of this work is to achieve a dynamically richer low-cost model, that is, one that adds only one more differential law to the classical SEIR model (without introducing new compartments in the population).

There is evidence in the literature that populations change their behavior when facing dangerous diseases, i.e., reacting to these by managing to modify the transmission rate. Our proposed model, the βSEIR model, is a serious candidate to contain the minimum aspects for disease transmission of a high impact infectious-contagious disease in populations that, while living with the urge to restore normal conditions, react to reduce favorable conditions for disease transmission. Additionally, the latter occurs in a setting where vaccines are not available (and hence the model just considers one susceptible class), the disease is new to the population, and there are no multiple variants. Specifically, the novelty of the βSEIR model we present is, that it incorporates a variation in the transmission rate, which occurs proportional to the transmission rate itself but also proportional to a tension given by the difference (fg) between a) reaction rate (g) to disease presence that may include other behavioral factors, such as compliance with mitigation strategies, and b) restoration rate (f) that aims to restore a certain intrinsic value of disease transmission, due to for instance socio-environmental elements.

Our results show an important gain in dynamic possibilities even in the case where Eq (2) in the βSEIR model given in Eq (5) is autonomous, i.e. f and g do not depend explicitly on time. Indeed, we can see in Figs 13 the appearance of several epidemic peaks and initial oscillatory dynamics, explained by the tension between reaction and restoration thrives of a population. In particular, we observe that high restoration coefficients ν– affecting the restoration rate f and representing a higher urge of the population to return to normal conditions– induce temporary stabilization of the transmission rate after an initial drop, being the duration of this plateau larger, the larger the reaction coefficients μ (affecting the reaction rate g) are, i.e. the higher the self-protective reaction is to disease presence. When considering small ν values (a small urge to return to normality), we observe oscillations in the transmission rate, with higher amplitudes for higher μ values, i.e. amplitudes are higher if individuals’ reaction to disease presence is higher. These oscillations in the transmission rate generate oscillations in the effective reproduction number, which lay around one for a period of time proportional to the value of the reaction coefficient μ. Regarding the curve of new cases, we show that in general higher restoration coefficients ν produce unimodal behavior, whereas lower ν values generate the appearance of a finite number of peaks with decreasing peak size, in a way that for large reaction coefficients μ the timing between peaks is smaller. These results already define interesting future mathematical challenges.

In general, our results show how the transmission rate is impacted by the reaction rate g and the restoration rate f. In particular, we observe that a small difference between reaction and restoration rates may produce a large impact in the transmission rate. This highlights the importance of individual behavior in a pandemic setting, where even the behavior of a small number of individuals could change the dynamics of a disease drastically.

Curves of new confirmed cases with two and up to three waves in a one year time frame have been observed, differing among countries in time and size, with peaks and valleys of different heights, proper to a pandemic still under development in a population without protective immunity [3]. We notice that the βSEIR model can capture such patterns at the cost of varying the reaction coefficient μ(t) and for further dynamic richness varying the restoration coefficient ν(t). We can justify a time-varying reaction coefficient μ(t) that considers two aspects: first, it follows the shape of the Stringency-Index [45]—that records for each country government mitigation measures–, second, it reflects the reduced compliance with or adherence to mitigation strategies observed with time. The βSEIR model shows even richer dynamics when introducing such a time-varying reaction coefficient (see Figs 46).

Additionally, our model is capable of capturing the time series of new confirmed cases of Chile and Italy when including– additionally to including a time-varying reaction coefficient– a time-dependent seasonal variation in the restoration coefficient ν(t), reflecting distinct temporal and possibly behavioral characteristics of two countries located at different hemispheres during their first pandemic year. Our results are, at first glance, good indicators for the richness that a model of such low structural complexity, as the one proposed here, can provide.

Acknowledgments

We thank Victor Osores from Facultad de Ciencias Básicas at Universidad Católica del Maule, Talca, Chile, for his valuable feedback on the numerical results section of this paper.

References

  1. 1. World Health Organization—Coronavirus disease (COVID-19) Pandemic;. https://www.who.int/emergencies/diseases/novel-coronavirus-2019.
  2. 2. CSSE. Coronavirus COVID-19 Global Cases by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University (JHU);. https://systems.jhu.edu/tracking-covid-19/.
  3. 3. CSSE. Coronavirus COVID-19 Global Cases by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University (JHU);. https://coronavirus.jhu.edu/data/new-cases.
  4. 4. Kermack WO, McKendrick AG. A contribution to the mathematical theory of epidemics. Proceedings of the royal society of london Series A, Containing papers of a mathematical and physical character. 1927;115(772):700–721.
  5. 5. Anderson RM, May RM. Infectious diseases of humans: dynamics and control. Oxford university press; 1992.
  6. 6. Keeling M, Danon L. Mathematical modelling of infectious diseases. British Medical Bulletin. 2009;92(1). pmid:19855103
  7. 7. Kolokolnikov T, Iron D. Law of mass action and saturation in SIR model with application to Coronavirus modelling. Infectious Disease Modelling. 2021;6:91–97. pmid:33225113
  8. 8. Ke R, Sanche S, Romero-Severson E, Hengartner N. Fast spread of COVID-19 in Europe and the US suggests the necessity of early, strong and comprehensive interventions. medRxiv. 2020;.
  9. 9. Law KB, Peariasamy KM, Gill BS, et al. Tracking the early depleting transmission dynamics of COVID‐19 with a time‐varying SIR model. Scientific Reports. 2020;10:21721. pmid:33303925
  10. 10. Towers S, Patterson-Lomba O, Castillo-Chavez C. Temporal variations in the effective reproduction number of the 2014 West Africa Ebola outbreak. PLoS currents. 2014;6. pmid:25642357
  11. 11. Cowling BJ, Lau MS, Ho LM, Chuang SK, Tsang T, Liu SH, et al. The effective reproduction number of pandemic influenza: prospective estimation. Epidemiology (Cambridge, Mass). 2010;21(6):842. pmid:20805752
  12. 12. Tariq A, Undurraga EA, Laborde CC, Vogt-Geisse K, Luo R, Rothenberg R, et al. Transmission dynamics and control of COVID-19 in Chile, March-October, 2020. PLoS neglected tropical diseases. 2021;15(1):e0009070. pmid:33481804
  13. 13. Santamaría L, Hortal J. COVID-19 effective reproduction number dropped during Spain’s nationwide dropdown, then spiked at lower-incidence regions. Science of the Total Environment. 2021;751:142257. pmid:33181975
  14. 14. Hwang J, Park H, Jung J, Kim SH, Kim N. Basic and effective reproduction numbers of COVID-19 cases in South Korea excluding Sincheonji cases. Medrxiv. 2020;.
  15. 15. Tang B, Bragazzi NL, Li Q, Tang S, Xiao Y, Wu J. An updated estimation of the risk of transmission of the novel coronavirus (2019-nCov). Infectious disease modelling. 2020;5:248–255. pmid:32099934
  16. 16. Taghvaei A, Georgiou T, Norton L, Tannenbaum A. Fractional SIR epidemiological models. Scientific Reports. 2020;10:20882. pmid:33257790
  17. 17. Wang X, Gao D, Wang J. Influence of human behavior on cholera dynamics. Mathematical biosciences. 2015;267:41–52. pmid:26119824
  18. 18. Bizet NGC, Pena DKM. Time-dependent and time-independent SIR models applied to the COVID-19 outbreak in Argentina, Brazil, Colombia, Mexico and South Africa. arXiv preprint arXiv:200612479. 2020;.
  19. 19. Córdova-Lepe F, Gutiérrez-Aguilar R, Gutiérrez-Jara JP. Number of COVID-19 cases in Chile at 120 days with data at 21/03/2020 and threshold of daily effort to flatten the epi-curve. Medwave. 2020;20(2):e7861–e7861. pmid:32225133
  20. 20. Gutiérrez-Aguilar R, Córdova-Lepe F, Muñoz-Quezada MT, Gutiérrez-Jara JP. Model for a threshold of daily rate reduction of COVID-19 cases to avoid hospital collapse in Chile. Medwave. 2020;20(3):e7871–e7871. pmid:32469855
  21. 21. Zhao S, Lin Q, Ran J, Musa SS, Yang G, Wang W, et al. Preliminary estimation of the basic reproduction number of novel coronavirus (2019-nCoV) in China, from 2019 to 2020: A data-driven analysis in the early phase of the outbreak. International journal of infectious diseases. 2020;92:214–217. pmid:32007643
  22. 22. Alimohamadi Y, Taghdir M, Sepandi M. Estimate of the basic reproduction number for COVID-19: a systematic review and meta-analysis. Journal of Preventive Medicine and Public Health. 2020;53(3):151. pmid:32498136
  23. 23. Eikenberry SE, Mancuso M, Iboi E, Phan T, Eikenberry K, Kuang Y, et al. To mask or not to mask: Modeling the potential for face mask use by the general public to curtail the COVID-19 pandemic. Infectious Disease Modelling. 2020;5:293–308. pmid:32355904
  24. 24. Bisin A, Moro A. Learning epidemiology by doing: The empirical implications of a spatial-sir model with behavioral responses. National Bureau of Economic Research; 2020.
  25. 25. Willis MJ, Díaz VHG, Prado-Rubio OA, von Stosch M. Insights into the dynamics and control of COVID-19 infection rates. Chaos, Solitons & Fractals. 2020;138:109937. pmid:32834573
  26. 26. d’Onofrio A, Manfredi P. Information-related changes in contact patterns may trigger oscillations in the endemic prevalence of infectious diseases. Journal of Theoretical Biology. 2009;256(3):473–478. pmid:18992258
  27. 27. Pedro SA, Ndjomatchoua FT, Jentsch P, Tcheunche JM, Anand M, Bauch CT. Conditions for a second wave of COVID-19 due to interactions between disease dynamics and social processes. medRxiv. 2020;.
  28. 28. Bauch CT. Imitation dynamics predict vaccinating behaviour. Proceedings of the Royal Society B: Biological Sciences. 2005;272(1573):1669–1675. pmid:16087421
  29. 29. Pharaon J, Bauch C. The Influence of Social Behavior on Competition Between Virulent Pathogen Strains. bioRxiv. 2018; p. 293936.
  30. 30. Inglesby T, Nuzzo J, O’Toole T, Henderson DA. Disease Mitigation Measures in the Control of Pandemic Influenza. Biosecurity and Bioterrorism. 2006;4(4):366–375. pmid:17238820
  31. 31. Milne G, Xie S. The Effectiveness of Social Distancing in Mitigating COVID-19 Spread: a modelling analysis. medRxiv. 2020; p. 2020.03.20.20040055.
  32. 32. Harper CA, Satchell D Liam P & Fido, D LR. Functional Fear Predicts Public Health Compliance in the COVID-19 Pandemic. Int J Ment Health Addiction. 2020;.
  33. 33. Malecki KMC, Malecki JA, Keating JA, Safdar N. Crisis Communication and Public Perception of COVID-19 Risk in the Era of Social Media. Clinical Infectious Diseases. 2021;72(4):697–702. pmid:32544242
  34. 34. Michaelis L, Menten ML. Die kinetik der invertinwirkung. Biochem z. 1913;49(333-369):352.
  35. 35. Vasconcelos GL, Brum AA, Almeida FAG, Macêdo AMS, Duarte-Filho GC, Ospina R. Standard and anomalous second waves in the COVID-19 pandemic. medRxiv. 2021; p. 2021.01.31.21250867.
  36. 36. Gronwall TH. Note on the derivatives with respect to a parameter of the solutions of a system of differential equations. Annals of Mathematics. 1919; p. 292–296.
  37. 37. Brauer F, Castillo-Chavez C. Mathematical models in population biology and epidemiology. vol. 2. Springer; 2012.
  38. 38. Heesterbeek J, Dietz K. The concept of Ro in epidemic theory. Statistica neerlandica. 1996;50(1):89–110.
  39. 39. Van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical biosciences. 2002;180(1-2):29–48. pmid:12387915
  40. 40. Kojima N, Klausner JD. Protective immunity after recovery from SARS-CoV-2 infection. The Lancet infectious diseases. 2022;22(1):12–14. pmid:34762853
  41. 41. Cabrera M, Córdova-Lepe F, Gutiérrez-Jara JP, Vogt-Geisse K. An SIR-type epidemiological model that integrates social distancing as a dynamic law based on point prevalence and socio-behavioral factors. Scientific Reports. 2021;11(1):1–16. pmid:33986347
  42. 42. Python Software Foundation;. https://www.python.org/.
  43. 43. Nishiura H, Chowell G. The effective reproduction number as a prelude to statistical estimation of time-dependent epidemic trends. In: Mathematical and statistical estimation approaches in epidemiology. Springer; 2009. p. 103–121.
  44. 44. Gumel AB, Iboi EA, Ngonghala CN, Elbasha EH. A primer on using mathematics to understand COVID-19 dynamics: Modeling, analysis and simulations. Infectious Disease Modelling. 2021;6:148–168. pmid:33474518
  45. 45. Hale T, Angrist N, Goldszmidt R, Kira B, Petherick A, Phillips T, et al. A global panel database of pandemic policies (Oxford COVID-19 Government Response Tracker). Nature Human Behaviour. 2021;5(4):529–538. pmid:33686204
  46. 46. Our World in Data, Oxford Martin School, University of Oxford;. https://ourworldindata.org/covid-stringency-index.
  47. 47. COVID-19 en América Latina y el Caribe: Panorama de las respuestas de los gobiernos a la crisis;.
  48. 48. Loayza NV. Costs and Trade-Offs in the Fight Against the COVID-19 Pandemic. Research and Policy Briefs From the World Bank Malaysia Hub;35.
  49. 49. Li Y, Li M, Rice M, Zhang H, Sha D, Li M, et al. The Impact of Policy Measures on Human Mobility, COVID-19 Cases, and Mortality in the US: A Spatiotemporal Perspective. Int J Environ Res Public Health. 2020;18:996.
  50. 50. Blanco F, Emrullahu D, Raimundo S. Do coronavirus containment measures work? Policy Research Working Paper. 2020;9490:996.
  51. 51. Goldstein P, Levy Yeyati L Eduardo Yeyati andSartorio. Lockdown fatigue: The diminishing effects of quarantines on the spread of COVID-19. Covid Economics. 2021;67:1–23.
  52. 52. Rozanova L, Temerev A, Flahault A. Comparing the Scope and Efficacy of COVID-19 Response Strategies in 16 Countries: An Overview. Int J Environ Res Public Health. 2020;17:9421. pmid:33339119
  53. 53. Wright L, Fancourt D. Do predictors of adherence to pandemic guidelines change over time? A panel study of 21,000 UK adults during the COVID-19 pandemic. medRxiv. 2020;20228403.
  54. 54. Nese M, Riboli G, Brighetti G, Sassi V, Camela E, Caselli G, et al. Delay discounting of compliance with containment measures during the COVID-19 outbreak: a survey of the Italian population. J Public Health. 2020;. pmid:32837837
  55. 55. Krpan D, MAKKI F, SALEH N, BRINK H S ND KLAUZNICER. When behavioural science can make a difference in times of COVID-19. Behavioural Public Policy. 2021;5(2):153–179.
  56. 56. Carlson CJ, Gomez ACR, Bansal S, Ryan SJ. Misconceptions about weather and seasonality must not misguide COVID-19 response. Nat Commun. 2020;11:4312. pmid:32855406
  57. 57. Malki Z, Atlam ES, Hassanien AE, Dagnew G, Elhosseini MA, Gad I. Association between weather data and COVID-19 pandemic predicting mortality rate: Machine learning approaches. Chaos, Solitons & Fractals. 2020;138:110137. pmid:32834583
  58. 58. Merow C, Urban MC. Seasonality and uncertainty in global COVID-19 growth rates. Proceedings of the National Academy of Sciences. 2020;117(44):27456–27464. pmid:33051302
  59. 59. Climate Data Worldwide;. https://en.climate-data.org/.
  60. 60. Ministerio de Ciencia, Tecnología, Conocimiento e Innovación, Base de Datos COVID-19;. https://www.minciencia.gob.cl/covid19/.
  61. 61. Statista. Daily new coronavirus (COVID-19) cases in Italy since February 2020, by date of report;. https://www.statista.com/statistics/1101690/coronavirus-new-cases-development-italy/.
  62. 62. Transmission of SARS-CoV-2: implications for infection prevention precautions;. https://www.who.int/news-room/commentaries/detail/transmission-of-sars-cov-2-implications-for-infection-prevention-precautions.
  63. 63. Rai B, Shukla A, Dwivedi LK. Incubation period for COVID-19: a systematic review and meta-analysis. Journal of Public Health. 2021; p. 1–8. pmid:33643779
  64. 64. Zhou F, Yu T, Du R, Fan G, Liu Y, Liu Z, et al. Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study. The lancet. 2020;395(10229):1054–1062. pmid:32171076
  65. 65. Tang B, Wang X, Li Q, Bragazzi NL, Tang S, Xiao Y, et al. Estimation of the transmission risk of the 2019-nCoV and its implication for public health interventions. Journal of clinical medicine. 2020;9(2):462. pmid:32046137
  66. 66. Setti L, Passarini F, De Gennaro G, Barbieri P, Perrone MG, Piazzalunga A, et al. The potential role of particulate matter in the spreading of COVID-19 in northern Italy: first evidence-based research hypotheses. MedRxiv. 2020;.
  67. 67. Wangping J, Ke H, Yang S, Wenzhe C, Shengshu W, Shanshan Y, et al. Extended SIR prediction of the epidemics trend of COVID-19 in Italy and compared with Hunan, China. Frontiers in medicine. 2020;7:169. pmid:32435645
  68. 68. Ke R, Romero-Severson E, Sanche S, Hengartner N. Estimating the reproductive number R0 of SARS-CoV-2 in the United States and eight European countries and implications for vaccination. Journal of theoretical biology. 2021;517:110621. pmid:33587929
  69. 69. May RM. Stability and complexity in model ecosystems. vol. 6. Princeton University Press; 1973.
  70. 70. Yodzis P. Introduction to theoretical ecology. HarperCollins College Division; 1989.