Transmission Dynamics and Short-Term Forecasts of COVID-19: Nepal 2020/2021

Nepal was hard hit by a second wave of COVID-19 from April–May 2021. We investigated the transmission dynamics of COVID-19 at the national and provincial levels by using data on laboratory-confirmed RT-PCR positive cases from the official national situation reports. We performed 8 week-to-week sequential forecasts of 10-days and 20-days at national level using three dynamic phenomenological growth models from 5 March 2021–22 May 2021. We also estimated effective and instantaneous reproduction numbers at national and provincial levels using established methods and evaluated the mobility trends using Google’s mobility data. Our forecast estimates indicated a declining trend of COVID-19 cases in Nepal as of June 2021. Sub-epidemic and Richards models provided reasonable short-term projections of COVID-19 cases based on standard performance metrics. There was a linear pattern in the trajectory of COVID-19 incidence during the first wave (deceleration of growth parameter (p) = 0.41–0.43, reproduction number (Rt) at 1.1 (95% CI: 1.1, 1.2)), and a sub-exponential growth pattern in the second wave (p = 0.61 (95% CI: 0.58, 0.64)) and Rt at 1.3 (95% CI: 1.3, 1.3)). Across provinces, Rt ranged from 1.2 to 1.5 during the early growth phase of the second wave. The instantaneous Rt fluctuated around 1.0 since January 2021 indicating well sustained transmission. The peak in mobility across different areas coincided with an increasing incidence trend of COVID-19. In conclusion, we found that the sub-epidemic and Richards models yielded reasonable short-terms projections of the COVID-19 trajectory in Nepal, which are useful for healthcare utilization planning.


Introduction
South Asia is the most populated region of the world, with a population of about 2 billion. This region includes eight of the world's 36 megacities with an urban population greater than 10 million and population density greater than 10,000 per square kilometer [1]. When the first confirmed case of COVID-19 in the region was reported in Nepal, the public health community feared for the worst given the high poverty and population density levels and a poorly resourced health care system [2]. Nevertheless, South Asia was not significantly affected by the first COVID-19 wave in 2020 compared to many developed countries. Yet, the situation took a different turn when neighboring India reported the emergence of the Delta variant (also called B.1.617) on 5 October 2020 [3], the most transmissible variant of COVID-19 reported as of September 2021 [3][4][5]. Since early March 2021, the Delta variant started to take hold of the Indian population and fueled a devastating second wave that engulfed South Asia. As of 30 November 2021, South Asia reported a total of 47.9 million confirmed cases (18.91% of global total) and 767,195 COVID-19 deaths (15.04% of global total) [6].
Nepal is one of the South Asian countries most heavily affected by the COVID-19 pandemic. By 30 November 2021, the Himalayan country of roughly 30 million people

Data
For short-term forecasting of the COVID-19 epidemic in Nepal, and for estimating reproduction number at the national and provincial level from the case incidence data, we utilized publicly available time series of laboratory-confirmed RT-PCR positive cases by dates of reporting that were obtained from the COVID-19 situation reports published by Ministry of Health and Population of the government of Nepal as of 30 November 2021 [23]. Similarly, to assess the mobility trend during the course of the pandemic in Nepal, we used Google's mobility data for Nepal [24]. Google's mobility data shows how visits to places, such as the grocery stores, parks, and recreation spots, are changing in each geographic region. A baseline day represents the normal value for that day of the week. The baseline day is the median value from the 5-week period from 3 January-6 February 2020. In this study, we analyzed Google's mobility data from 15 February 2020 to 30 November 2021.

Modeling Framework for Forecast Generation
We utilized three dynamic phenomenological growth models to generate short-term forecasts for Nepal (10-days and 20-days ahead). These models have been applied to various infectious diseases including SARS, foot and mouth disease, Ebola [25], and the current COVID-19 outbreak [26,27]. Phenomenological growth models applied in this study include the following differential equation models: the generalized logistic growth model [25], the Richards growth model [28], and the sub-epidemic model [29]. The forecasts obtained from these dynamic growth models can provide an assessment of the potential scope of the outbreak in near real-time and insight on the impact of the implementation and relaxation of control interventions, which help guide public health policy. The description of these models is provided as follows.

Generalized Logistic Growth Model
The generalized logistic growth model (GLM) [25] displays a range of epidemic growth profiles, including the polynomial and exponential growth patterns. GLM characterizes epidemic growth by estimating three parameters: (i) the intrinsic growth rate, r (ii) a dimensionless "deceleration of growth" parameter, p and (iii) k 0 , the epidemic size. The varied epidemic growth patterns are observed by the modulation of the deceleration of growth parameter resulting in the exponential growth dynamics (p = 1), the sub-exponential growth (0 < p < 1), or the constant incidence (p = 0) patterns. The GLM model is given by the following differential equation: where C(t) denotes the cumulative number of cases over time t, and dC(t) dt describes the incidence curve over time t.

Richards Growth Model
The Richards model [28] extends the simple logistic growth model by incorporating a scaling parameter, a, that measures the deviation from the symmetric simple logistic growth curve [28,30,31]. The Richards model is given by the differential equation: where C(t) represents the cumulative case count at time t, r is the growth rate, a is a scaling parameter, and k 0 is the final epidemic size.

Sub-Epidemic Model
The sub-epidemic model [29] can be used to model an epidemic wave comprising of multiple overlapping sub-epidemics where each sub-epidemic is modeled using a GLM model denoted by the following differential equation: We model an epidemic wave comprising of n overlapping sub-epidemics using a system of coupled differential equations, as follows, where C i (t) is the cumulative cases for the ith sub-epidemic, and k i is the size of ith subepidemic where i = 1 . . . n. Hence when n = 1, the sub-epidemic model reduces to the simple logistic-type model. To model the onset timing of the (i + 1)th sub-epidemic, we use an indicator variable A i (t), making sure that the sub-epidemics comprising an epidemic wave follow a regular structure. Therefore, where 1 ≤ C thr < k o and A 1 (t) = 1 for the sub-epidemic 1. For the subsequent subepidemics, the size of the ith sub-epidemic (k i ) declines exponentially at a rate q. This can happen due to multiple factors such as the seasonal transmission effect, effect of interventions and population behavior changes. If q = 0, then the sub-epidemic model predicts an epidemic wave composed of equal sized sub-epidemics [29]. If we assume that the subsequent sub-epidemic sizes decline exponentially, we get k i = k 0 e −q(i−1) , where k 0 is the size of initial sub-epidemic i.e., k 1 = k 0 . Therefore, when q > 0, the total number of sub-epidemics supported by the model depends on C thr , q and k 0 because the (i + 1)th sub-epidemic is triggered only if C thr > C i (t) [29].

Model Calibration and Forecasting Approach
We conducted 8 week-to-week subsequential 10-days and 20-days ahead short-term forecasts at the national level utilizing the data retrieved from Ministry of Health and Population, Government of Nepal by the dates of report [23]. Table 1 shows the calibration and forecast period for each of the three dynamic models used. For each of the models, we estimated the best fit solution using non-linear least square fitting procedure [32]. This process minimizes the sum of squared errors between the model fit f (t,Θ) and the smoothed data estimates y t , and yields the best set of parameter estimates Θ = (θ 1 , θ 2 , . . . , θ m ), where the smoothed data was obtained by using smoothfactor = 7 in Matlab. The smooth function in Matlab reduces the noise within a data set by using a moving average method. LetΘ = argmin ∑ n t=1 f (t,Θ − y t ) 2 denote the best fit estimates. HereΘ = (r, p, k o ) corresponds to the set of estimates for the parameters of the GLM modelΘ = (r, k o , a) corresponds to estimates of parameters of the Richards model, and Θ = (r, p, k o , q, C thr ) corresponds to the estimates of parameters of the sub-epidemic model [30]. While for the sub-epidemic model, we provided the initial best guesses of the parameter estimates, for GLM and Richards growth model, we initialized the parameters for the nonlinear least squares' method [32] over a wide range of plausible parameters sampled from a uniform distribution. This allowed us to test the uniqueness of the best fit solution. Moreover, the initial conditions were set at the first data point for both models [30]. Uncertainty bounds around the best-fit solution were generated using parametric bootstrap approach assuming a negative binomial error structure for each of the models, where the variance was calculated by averaging mean to variance ratio in the data. The details of this method are provided in ref [30].
For each model, we generated M = 300 datasets by the bootstrap approach during the calibration phase, refitted the model to each generated dataset, and used the M sets of parameters estimates to construct the 95% confidence intervals for each parameter. Further, each M best fitted model was used to generate m = 30 additional data points with negative binomial error structure extended through the forecasting period. For the forecasting period, we constructed the 95% prediction intervals with these 9000 (M × m) curves. Detailed description of the methods of parameter estimation can be found in references [30,33,34].

Performance Metrics
We utilized the following five performance metrics to assess the quality of our model fit and the 10-days and 20-days ahead short-term forecasts: root mean squared error (RMSE) [35], the mean absolute error (MAE) [36], the mean interval score (MIS) [35], the coverage of the 95% prediction intervals [35], and the weighted interval score (WIS) [37] for each of the three models. For calibration performance, we compared the model fit to the smoothed incidence data fitted to the model, whereas for the performance of forecasts, we compared our forecasts with the incidence data for the time-period of the forecast.
The root mean squared error (RMSE) and the mean absolute error (MAE) assess the average deviations of the model fit to the observed data in L 2 and L 2 norm, respectively. The root mean squared error (RMSE) is given by and the mean absolute error (MAE) is given by where y t i is the time series of cases by date of onset, t i is the time stamp, andΘ is the set of estimated parameters. For the calibration period, n equals the number of data points used for calibration, and for the forecasting period, n = 10 and 20 for the 10-day and 20-day ahead short-term forecast respectively. Moreover, to assess the model uncertainty and performance of prediction interval, we used the 95% PI and MIS. The prediction coverage is defined as the proportion of observations that fall within 95% prediction interval and is calculated as where y t i are the case incidence data, L t i and U t i are the lower and upper bounds of the 95% prediction intervals, respectively, n is the length of the period, and I is an indicator variable that equals 1 if value of Y t i is in the specified interval and 0 otherwise. The mean interval score addresses the width of the prediction interval as well as the coverage. The mean interval score (MIS) is given by In this equation, L t i , U t i , n, and I are as specified above for PI coverage. Therefore, if the PI coverage is 1, the MIS is the average width of the interval across each time point. For two models that have an equivalent PI coverage, a lower value of MIS indicates narrower intervals.

Reproduction Number
Reproduction number (R t ) characterizes the average number of secondary cases generated by a primary case at calendar time t during an outbreak when control measures are in place [39]. R t is an important public health indicator of effectiveness of public health interventions during an epidemic [39]. While R t estimates of more than 1 indicate continuation of widespread disease transmission, R t less than 1 indicates that sustained disease transmission is unlikely and the outbreak is under control [40].

Estimating Reproduction Number (R t ) Using GGM
We estimated the effective reproduction number (R t ) for the early ascending phase of the COVID-19 epidemic in Nepal using a generalized growth model (GGM). We estimated reproduction number (R t ) at the national level for first 30 days of the phase during initial rise in cases (25 May to 23 June 2020), and first wave (1 August to 30 August 2020) of COVID-19 pandemic. For the second wave, we estimated reproduction number for the first 30 days, at both the national and at provincial level. Based on the growth of incidence curve, we estimated R t for the period from 12 April to 11 May 2021 for national level and for six provinces, except for Karnali province for which we took data for the period (20 April to 19 May 2021). Similarly, we also estimated R t for the resurgence wave after the second wave for the first 30 days period (10 July to 8 August 2021) at the national level.
We modeled the generation interval of SARS-CoV-2 assuming gamma distribution with a mean of 5.2 days and a standard deviation of 1.72 days [41]. We estimated the growth rate parameter r, and the deceleration of growth parameter p from the GGM. The GGM model is used to simulate the progression of local incidence cases I i at calendar time t i . This is followed by the application of the discretized probability distribution of the generation interval, denoted by ρ i , to the renewal equation to estimate the reproduction number at time t i [42][43][44] 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 generation interval. The factor α assesses the relative contribution of imported cases to secondary disease transmission. The factor α was set at 0. The numerator represents the total new cases I i , and the denominator represents the total number of cases that contribute (as primary cases) to generating the new cases I i (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 R t are derived directly from the uncertainty associated with the parameter estimates (r, p) obtained from the GGM. We estimated R t for 300 simulated curves assuming a negative binomial error structure where the variance was estimated by averaging mean to variance ratio for the local case incidence [30].

Estimating Instantaneous Reproduction Number (R t )
The instantaneous reproductive number is the expected number of secondary infections occurring at time t, divided by the number of infected individuals, each scaled by their relative infectiousness at time t. We estimated the instantaneous reproduction number, R t , at the national and regional level using the case incidence date by dates of report, using the method of Cori et al. [45] in which the instantaneous R t is estimated as where I t is the number incident infections on day t and w s is the generation interval at time since infection s. An individual's relative infectiousness depends on generation interval and time since infection (s) [45]. Here the term ∑ t s=1 I t−s w s describes the sum of infection incidence up to time step (t − 1) weighted by the current infectiousness (w s ) of individuals who became infected s days in the past, and who may be shedding the virus now. The standard assumption is that w s follows a discretized gamma distribution [45]. However, since the infection time may not be accurately observed, the measurement of generation time becomes difficult [46] and so the time of symptom onset is usually used to estimate the distribution of serial interval (SI), which is the time interval between dates of symptom onset among two successive cases in a disease transmission chain [45]. We assumed the generation interval as equal to SI and that SI follows a gamma distribution with a mean of 5.2 days and a standard deviation of 1.72 days [41]. We then obtained the R t estimates over weekly time interval within a Bayesian framework using EpiEstim R package in R language [45]. We reported the mean and 95% credible interval (CrL).

Results
As of 30 November 2021, Nepal reported a total of 821,366 RT-PCR positive COVID-19 cases and a total of 11,526 deaths. Figure 1 shows the incidence curve of COVID-19 cases and deaths in Nepal as of 30 November 2021. The COVID-19 epidemic curve shows three distinct phases: an initial growth phase occurring during late-May to late-June 2020, the first wave occurring from early-August 2020 to mid-January 2021, and the second wave taking off in mid-April 2021 ( Figure 1). The first wave started to peak after the lifting of the first national lockdown, which lasted about four months. The second wave was preceded by the second wave of COVID-19 in India, a celebration of Hindu massive festivals, and the start of wedding season in Nepal. The second wave came to a decline in late June 2021, but it resurged a few weeks later ( Figure 1). Epidemiologia 2021, 2, FOR PEER REVIEW 8      In Table 2, we present the comparison of performance metrics of the three models for 10-days ahead forecasts. For 10-days ahead forecast, GLM model performed better in RMSE and MIS, while sub-epidemic model performed better in PI-Coverage. Overall, we saw that in the first three forecast periods, the sub-epidemic model out-performed the GLM and Richards models in all the metrics. Likewise, the Richards model out-performed GLM and the sub-epidemic model in all the metrics in the fifth and sixth forecast period, while the GLM model performed better in the eighth forecast period ( Table 2). Dates in each column indicate date of the end of 30-days calibration period. Table 3 presents the comparison of performance metrics of the three models for 20-days ahead forecast. For 20-days ahead forecast, GLM model performed better in RMSE, and sub-epidemic model performed better in MIS, and PI-coverage in the majority of the forecast periods. Similar to the 10-days ahead forecast, the 20-days ahead forecast also observed the Richards model performing better in most of the metrics.

Estimate of Reproduction Number, from Case Incidence Data Using GGM
The reproduction number for the early ascending growth phase of the initial mild wave of epidemic from the case incidence data (25 May to 23 June 2020) using GGM was estimated at ~1.1 (95% CI: 1.1, 1.2) for the national data. The growth rate was 13 (95%

Estimate of Reproduction Number, R t from Case Incidence Data Using GGM
The reproduction number for the early ascending growth phase of the initial mild wave of epidemic from the case incidence data (25 May to 23 June 2020) using GGM was estimated at R t~1 .1 (95% CI: 1.1, 1.2) for the national data. The growth rate was 13 (95% CI: 7, 20) and the deceleration of growth rate parameter (p) was estimated at 0.41 (95% CI: 0.35, 0.48). Likewise, for the first 30 days of the first wave from 1 August to 30 August 2020, the reproduction number was estimated at 1.1 (95% CI: 1.1, 1.2), the growth rate was estimated at 15 (95% CI: 8.5, 20), and the deceleration of growth rate parameter was estimated at 0.43 (95% CI: 0.39, 0.49). The deceleration of growth rate parameter for the initial mild wave and the first wave (p~0.41-0.43) indicate almost a linear pattern of epidemic trajectory of COVID-19 in Nepal.
For the second major wave, we estimated the reproduction number (R t ) at national and provincial level. For the early phase of the second wave (first 30 days), R t at national level was 1.3 (95% CI: 1.3, 1.4) with the growth rate of 7.9 (95% CI: 5.8, 11) and p of 0.61 (95% CI: 0.58, 0.64) (Figure 9, Table 4), indicating a sub-exponential growth dynamic of COVID-19 pandemic. At provincial level, while province 1, Gandaki, and Sudurpaschim had the highest R t at 1.5, Lumbini and Karnali regions observed the lowest R t at 1.2 ( Figure 10, Table 4). For province 1, Bagmati, Gandaki, and Sudurpaschim, the deceleration of growth rate parameter indicated sub-exponential growth dynamics with estimated between p~0.61-0.72. For other provinces, p~0.51-0.58 indicated an almost linear trajectory of the COVID-19 epidemic (Table 4).  (Table 4).     We also estimated the reproduction number for the 30 days period from 10 July-8 August 2021, for the early ascending phase of the resurgence that started to peak from early July 2021. The reproduction number for the early ascending growth phase of the resurgence phase after the second wave from the case incidence data (10 July to 8 August 2021) using GGM was estimated at Rt~1.2 (95% CI: 1.2, 1.2) for the national level. The growth rate was 18 (95% CI: 12,20) and the deceleration of growth rate parameter (p) was estimated at 0.47 (95% CI: 0.46, 0.52). Figure 11 shows the weekly instantaneous reproduction number at national and provincial level for the period from 1 October 2020 to 30 September 2021. The reproduction number was estimated at, ~2.5 (95% CrI: 2.4, 2.5) on 16 April 2021, which declined to 0.99 (95% CI: 0.98, 1.0) on 10 May 2021, and increased to 1.04 (95% CI:1.02, 1.06) on 9 July 2021 at national level (Figure 11). At the provincial level, highest was observed in Karnali region with an estimate of ~5.08 (95% CrI: 2.04, 9.48) on 26 February 2021. The instantaneous reproduction number fluctuated at around ~1 from December 2020 to January 2021 for all the regions. Similarly, instantaneous reproduction has remained consistently below 1 for national level and fluctuated between 0-1.5 at provincial level since 8 August We also estimated the reproduction number for the 30 days period from 10 July-8 August 2021, for the early ascending phase of the resurgence that started to peak from early July 2021. The reproduction number for the early ascending growth phase of the resurgence phase after the second wave from the case incidence data (10 July to 8 August 2021) using GGM was estimated at R t~1 .2 (95% CI: 1.2, 1.2) for the national level. The growth rate was 18 (95% CI: 12,20) and the deceleration of growth rate parameter (p) was estimated at 0.47 (95% CI: 0.46, 0.52). Figure 11 shows the weekly instantaneous reproduction number at national and provincial level for the period from 1 October 2020 to 30 September 2021. The reproduction number was estimated at, R t~2 .5 (95% CrI: 2.4, 2.5) on 16 April 2021, which declined to 0.99 (95% CI: 0.98, 1.0) on 10 May 2021, and increased to 1.04 (95% CI:1.02, 1.06) on 9 July 2021 at national level (Figure 11). At the provincial level, highest R t was observed in Karnali region with an estimate of R t~5 .08 (95% CrI: 2.04, 9.48) on 26 February 2021. The instantaneous reproduction number fluctuated at around~1 from December 2020 to January 2021 for all the regions. Similarly, instantaneous reproduction has remained consistently below 1 for national level and fluctuated between 0-1.5 at provincial level since 8 August 2021.

Estimate of Instantaneous Reproduction Number, R t
Epidemiologia 2021, 2, FOR PEER REVIEW 17 Figure 11. Instantaneous reproduction number for the period from 1 October 2020-30 November 2021.

Analysis of Mobility Data
The curves of mobility from Google data tracked in the form of visits to retail and recreation, grocery and pharmacy, parks, and workplaces all follow the same pattern, showing inclining trends from late May to late June 2020, which correspond to the peak of initial rise in cases of COVID-19. The mobility curve started to incline again in August, which corresponded to the early phase of the first wave. The second wave started soon after there was a peak in mobility around parks and grocery and pharmacy. During the peak in the second wave, mobility declined rapidly in all the areas except residential, which corresponds to the lockdown by government as an intervention to contain the epidemic ( Figure 12).

Analysis of Mobility Data
The curves of mobility from Google data tracked in the form of visits to retail and recreation, grocery and pharmacy, parks, and workplaces all follow the same pattern, showing inclining trends from late May to late June 2020, which correspond to the peak of initial rise in cases of COVID-19. The mobility curve started to incline again in August, which corresponded to the early phase of the first wave. The second wave started soon after there was a peak in mobility around parks and grocery and pharmacy. During the peak in the second wave, mobility declined rapidly in all the areas except residential, which corresponds to the lockdown by government as an intervention to contain the epidemic (Figure 12).

Discussion
Reliable forecasting tools are of critical importance in guiding public health policy during the period of an ongoing pandemic [47]. Appropriate short-term forecasts not only help assess the impact of interventions but also help guide the distribution of limited resources. In this study we compared the performance of three different phenomenological models for short term forecast of COVID-19 cases based on the estimates derived using case incidence data at national level. We performed eight different week-to-week sequential forecasts of 10-days and 20-days. Our findings indicate the better performance of subepidemic model based on five performance metrics during the model calibration phase. However, for the 10-days and 20-days ahead forecast performance, we found the subepidemic model performing better in the initial forecast periods and the Richards model performing better in the later forecast periods. Our finding is in line with previous studies indicating that the sub-epidemic model frequently outperforms the GLM and Richards model [48]. Our forecast estimates from the three phenomenological models indicate a declining trend of COVID-19 cases in Nepal as of June 2021.
The early estimate of in different epidemic waves reflect a sustained transmission of COVID-19 in Nepal with above 1. We report almost a linear pattern of COVID-19 incidence during the initial transmission phase and the first waves in Nepal (deceleration of growth parameter p~0. 41-0.43) with a reproduction number at 1.1 (95% CI: 1.1, 1.2). We also report a sub-exponential growth pattern in the second wave with a deceleration of

Discussion
Reliable forecasting tools are of critical importance in guiding public health policy during the period of an ongoing pandemic [47]. Appropriate short-term forecasts not only help assess the impact of interventions but also help guide the distribution of limited resources. In this study we compared the performance of three different phenomenological models for short term forecast of COVID-19 cases based on the estimates derived using case incidence data at national level. We performed eight different week-to-week sequential forecasts of 10-days and 20-days. Our findings indicate the better performance of subepidemic model based on five performance metrics during the model calibration phase. However, for the 10-days and 20-days ahead forecast performance, we found the subepidemic model performing better in the initial forecast periods and the Richards model performing better in the later forecast periods. Our finding is in line with previous studies indicating that the sub-epidemic model frequently outperforms the GLM and Richards model [48]. Our forecast estimates from the three phenomenological models indicate a declining trend of COVID-19 cases in Nepal as of June 2021.
The early estimate of R t in different epidemic waves reflect a sustained transmission of COVID-19 in Nepal with R t above 1. We report almost a linear pattern of COVID-19 incidence during the initial transmission phase and the first waves in Nepal (deceleration of growth parameter p~0.41-0.43) with a reproduction number at 1.1 (95% CI: 1.1, 1.2). We also report a sub-exponential growth pattern in the second wave with a deceleration of growth parameter p at 0.61 (95% CI: 0.58, 0.64) and a corresponding reproduction number at 1.3 (95% CI: 1.3, 1.3). We found much variation in the reproduction number, growth rate, and deceleration of growth rate parameters at the provincial level, with province 1, Gandaki and Sudurpaschim provinces having R t at 1.5 during the initial phase of the major second wave. Our estimate of reproduction number is less than that of India [49,50] and comparable with early estimates of reproduction number reported from Mexico [48] and South Korea [51], and higher than that of the early estimates reported from Singapore [52].
As of 24 August 2021, the SARS Cov-2 variants of concern (VOC) reported from Nepal include Alpha and Delta variants [53]. The Delta variant is at least twice as contagious as previous variants [54]. The national public health laboratory of Nepal confirmed the detection of the Delta variant in 47 of the 48 swab samples collected from infected individuals from different parts of the county from 9 May to 3 June 2021, and in nine of these samples, the highly infectious new sub-lineage K417N (also called AY.1) was detected [55], indicating the presence of Delta variant since at least May 2021. However, unlike expected higher value of reproduction number for the period with the circulating Delta variant, our estimate of instantaneous reproduction number did not show an increase in R t above 1, both at national level and provincial level since 25 May-12 June 2021 ( Figure 11). This could be due to the vaccination efforts by the country during early/mid 2021 [56,57] and also due to inadequate COVID-19 testing rates around that time, due to a large number of cases that may have been undetected [9].
Our analysis of mobility trend shows the peak in mobility across different areas relative to the trajectory of the COVID-19 incidence (Figure 12). The increase in mobility before the third wave corresponds to the start of wedding season [14], celebration of Hindu massive festivals in Nepal, crowded movement of people in the country, a lack of adherence to social distancing guidelines [13], and the start of second wave in India (Figure 1). Our findings from google mobility data are in line with the country's severity index of public health and social measures (PHSM) published by the WHO regional office for South East Asia ( Figure 13) [58]. Figure 13 shows that decline in the severity index for PHSM was followed by the surge in daily COVID-19 cases. The lifting of bans on businesses, gatherings, stay at home, and public transport during February and March 2021 preceded the start of second wave during April 2021 in Nepal (Figure 1, Figure 13).
Epidemiologia 2021, 2, FOR PEER REVIEW 19 growth parameter p at 0.61 (95% CI: 0.58, 0.64) and a corresponding reproduction number at 1.3 (95% CI: 1.3, 1.3). We found much variation in the reproduction number, growth rate, and deceleration of growth rate parameters at the provincial level, with province 1, Gandaki and Sudurpaschim provinces having at 1.5 during the initial phase of the major second wave. Our estimate of reproduction number is less than that of India [49,50] and comparable with early estimates of reproduction number reported from Mexico [48] and South Korea [51], and higher than that of the early estimates reported from Singapore [52].
As of 24 August 2021, the SARS Cov-2 variants of concern (VOC) reported from Nepal include Alpha and Delta variants [53]. The Delta variant is at least twice as contagious as previous variants [54]. The national public health laboratory of Nepal confirmed the detection of the Delta variant in 47 of the 48 swab samples collected from infected individuals from different parts of the county from 9 May to 3 June 2021, and in nine of these samples, the highly infectious new sub-lineage K417N (also called AY.1) was detected [55], indicating the presence of Delta variant since at least May 2021. However, unlike expected higher value of reproduction number for the period with the circulating Delta variant, our estimate of instantaneous reproduction number did not show an increase in above 1, both at national level and provincial level since 25 May-12 June 2021 ( Figure 11). This could be due to the vaccination efforts by the country during early/mid 2021 [56,57] and also due to inadequate COVID-19 testing rates around that time, due to a large number of cases that may have been undetected [9].
Our analysis of mobility trend shows the peak in mobility across different areas relative to the trajectory of the COVID-19 incidence (Figure 12). The increase in mobility before the third wave corresponds to the start of wedding season [14], celebration of Hindu massive festivals in Nepal, crowded movement of people in the country, a lack of adherence to social distancing guidelines [13], and the start of second wave in India (Figure 1). Our findings from google mobility data are in line with the country's severity index of public health and social measures (PHSM) published by the WHO regional office for South East Asia ( Figure 13) [58]. Figure 13 shows that decline in the severity index for PHSM was followed by the surge in daily COVID-19 cases. The lifting of bans on businesses, gatherings, stay at home, and public transport during February and March 2021 preceded the start of second wave during April 2021 in Nepal (Figure 1, Figure 13).  [58]. Indicators and PHSM severity index are assigned a color. The shading of that color is based on the severity score. The darkest shade represents the most severe measure/s (value of 100) and a lack of color represents no measure (value of 0). Nepal depends on India for most of its supplies, including medical equipment, liquid oxygen, and vaccines. Therefore, when the situation deteriorated in India's health system with a total of 26,031,991 COVID-19 cases and 291,331 deaths, as of 21 May 2021, Nepal struggled to find alternate sources of supplies beyond India [15]. India's surge in COVID-19 cases impacted the availability of COVID-19 vaccines in Nepal. Nepal launched its COVID-19 vaccination program on 27 January 2021 after receiving a first batch of one million doses of AstraZeneca (Covishield) vaccines produced by the Serum Institute of India [15]. Nepal also received doses from Covax and from China [15]. However, Nepal's attempt to buy additional doses of AstraZeneca vaccine from India failed due to the increase in demand in that country [16]. As of 30 November 2021, 27.8% of the total population in Nepal have been fully vaccinated and 32.5% have received at least one dose [7]. Given the vaccination rate, fragile health system, direct impact of the COVID-19 trajectory in India (in terms of cases flow due to open border as well as in terms of management as seen during the vaccination), and a sustained COVID-19 transmission, there is a need to implement non-pharmaceutical measures in order to reduce the reproduction number.
Our study has several public health implications. Short-term forecasts using the sub-epidemic and Richards models yield the most reliable short-term forecasts. This can be useful in planning the interventions during the early growth phase of the epidemic at present and for similar epidemics in the future. Furthermore, the analysis of the subnational estimates of the reproduction numbers could help prioritize subpopulations on interventions, especially for a resource-scarce country such as Nepal.
Our study has some limitations. The case positivity rate fluctuated around 10% during the second week of April 2021 and rose to 73% between mid-April and mid-May [16], suggesting that only severe symptomatic cases were being tested and a much larger population could have been infected. Therefore, our estimates of reproduction number may be underestimated. Similarly, we used only RT-PCR positive case counts in our study while the national situation reports provide the total number of cases includes both RT-PCR and antigen positive cases. Of note, the descriptive analysis of the data in most of situational reports included only the RT-PCR positive cases. Considering only RT-PCR cases will make our findings comparable with other countries [59]. Moreover, the cases incidence and mortality data used in this study are based on the reporting date, rather than date of symptom onset or date of death. Therefore, the difference in the report date and actual date of event occurrence, including the differences in testing rates and reporting delay, might affect our model estimates and estimates of reproduction number.

Conclusions
The reproduction number in Nepal has been fluctuating around 1 since January 2021 indicating a sustained virus transmission in the country. We report a sub-exponential growth dynamic during the early phase of the second wave with a reproduction number at 1.3 at the national level with province 1, Gandaki and Sudurpaschim provinces having higher reproduction numbers (R t at 1.5) during the initial phase of the major second wave. The increase in mobility in the country appears to be positively correlated with COVID-19 transmission. The sub-epidemic and Richards model provide reasonable short-term projections of the COVID-19 trajectory in Nepal and indicate a declining trend of COVID-19 cases in the country until June 2021. Simple mathematical models can provide reliable short-term projections, which are crucial for public health planning.

Conflicts of Interest:
The authors declare no conflict of interest.