Data

Six sources of data are analyzed in this study. A brief description of the data sets is provided as follows:

1. (i) Case incidence data

Case incidence data based on the dates of symptom onset is used to generate the epidemic curve, scaled incidence curve (see S1 Text), short-term forecasts of the COVID-19 pandemic, assess the spatial heterogeneity in Colombia and estimate the national and regional reproduction numbers. Specifically, we utilize publicly available time series of laboratory-confirmed cases based on the dates of symptom onset for the national and regional COVID-19 case incidence. These data are retrieved from the Colombian Ministry of Health as of October 31, 2021 [23]. Colombian departments and their four main districts; Barranquilla, Bogotá, Cartagena, and Santa Marta can be grouped into five regions as follows: the Amazon, the Andean, the Orinoquía, the Caribbean and the Pacific [41]. To evaluate the spatial heterogeneity and estimate the regional reproduction numbers, for simplicity we assume a regional level distribution as follows (S1 Fig):

1. Amazon region: Amazonas, Caquetá, Guainía, Guaviare, Putumayo and Vaupés.
2. Andean region: Antioquia, Boyacá, Bogotá, Caldas, Cauca, Cundinamarca, Huila, Norte De Santander, Quindío, Risaralda, Santander, and Tolima.
3. Orinoquía region: Arauca, Casanare, Meta, and Vichada
4. Caribbean region: Atlántico, Bolívar, Cesar, Córdoba, La Guajira, Magdalena, San Andrés, Sucre, Barranquilla, Cartagena, and Santa Marta
5. Pacific region: Chocó, Nariño, and Valle Del Cauca

Based on the geographical distribution of the departments, Nariño, Valle Del Cauca, Guainía, and Cauca appear in two regions each. However, we include them in one region as stated above, to avoid the duplication of data for our analysis.

1. (ii) Mortality data

COVID-19 mortality data based on the date of death is retrieved from the Colombian Ministry of Health as of October 31, 2021 [23]. Mortality data is used to generate the national short-term forecast of the COVID-19 pandemic and estimate the national reproduction number.

1. (iii) Genomic data

To estimate the reproduction number from the genomic data, from the 1855 samples available as of June 30, 2021, 136 SARS-CoV-2 genome samples for Colombia are obtained from the “Global Initiative on Sharing Avian Influenza Data” (GISAID) repository [42] between February 27 and April 5, 2020 to estimate the reproduction number.

1. (iv) Mobility trends data

Two sources of mobility data are analyzed at the national level from March 6, 2020, to October 31, 2021. These data were retrieved as of October 31, 2021.

Apple’s mobility data is a set of data that summarizes the number of requests made to Apple Maps for directions. The data are categorized by countries/regions, sub-regions, and cities. Apple has released the data for the three modes of human mobility: driving, walking and public transit. These data are published publicly by Apple’s mobility trends reports [43]. These aggregated and anonymized data are updated daily and includes the relative volume of directions requests per country compared to the baseline volume on January 13, 2020. To emphasize on the change rather than the gross number of requests, Apple has normalized all the data. The number of direction requests for each country or city and the type of direction request (driving, walking and public transit) has been set to 100 for January 13, 2020. All subsequent data are relative to the starting point.

Google’s mobility data [44] is retrieved from users who have turned on the location history in their Google account. Google’s mobility data shows the changes in mobility in each geographic and administrative region for visits to different places, such as grocery stores, parks, and recreation spots. The baseline day represents the normal value on the day of the week. The baseline day is the median value from the 5-week period from January 3, 2020, to February 6, 2020.

1. (v) Twitter data

We retrieved data from the publicly available Twitter data set of the COVID-19 chatter from March 12, 2020 to October 31, 2021 [45]. A detailed description of Twitter data acquisition and processing is provided in a prior study [45].

Modeling framework for forecast generation

We utilize and compare three dynamic phenomenological growth models to generate short-term (i.e., 30-day ahead) forecasts for Colombia. These models have been applied to various infectious diseases including SARS, foot and mouth disease, Ebola [46–48] and the current COVID-19 outbreak [49–51]. The phenomenological growth models applied in this study include the two models defined by one scalar differential equation namely the generalized logistic growth model (GLM) [47] and the Richards growth model [52]. Additionally we apply the sub-epidemic wave model [46] which captures diverse epidemic trajectories such as the multimodal outbreaks. The forecasts obtained from these dynamic growth models can assess the potential scope of the pandemic in near real-time, provide insights on the contribution of disease transmission pathways, predict the impact of control interventions, and evaluate optimal resource allocation to inform public health policies. A detailed description of these models is provided in S2 Text.

Model calibration and forecasting approach

We utilize the national and regional level case incidence data retrieved from the Colombian Ministry of Health based on the dates of symptom onset as of October 31, 2021. We also utilize the national level mortality data based on the date of death as of October 31, 2021. Each forecast is fitted to the daily case counts based on the dates of symptom onset and daily death counts based on the date of death between July 4, 2021, and October 1, 2021 (90 days calibration period) to conduct a 30-day ahead short-term forecast for each model. The data from October 2, 2020, to October 31, 2021, is utilized to assess the performance of our 30-day ahead short-term forecasts. Details of the model calibration and forecasting approach are provided in S3 Text.

Performance metrics

The following five performance metrics are used to assess the quality of our model fit and the 30-day ahead short-term forecasts: the mean absolute error (MAE) [53]; the root mean squared error (RMSE) [54]; the coverage of the 95% prediction intervals (95% PI) [54]; the mean interval score (MIS) [54] and the weighted interval score (WIS) [55]; for each of the three models: GLM; Richards model; and the sub-epidemic wave model. A detailed description of the performance metrics is provided in S4 Text.

Reproduction number

In order to understand the transmission dynamics of COVID-19, we estimate the effective reproduction number, R_t, for the early ascending phase from February 27, 2020, to March 27, 2020, and the effective reproduction number R_t throughout the pandemic for the national and regional COVID-19 epidemic curves. The effective reproduction number, R_t, is the key parameter that characterizes the average number of secondary cases generated by a primary case at calendar time (t) during the outbreak. This quantity is crucial for identifying the magnitude of public health interventions required to contain an epidemic [56–58]. The estimates of R_t indicate whether widespread disease transmission continues (R_t>1) or disease transmission declines (R_t<1). Therefore, to contain an outbreak, it is vital to maintain R_t<1. We utilize the case incidence, mortality, and the genomic data to estimate the reproduction numbers using three rigorous methods to compare the estimates of R_t obtained from each method and data source.

Effective reproduction number R_t, using the generalized growth model (GGM)

We estimate the national and regional reproduction numbers by calibrating GGM to the early growth phase of the pandemic [59], from February 27, 2020 to March 27, 2020. The description of the GGM is provided in S2 Text. The generation interval of SARS-CoV-2 is modeled with the assumption of a gamma distribution with a mean of 5.2 days and a standard deviation of 1.72 days [60]. We first characterize the daily incidence of local cases using the GGM. The progression of local incidence cases by dates of symptom onset, I_i, is simulated using the calibrated GGM model and account for the daily series of imported cases by dates of symptom onset, J_i, into the renewal equation to estimate the effective reproduction, as

The factor J_i represents the imported cases at time t_i, I_i denotes the local case incidence at calendar time t_i and ρ_j represents the discretized probability distribution of the generation interval. The factor 0≤α≤1 represents the relative contribution of imported cases to secondary disease transmission. We perform sensitivity analysis by setting α = 0.15, α = 1.00 and α = 0.00 to assess the relative contribution of the imported cases to the secondary disease transmission [40]. The numerator represents the total new cases I_i, and the denominator represents the total number of primary cases that contributes to the generation the new cases I_i (i.e., as secondary cases) at time t_i. Hence, R_t, represents the average number of secondary cases generated by a single case at calendar time t. The uncertainty bounds around the curve of are derived directly from the 300 bootstrap samples generated from the GGM model with estimated parameters () where i = 1,2, …S. We assume a negative binomial error structure where the variance is assumed to be three times the mean based on the noise of the data [61]. This method is utilized to derive the early estimates of reproduction number and has been employed in several prior studies [50,61,62].

Since the national and regional epidemic curves have a distinct date of onset according to the seeding of the first local case, therefore, R_t for the national and regional level is estimated based on the date of onset of the first local case. For the national and regional epidemic curves, we estimate R_t for the first 30 days (Table 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Dates for R_t estimation for the national and regional epidemic curves for the first 30 epidemic days.

https://doi.org/10.1371/journal.pntd.0010228.t001

Effective reproduction number R_t, using the Cori et al. method

We estimate the regional and national effective reproduction numbers (also known as the time dependent instantaneous reproduction numbers) using the local and imported case incidence data based on the date of symptoms onset as of October 31, 2021. The total number of incident cases arising at time step t, I_t, is the sum of the number of local incident cases () and imported incident cases (). Hence,

We assume that if the imported cases exist, they can be distinguished from the local cases through epidemiological investigations so that and can be observed at each time step [63,64].

The time dependent (instantaneous) R_t is defined as the ratio of the number of new locally infected cases generated at calendar time t (), and the total infectiousness across all infected individuals at time t given by [65,66]. Hence R_t can be written as

In this equation, the factor w_s represents the infectivity function, which is the infectivity profile of the infected individual. The infectivity function is dependent on the time since infection (s), but is independent of the calendar time (t) [67,68].

The term describes the sum of infection incidence up to time step t − 1, weighted by the values of infectivity function w_s. The distribution of the generation time can be applied to approximate w_s, however, since the time of infection is a rarely observed event, measuring the distribution of generation time becomes difficult [65]. Therefore, the time of symptom onset is usually used to estimate the distribution of the serial interval (SI), which is defined as the time interval between the dates of symptom onset among two successive cases in a disease transmission chain [69].

The infectiousness of a case is a function of the time since infection, which is proportional to w_s if the timing of infection in the primary case is set as time zero of w_s and we assume that the generation interval equals the serial interval (SI). The SI is assumed to follow a gamma distribution with a mean of 5.2 days and a standard deviation of 1.72 days [60]. We also estimate R_t for the mortality curve by the date of death as of October 31, 2021. The analytical estimates of R_t are obtained within a Bayesian framework using the Cori et al. method, employing the EpiEstim package in R version 4.0.3. [69]. The values for R_t are estimated using 1-weekly sliding window. We report the median and 95% credible interval (CrI).

Reproduction number, R, from the genomic analysis

To estimate the reproduction number for SARS-CoV-2 between February 27, 2020 and April 5, 2020 from the genomic data, 136 SARS-CoV-2 genome samples from Colombia and their sampling times are obtained from the GISAID repository [42,70]. Short sequences and sequences with significant number of gaps and non-identified nucleotides are removed. For clustering, they are complemented by the sequences from other geographical regions across the globe, downsampled to n = 1470 representative sequences using a country specific sub-sampling algorithm of Nextstrain. Briefly, the algorithm first chooses the sequences from the focal country, and then complements them with the sequences obtained from other geographical regions, prioritizing sequences from the geographical neighbors and the sequences that are genetically close to the focal sequences. This decreases the size of the genomic data set to allow for feasible phylogenetic inference. We use the sequence subsample from Nextstrain (www.nextstrain.org) global analysis as of July 4, 2021. These sequences are aligned to the reference genome obtained from the literature [71] using Muscle [72] and trimmed using Molecular Evolutionary Genetics Analysis version: 7 (Mega 7) [73].

The largest Colombian cluster that possibly corresponds to within-country transmissions has been identified using hierarchical clustering of sequences. Phylodynamics analysis of that cluster is carried out using Bayesian Evolutionary Analysis by Sampling Trees (BEAST) [74] version 2. We use a strict molecular clock and a tree prior with exponential growth coalescent. Markov Chain Monte Carlo (MCMC) sampling is run for 10,000,000 iterations, and the parameters are sampled every 1000 iterations. The results are acceptable if the effective sample sizes are above 200 for all parameters. The exponential growth rate, f estimated by BEAST is used to calculate the basic reproduction number, R₀. We utilize the standard assumption that SARS-CoV-2 generation intervals (times between infection and onward transmission) are gamma-distributed [75]. In this case, R₀ can be estimated as where μ and σ are the mean and standard deviation of the gamma distribution, respectively. The values of μ = 5.2 days and σ = 1.72 days are taken from a prior study [60]. The formula for R₀ is a strictly monotone transformation of , therefore it also transforms the 95% highest posterior density (HPD) intervals for f into those for R₀.

Spatial analysis

We quantify and statistically analyze the shapes of the incidence rate curves. Shape analysis of the national and regional incidence rate curves is performed following the analytical methodology as described in previous literature [76] including steps to pre-process the daily cumulative COVID-19 case data based on the dates of symptom onset. Then, we analyze the shapes of these rate curves to compare, cluster and summarize the incidence rates. The pre-processing steps applied to the COVID-19 daily count data are as follows:

1. Smoothing: Cumulative case curves are smoothed over 10 days using the smooth function in the MATLAB. The smooth function returns a moving average of the elements of a vector using a fixed window length which is determined heuristically, and in the case of spatial analysis, it is 10 days. The window slides down the length of the vector, computing an average over the elements within each window.
2. Time differencing: If f_i(t) denotes the given cumulative number of confirmed cases for department i on day t, then the growth rate per day at time t is given by g_i(t) = f_i(t)−f_i(t−1).
3. Rescaling: Each curve is rescaled by dividing each g_i by the total number of confirmed cases.
4. Smoothing: We then smooth the normalized curves again over a 5-day span using the smooth function in the MATLAB.

To identify the clusters by comparing the curves, we used a simple metric. For any two rate curves, h_i and h_j, we compute the norm ||h_i −h_j||, where the double bars denote the L² norm of the difference function, that is,

To perform clustering of thirty six curves into smaller groups, we apply the dendrogram function in the MATLAB using the “ward” linkage as explained in reference [76]. The number of clusters is determined empirically based on the display of overall clustering results. After clustering the departments and administrative divisions into different groups, we derive the average curve for each cluster after using a time wrapping algorithm as performed in previous studies [76,77]. All analysis is conducted in MATLAB R2020b.

Twitter data analysis

To observe any relationship between the COVID-19 cases based on the dates of symptom onset and the frequency of tweets indicating stay-at-home orders we used a public dataset of 9,233,913,414 unique tweets of the COVID-19 chatter [45]. The frequency of tweets indicating stay-at-home order is used to infer people’s compliance with the stay-at-home-orders implemented to avoid spread of the virus by practicing social distancing. Tweets indicate the number of the people in favor of lockdowns and show how these numbers have varied throughout the pandemic. We removed all retweets and tweets not in the Spanish language to obtain the plotted data. We also filtered the tweets using the following hashtags: #quedateencasa (stay-at-home), and #trabajardesdecasa (work-from-home), which are two of the most used hashtags when users refer to the COVID-19 pandemic and their engagement with health measures. Lastly, our analysis was restricted to the tweets that originated from Colombia. A set of 493,276 unique tweets were gathered between March 12, 2020, and October 31, 2021. The usage of Twitter in Colombia is less compared to other Latin American countries like Mexico and Chile [78], hence we observe a lower proportion of filtered tweets. We overlay the curve of filtered tweets and the epidemic curve in Colombia to observe any relation between the shape of the epidemic trajectory and the shape of curve for the frequency of filtered tweets during the established time-period. The correlation coefficient between the daily cases based on the dates of symptom onset and frequency of filtered tweets is also estimated using the "corrcoef” function in MATLAB version R2020b. This coefficient explains the linear dependency of the daily cases based on the dates of symptom onset and frequency of filtered tweets.

Mobility data analysis

We utilize the R code developed by Healy [79] to analyze the time series data for Colombia from March 6, 2020 to October 31, 2021, for two modes of mobility: driving and walking derived from Apple’s mobility reports. We analyze the mobility trends to discover any commonality in pattern of the mobility curves with the epidemic curve of COVID-19. The time series for mobility requests is decomposed into trends, weekly and remainder components. The trend is obtained by applying locally weighted regression fitted to the data and the remainder is any residual leftover on any given day after accounting for the underlying trend and normal daily fluctuations. The analysis is performed using ‘covmobility’ package in R software, version 4.0.3.

We analyze the Google mobility trends from March 6, 2020, to October 31, 2021 for visits to the grocery stores and pharmacy, parks, transit stations, workplaces, and residential areas compared to the baseline which is the median value, for the corresponding day of the week, during the five weeks from January 3, 2020 to February 6, 2020. For each region category, the baseline is not a single value; it comprises of seven individual values. The same number of visitors on two different days of the week results in different percentage changes. The analysis is performed using ‘tidycovid19’ package and R software, version 4.0.3.

Model calibration and forecasting performance

We compare the results from model calibration and 30-day ahead forecasting across three models: GLM, Richards and the sub-epidemic wave model. Model calibration, using 90 days of the epidemic data, between July 4, 2021 and October 1, 2021 shows that at the regional and national levels, the sub-epidemic model outperformed the GLM and Richards growth model in terms of the all five performance metrics (RMSE, MAE, MIS, WIS and 95% PI coverage) (Table 2) using the case incidence data. Therefore, the sub-epidemic wave model can be declared the most accurate model for the calibration period. The model fits exhibited sub-exponential growth dynamics [48,59] for three models (p~0.6–0.8) at the national and regional levels. The calibration performances for each region and the national data are listed in Table 2. Calibrating the three models to the mortality data also shows that the sub-epidemic model outperforms the other two models (S1 Table) and exhibits sub-exponential growth dynamics.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Comparison of model performance metrics by calibrating the GLM, Richards and the sub-epidemic model for 90 epidemic days (July 4, 2021 to October 1, 2021) at the national and regional level.

Higher 95% PI coverage and lower RMSE, MAE, WIS and MIS represent better performance. Best performing model is given in bold with the superscript “a”.

https://doi.org/10.1371/journal.pntd.0010228.t002

In terms of the forecasting performance, again the sub-epidemic wave model performed better than the GLM and Richards model for the national and regional case incidence data (Table 3). The sub-epidemic model also outperformed the GLM and Richards model in the forecasting phase of the mortality curve (S2 Table). Hence, the sub-epidemic model can be considered the most accurate model to forecast the epidemic trajectory.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Comparison of 30-day ahead forecasting performance (October 2, 2021 to October 31, 2021) by calibrating the GLM, Richards and the sub-epidemic model for 90 epidemic days (July 4, 2021 to October 1, 2021) at the national and regional level.

Higher 95% PI coverage and lower RMSE, MAE, WIS and MIS represent better performance. Best performing model is given in bold with the superscript "a”.

https://doi.org/10.1371/journal.pntd.0010228.t003

30-day ahead forecasts

Calibrating our models from July 04, 2021, to October 1, 2021 and generating the 30-day ahead forecasts from October 2, 2021 to October 31, 2021 utilizing the GLM and Richards growth model indicates a downwards trend for the national and regional case incidence data (Figs 3 and 4). On the other hand, the sub-epidemic wave model captures the multiple sub-epidemics comprising the coarse aggregated COVID-19 epidemic wave of Colombia. The sub-epidemic model predicts a downward trend for the Amazon and Orinoquía region, consistent with the findings of the GLM and Richards model. Whereas, for the national, Andean, Caribbean and the Pacific region, the sub-epidemic model predicts a stable case incidence pattern (Fig 5). According to the GLM and Richards model, the COVID-19 pandemic in Colombia would decline to zero during the month of October whereas the sub-epidemic wave model predicts an accumulation of 24,525 (95% PI: 13,677, 44,515) cases at the national level for the month of October 2021.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. 30-days ahead forecasts of the national and regional COVID-19 epidemic curves in Colombia by calibrating the GLM model from July 4, 2021 to October 1, 2021.

The blue circles correspond to the data points; the solid red line indicates the best model fit, and the red dashed lines represent the 95% prediction interval. The vertical black dashed line represents the time of the start of the forecast period.

https://doi.org/10.1371/journal.pntd.0010228.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. 30-days ahead forecast of the national and regional COVID-19 epidemic curves in Colombia by calibrating the Richards model from July 4, 2021 to October 1, 2021.

The blue circles correspond to the data points; the solid red line indicates the best model fit, and the red dashed lines represent the 95% prediction interval. The vertical black dashed line represents the time of the start of the forecast period.

https://doi.org/10.1371/journal.pntd.0010228.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. 30-days ahead forecast of the national and regional COVID-19 epidemic curves in Colombia by calibrating the sub-epidemic wave model from July 4, 2021 to October 1, 2021.

The blue circles correspond to the data points; the solid red line indicates the best model fit, and the red dashed lines represent the 95% prediction interval. The vertical black dashed line represents the time of the start of the forecast period.

https://doi.org/10.1371/journal.pntd.0010228.g005

At the regional level, the Richards and the GLM model predict zero cases for the month of October (Figs 3 and 4). However, the sub-epidemic model predicts 59 (95% PI: 0, 507) cases for the Amazon region in the month of October. The sub-epidemic model also predicts 12,186 (95% PI: 7107, 19,235) cases for the Caribbean region, 18,165 (95% PI: 9205, 31,205) cases for the Andean region, 5039 (95% PI: 2517, 8725) cases for the Pacific region and 510 (95% PI: 87, 1456) cases for the Orinoquía region (Fig 5). The 30-day ahead forecast of the national mortality data generated by the GLM and Richards model indicate a decline in the number of deaths whereas the sub-epidemic wave model indicates an upward trend in the mortality curve with ~2893 (95% PI: 1860, 5325) deaths that can accumulate in the month of October 2021 (S3 Fig).

Reproduction number

Estimate of the effective reproduction number, R_t from case incidence data.

The reproduction number for the early ascending growth phase of the epidemic from the case incidence data (February 27, 2020 to March 27, 2020) using GGM was estimated at R_t~1.30 (95% CI:1.20, 1.50) at α = 0.15 for the national data. The growth rate parameter, r, was estimated at 1.40 (95% CI: 0.91, 2.0) and the deceleration of growth parameter, p, was estimated at 0.64 (95% CI: 0.56, 0.71), indicating early sub-exponential growth dynamics of the COVID-19 pandemic in Colombia (Fig 6). Simultaneously, the estimates of R_t for all the regions remained consistently above 1.20 (between ~1.20–2.22) for the early ascending phase of the pandemic with the Amazon region showing the highest estimate of reproduction number, at R_t ~2.2 followed by the Orinoquía region with R_t~1.8. The estimates of R_t for the Andean, Pacific and Caribbean regions remained between R_t~1.2–1.4. All regions except the Amazon region depict sub-exponential growth dynamics for the COVID-19 pandemic in Colombia with the deceleration of the growth parameter, p, estimated between 0.54–0.86. The Amazon region shows almost exponential growth dynamics with the deceleration of growth parameter, p ~0.95 (95% CI: 0.74, 1.00). In contrast, the Andean and the Caribbean regions show an almost linear pattern of the epidemic trajectory with the deceleration of growth parameter estimated at, p ~0.58 (95% CI: 0.48, 0.69) and p ~0.59 (95% CI: 0.34, 0.87) respectively. The sensitivity analysis shows that the estimates of reproduction numbers do not vary significantly at α = 0.15 and α = 1.00 (Tables 4 and 5). Moreover, setting α = 0.00, thereby assuming that there is no contribution of imported cases to the disease transmission process or generation of secondary cases also does not substantially change the estimates of reproduction number.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6.

Upper panel: Reproduction number for Colombia with 95% CI estimated using the GGM model. The estimated reproduction number of the COVID-19 epidemic in Colombia as of March 27, 2020, is 1.30 (95% CI: 1.20, 1.50). The growth rate parameter, r, is estimated at 1.40 (95%CI: 0.91, 2.0) and the deceleration of the growth parameter, p, is estimated at 0.64 (95%CI: 0.56, 0.71) at α = 0.15. Lower panel: The lower panel shows the GGM fit to the case incidence data for the first 30 days from February 27, 2020 to March 27, 2020. The blue circles correspond to the data points; the solid red line indicates the best model fit, and the red dashed lines represent the 95% confidence interval. The cyan lines are the model fits obtained via bootstrapping.

https://doi.org/10.1371/journal.pntd.0010228.g006

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Estimates of reproduction number (R), growth rate parameter (r) and deceleration of the growth parameter (p) obtained from the renewal equation method utilizing the GGM for the early ascending phase of the epidemic (30 days) at the national and regional level at α = 0.15.

https://doi.org/10.1371/journal.pntd.0010228.t004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Estimates of reproduction number (R), growth rate parameter (r) and deceleration of the growth parameter (p) obtained from the renewal equation method utilizing the GGM for the early ascending phase of the epidemic (30 days) at the national and regional level at α = 1.00.

https://doi.org/10.1371/journal.pntd.0010228.t005

Estimate of effective (instantaneous) reproduction number, R_t, from the Cori et al. method.

The early estimate of effective (instantaneous) reproduction number at the national level for the first 30 epidemic days exhibits a steep decline. The reproduction number was estimated at, R_t~2.85 (95% CrI: 1.91, 4.09) on March 6, 2020, which declined to R_t~0.81 (95% CrI: 0.76, 0.85) by March 27, 2020 (Figs 7 and 8). At the national level the effective reproduction number has remained consistently above R_t~1.00 between April and July 3, 2020, after which the reproduction number declined to less than 1.00 (R_t~0.9) until September 22, 2020. Since the end of September 2020, R_t has fluctuated around ~1.00 (R_t~ 0.8–1.4) until the end of April 2021. From May till mid-June 2021, the estimates of R_t remained above ~1.0 indicating continued disease transmission. Since then, the estimates of R_t have mostly remained below ~1.0 with the most recent estimate of R_t~0.84 (95% CrI: 0.82, 0.86) as of October 31, 2021 (Fig 8).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Effective (instantaneous) reproduction number with 95% CrI estimated using the Cori et al. method and the estimate of reproduction number with 95% CI utilizing the GGM for simulating the parameters employed in the renewal equation method for the COVID-19 epidemic in Colombia as of March 27, 2020 (first 30 days).

The red solid line is the mean reproduction number estimated from the renewal equation method and the red dashed lines are the 95% CI coverage around the mean reproduction number. The green solid line is the median effective reproduction number estimated from the Cori et al. method and the green dashed lines are the 95% CrI coverage around the median instantaneous reproduction.

https://doi.org/10.1371/journal.pntd.0010228.g007

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8.

Upper panel: Effective reproduction number estimated from the Cori et al. method with 95% credible intervals for the COVID-19 pandemic in Colombia as of October 31, 2021. The red solid line represents the mean reproduction number for Colombia and the red shaded area represents the 95% credible interval around it. The blue solid line represents the mean reproduction number in the Caribbean region and the blue shaded area represents the 95% credible interval around it. The yellow solid line represents the mean reproduction number in the Pacific region and the yellow shaded area represents the 95% credible interval around it. The brown solid line represents the mean reproduction number in the Orinoquía region and the brown shaded area represents the 95% credible interval around it. The pink solid line represents the mean reproduction number in the Andean region and the pink shaded area represents the 95% credible interval around it. The green solid line represents the mean reproduction number in the Amazon and the green shaded area represents the 95% credible interval around it. Lower panel: The percentage of the COVID-19 cases in Colombia, nationally and regionally as of October 31, 2021. The red solid line represents the percentage of cases in Colombia, the blue solid line represents the percentage of cases in the Caribbean region, the yellow solid line represents the percentage of cases in the Pacific region, the brown dotted line represents the percentage of cases in the Orinoquía region, the pink dashed line represents the percentage of cases in the Andean region and the green dotted line represents the percentage of cases in the Amazon region.

https://doi.org/10.1371/journal.pntd.0010228.g008

At the regional level, majority of the fluctuations in R_t were observed during the first 130 epidemic days (until July 5, 2020) after which estimates of R_t remain at around ~1.0. For the Caribbean region the early estimate of reproduction number was R_t~3.42 (95% CrI: 1.85, 5.69) on March 6, 2020. The reproduction number peaked on March 16, 2020, with an estimate of R_t~3.46 (95% CrI: 2.69, 4.36) and declined thereafter to remain consistently above 1.00 (R_t~1.00–3.21) until July 7, 2020. The reproduction number then fluctuated at ~1.00 until February 11, 2021. Since then, the estimates of R_t have remained between ~0.8–1.5 with the most recent estimate of R_t~0.78 (95% CrI: 0.75, 0.81) as of October 31, 2020. In the Pacific region, the early estimate of reproduction number as estimated on March 5, 2020, was R_t~3.04 (95% CrI: 1.28, 5.96). The reproduction number remained above 1.50 in the next 15 days. This was followed by the estimates of R_t consistently above 1.00 until the end of July 20, 2020, after which R_t fluctuated ~1.00 until February 2021. From March to October 31, 2021, the estimates of R_t remained between 0.9–1.5 with the most recent estimate of R_t~1.01 (95% CrI: 0.09, 1.08) as of October 31, 2021. The reproduction number in the Orinoquía region was estimated at R_t~3.46 (95% CrI: 0.25, 14.97) from March 5, 2020, to March 11, 2020, and declined thereafter to fluctuate ~1.00 until the end of May 2020. From June 1, 2020, till the end of August 2020 the estimates of R_t remained above 1.00 followed by fluctuations in R_t around 1.00. From April 9, 2021, the estimates of R_t have remained below 1.0 with the most recent estimate of R_t~0.73 (95% CrI: 0.59, 0.88) as of October 31, 2021.

The reproduction number in the Andean region peaked at R_t~3.28 (95% CrI: 2.27, 4.55) on March 6, 2020. This was followed by the estimates of reproduction number above 1.00 until July 27, 2020. Since then, R_t has fluctuated ~1.00 until mid-June 2021 after which the estimates of R_t remained consistently below 1.0 until October 10, 2021. The most recent estimate of R_t is at ~0.86 (95% CrI: 0.83, 0.88) as of October 31, 2021. In contrast, the reproduction number in the Amazon region was estimated at R_t~3.46 (95% CrI: 0.25, 14.9) on March 5, 2020. This was followed by the estimates of R_t consistently above 1.00 until May 11, 2020 which have since fluctuated ~1.00 until July 2021. The estimates of R_t remained consistently below 1.00 from July 10, 2021, until September 30, 2021, with the most recent estimate of R_t~0.64 (95% CrI: 0.07, 0.77) as of October 31, 2021 (Fig 8). The estimates of R_t obtained from the mortality curve remain higher than 1.00 for the first fifty days after which R_t fluctuates around 1.00 till October 31, 2021 (S4 Fig). As can be observed in both, the national and the regional level case incidence data and the national-level mortality data, the Cori et al. method tends to overestimate R_t early in the time series as infections occurring before the first date in the time series remain missing terms in the denominator [80].

Estimate of reproduction number, R from genomic data analysis.

Between February and April 2020, 136 analyzed Colombian sequences were sampled. These sequences were spread along the whole global SARS-CoV-2 phylogeny (Fig 9) and split into multiple clusters. This indicates multiple introductions of SARS-CoV-2 to the country during the initial pandemic stage (February 27, 2020 to April 5, 2020). For the largest cluster of size 35, growth rate for the early ascending phase of the COVID-19 pandemic in Colombia was estimated to be ~12.718. Subsequently, the reproduction number was estimated at, R~1.2 (95% HPD interval:[1.03,1.44]) indicating sustained disease transmission during the initial phase of the pandemic.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Global Maximum Likelihood (ML) tree for SARS-CoV-2 genomic data analyzed in the paper.

The leaves of the tree are Colombian sequences of SARS-CoV-2 sampled between February 27, 2020, and April 5, 2020 (highlighted in red) and background sequences sampled from around the globe for genomic context. Colombian sequences are distributed among different lineages, indicating multiple virus introductions. The largest estimated intra-Colombian cluster (i.e., a monophyletic clade consisting solely of Colombian sequences) was analyzed using an exponential growth coalescent model to estimate the intra-country basic reproduction number.

https://doi.org/10.1371/journal.pntd.0010228.g009

Spatial analysis.

Fig 10 shows the result from pre-processing COVID-19 data into growth rate functions. The results of clustering are shown in Fig 11. The four predominant clusters identified include the following departments and districts:

1. Cluster 1: Antioquia, Caquetá, Cesar, Choco, Magdalena, Norte Santander, Putumayo, Quindío, Risaralda, Sta Marta D.E., Valle del Cauca
2. Cluster 2: Arauca, Atlántico, Barranquilla, Bogotá, Bolívar, Boyacá, Caldas, Cartagena, Casanare, Cauca, Córdoba, Cundinamarca, La Guajira, Guaviare, Huila, Meta, Nariño, Santander, Sucre, and Tolima
3. Cluster 3: Vichada
4. Cluster 4: Amazonas, Guainía, San Andrés, and Vaupés

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Pre-processing COVID-19 data into incidence rate functions.

From left to right: the original lab-confirmed COVID-19 cases, curve of daily new cases, smoothed and scaled rate curves, and average of rate curves before scaling and smothing.

https://doi.org/10.1371/journal.pntd.0010228.g010

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11.

Panel A. Clustering of states according to the shapes of their rate curves. Cluster 1 is shown in blue, cluster 2 is shown in red, the smallest cluster, cluster 3 is shown in black whereas cluster 4 is shown in green. One can see that states with similar shapes of rates curves are geographically close to each other. The panel shows the geographic distribution of the clusters created in PaintMaps.com, https://paintmaps.com/map-charts/51/Colombia-map-chart [81]. Panel B. Shows the average shape of the incidence rate curves in each cluster and the overall average. Panel C. Shows the dendrogram plot- a hierarchal clustering of states, which shows that there are four predominant clusters of states.

https://doi.org/10.1371/journal.pntd.0010228.g011

Fig 12 shows the mean growth rate curves and one standard-deviation bands around it, in each cluster. The average incidence patterns in these four clusters are very distinct and clearly visible (Fig 12). For cluster 1, the rate increases rapidly from April 2020 to July 2020 and then declines to increase again at the end of November 2020 followed by another incline at the start of January 2021. The most recent upsurge in growth rate can be observed by the two-peak shape from May-July 2021. For cluster 2, we can observe a three modal growth pattern. There is rapid increase in growth rate in July 2020 followed by a decline and continuous fluctuations. This pattern is followed by two more peaks, one in January 2021 and the other peak can be observed in March and April 2021. Cluster 3 shows a slow growth rate until July 2020 followed by a rapid rise in growth curve until mid-September 2020. This is followed by a gradual decline of the growth curve with fluctuations and another large peak between June to July 2021. For cluster 4, the slow growth rate is followed by a spike in September 2020 and multiple fluctuations. This is followed by another peak in June 2021 (Fig 12).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. Average shapes of the COVID-19 incidence rate curves, along with a one standard-deviation band around the average, in each of the clusters.

https://doi.org/10.1371/journal.pntd.0010228.g012

Mobility data

The curve of mobility trends from Apple tracked in the form of driving and walking declined in April 2020 to less than 60% corresponding to multiple interventions focused on the social distancing mandates as implemented by the government of Colombia. Mobility then started to increase in June, peaking in July 2020. Soon after, a decline in mobility was observed in August and September 2020 which was followed by fluctuations in mobility trends around the baseline until the end of June 2021. Thereafter, fluctuations in mobility can be seen trending away from the baseline until October 31, 2021 (Fig 13). Correspondingly it can be observed that as the mobility increased in July 2020, the number of cases also peaked.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 13.

Panel A: Apple mobility trends data showing the variation in walking and driving among the population in Colombia. Orange curve shows the driving trend, and the blue curve shows the walking trend. Panel B: Google mobility trends data showing the variation in movement of people in Colombia across the following categories: grocery and pharmacy (green solid curve), parks (pink solid curve), transit stations (purple solid curve), workplaces (red solid curve) and residential areas (golden solid curve). Panel C: The COVID-19 incidence curve in Colombia by the dates of symptoms onset as of October 31, 2021.

https://doi.org/10.1371/journal.pntd.0010228.g013

The curves of mobility from Google data tracked in the form of visits to grocery and pharmacy, parks, transit stations and workplaces all follow the same pattern; declining in April 2020 to less than 60% incoherence with the Apple mobility data and corresponding to the government’s interventions to contain the spread of the virus. The mobility started to increase gradually thereafter but remained below the baseline until December 2020. The greatest fluctuations in mobility were observed in January 2021, the same time when the cases also peaked. Mobility has since fluctuated around the baseline with an inclining trend in the mobility for transit stations and grocery and pharmacy. The residential category is measured in duration rather than the number of visitors, thus we do not compare it with other five categories. Moreover, since people spend a larger amount of time at home, changes in the residential category remain inconspicuous (Fig 13B).

Twitter data analysis: Social media trends

The epidemic curve for Colombia is overlaid with the curve of tweets indicating the stay-at-home orders in Colombia. In Colombia the engagement of people with the #quedateencasa hashtag (stay-at-home order hashtag) fluctuated at a steady level during the course of the pandemic with a slight decline in the frequency of tweets observed after June 2020 while the number of cases continued to increase, which could point towards the apathy or frustration of the public towards the lockdowns and restrictions. This could also imply that the population does not follow the government’s stay-at-home orders. The decline in the number of tweets indicating stay-at-home orders was followed by a sharp increase in the frequency of tweets in November 2020 which gradually declined until April 2021. A short spike in the frequency of tweets was observed in April 2021 which have since remained around ~500 tweets per day. On the contrary, the cases peaked in July 2020 levelling off at ~7000 cases per day until December 2020 (Fig 14). Another peak in case incidence was observed in January 2021 followed by a sharp decline. This was followed by two upsurges in case incidence in November 2020 and April 2021. The correlation coefficient between the epidemic curve of cases based on the dates of symptom onset and the curve of tweets representing the stay-at-home orders was estimated at R = 0.1797 from March 12, 2020, to October 31, 2021.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 14.

The daily number of tweets indicating stay at home orders (blue) and the daily number of new COVID-19 cases based on the dates of symptom onset (orange) as of October 31, 2021.

https://doi.org/10.1371/journal.pntd.0010228.g014