key: cord-0793368-7pq2gbxm
authors: Albani, Vinicius V. L.; Velho, Roberto M.; Zubelli, Jorge P.
title: Estimating, monitoring, and forecasting COVID-19 epidemics: a spatiotemporal approach applied to NYC data
date: 2021-04-27
journal: Sci Rep
DOI: 10.1038/s41598-021-88281-w
sha: 018d876b39966db2aa03e205f8acc6666c727ce8
doc_id: 793368
cord_uid: 7pq2gbxm

We propose a susceptible-exposed-infective-recovered-type (SEIR-type) meta-population model to simulate and monitor the (COVID-19) epidemic evolution. The basic model consists of seven categories, namely, susceptible (S), exposed (E), three infective classes, recovered (R), and deceased (D). We define these categories for n age and sex groups in m different spatial locations. Therefore, the resulting model contains all epidemiological classes for each age group, sex, and location. The mixing between them is accomplished by means of time-dependent infection rate matrices. The model is calibrated with the curve of daily new infections in New York City and its boroughs, including census data, and the proportions of infections, hospitalizations, and deaths for each age range. We finally obtain a model that matches the reported curves and predicts accurate infection information for different locations and age classes.

Main findings. After smoothing out the daily curves through 7-day moving averages, we estimate the model parameters. The predicted curve by the model for daily new infections shows good agreement with the averaged data curve. Furthermore, the predictions of hospitalizations and deaths match well the reported values based on NYC data and its boroughs.

We observe a dramatic change in the pattern of disease transmission on 19-Mar-2020, identifying the effectiveness of containment measures imposed a week earlier, when a state of emergency was declared and people were asked to stay home. We also observe a considerable decrease in the time-dependent transmission coefficient and the time-dependent basic reproduction rate (obtained via the next-generation matrix technique 34 ). We also noticed this phenomenon in the dynamics of the time-dependent transmission coefficients associated with the NYC boroughs.

In the analysis of the datasets, the patterns of daily hospitalizations and deaths changed consistently between the end of February and the end of August, especially the rate of hospitalization, which has decreased systematically since the end of March. To account for this feature, we allow the model rates of hospitalization and death to be time-dependent. This approach produces model predictions in agreement with the datasets.

Short-term forecasts, with calibrated parameters, were also tested in two different situations, namely, during the transmission regime change and after the spread containment. In both cases, the model accurately predicted the observed scenarios.

Moreover, different reopening scenarios were simulated, considering the impacts of reopening the entire NYC, the borough of Staten Island only, or schools only. In all such cases, unless strict social distancing measures were maintained, the model predictions indicate new infection waves that affect the population of the entire city (Figs. 5, 6, 7) . Such findings are corroborated by recent news, with new infection waves identified in Europe, New Zealand, and China [35] [36] [37] , as well as the reports of COVID-19 spread among a youth population in an overnight camp in Georgia (United States) and in schools in Israel 38, 39 . 

This section presents the epidemiological model and the procedure to estimate the model parameters from the reported data. The estimation of the parameters is performed by minimizing a log-posterior density with a gradient-descent technique. Bootstrap sampling is used to test parameter sensitivity as well as to provide 90% confidence intervals 40 .

The epidemiological model. The SEIR-type model considered here accounts for disease severity, age range, sex, and geographic distribution of some predefined group or population. For simplicity, we postpone the inclusion of sex and (spatial) location dependence to the end of the present section. A number n of age ranges is assumed, each one represented by the superscript i = 1, . . . , n and distributed in seven epidemic categories: susceptible ( S i ), exposed but not yet infective ( E i ), infective in mild conditions ( I i M ), infective in severe condition or hospitalized ( I i H ), infective in critical condition or in an intensive care unit (ICU), denoted by I i I , recovered ( R i ), and deceased ( D i ). Following 10 , we assume the following forms are synonyms: in severe condition and hospitalized. The same applies for the forms of in critical conditions and in ICU. Each individual in the first two infective compartments, mildly infective (M) and hospitalized (H), can recover, die, or develop a more severe disease outcome. The individuals in ICU can only recover or die.

To describe the model, we introduce the following notation: Define the vector where the superscript T denotes the transposed vector, and similarly for E , I M , I H , I I , R , and D . Define also the tensor product:

Then, the epidemiological model can be written as:

A schematic representation of the model in Eqs. (1)- (7) is shown in Fig. 1 . The matrices β M , β H , and β I contain the transmission parameters for the infective individuals in the mild, hospitalized, and ICU classes, respectively. Such parameters are time-dependent and, if well calibrated, may be used to address the effectiveness of the containment measures, to verify the changes in the transmission pattern, or to track the impact of the suspension of a lockdown. It is important to mention that, depending on the information available, simplifying assumptions on the structure of such matrices must be made. In our study, β M , β H , and β I assume the following form:

where the symmetric n × n matrix B is of the form: a 1 b n−1 · · · a n−2 b 2 a n−1 b 1 a n www.nature.com/scientificreports/ β(t) is a time-dependent scalar parameter that controls the epidemic dynamics, and the parameters a and b are scale factors between 0 and 1. The matrix B , which describes the mixing between age ranges, depends on 2n − 1 parameters, namely, a 1 , . . . , a n , containing the observed rates of infections in the n age ranges reflecting the observed heterogeneity of the infectiousness of COVID-19 into the model, and b 1 , . . . , b n−1 , which will be estimated from the available data. The n-dimensional vectors ν M , ν H , and ν I contain the recovery rates for each age range in the infective classes M, H, and I, respectively. Similarly, µ M , µ H , and µ I contain the mortality rates for mild, hospitalized, and in ICU infective individuals. The mean time of incubation is 1/σ . The rates γ M and γ H represent the hospitalization and ICU admission, respectively. The rates of recovery, mortality, hospitalization, and ICU admission are inversely proportional to the corresponding mean times of disease evolution and are directly proportional to the probabilities of moving on to other compartments. All quantities defining such rates are based on the results of references 10, [12] [13] [14] [27] [28] [29] [30] [31] .

The time-dependent transmission parameters β M , β H , and β I , as well as the initial number of infective cases, are unknowns and will be estimated from the recorded data of daily new infections. The available census data are used to determine the population size and the proportions of susceptible population on each age range.

Moreover, whenever the data from daily reports of new cases (infections, hospitalizations, and deaths) include different age ranges, the model can be generalized to incorporate such information. In this case, the entries of β M are as follows:

with β j (t) , j = 1, . . . , n , time-dependent scalar coefficients, and B ij the entries of the matrix B defined above. The other transmission parameters are still of the form β H = aβ M and β I = bβ M .

Since for the NYC datasets only the accumulated numbers of infections, hospitalizations and deaths are agestructured, we assume that the daily reported cases are not age-structured. To introduce more realistic death and hospitalization rates, we adjust µ I and γ M by appropriate delayed ratios from daily reports. More precisely, if γ M and µ I represent the mean rates of hospitalization and death, respectively, for each age range i = 1, . . . , n , the constant rates γ i M and µ i I are replaced by where Ĩ M , Ĩ H , and D are the time series from the daily reports of new infections, hospitalizations, and deaths, respectively. In addition, τ M is the mean time of onset to hospitalization and τ D is the mean time from hospitalization to death. We set τ M = 1 , approximating the median value found in 28 , and the parameter τ D is set to τ D = 1 , obtained empirically in the numerical tests. Note that we do not consider the curves of daily reports of ICU admissions because these data are not available in the NYC dataset.

Including sex in the model. COVID-19 affects male and female individuals differently. Depending on the age range, the case fatality ratio is much larger for male individuals 15, 16 . To take sex variance into account in the model, the transmission parameters ( β M , β H , and β I ) are generalized. The transmission matrix for the mild class assumes the following form: Note that the transmission between sexes is accounted by the mean value of the transmission inside sexes. Intuitively, this means that a female individual who maintains social distancing with other females will also continue such containment measures with male individuals. This assumption also applies to male individuals.

Monitoring and forecasting the disease spread and the effectiveness of containment measures in large regions, such as metropolitan areas, states, and countries, are difficult tasks. Heterogeneous distribution of population and differences in the implementation of social restrictions may lead to quite different disease dynamics among different locations. Moreover, people moving between regions can cause new infection waves. Thus, to account for these aspects, an epidemiological model must include the geographical distribution of the population. A number of approaches have been proposed, and a review on this subject can be found in 20 . In particular, SEIR-type models have been repeatedly used to describe the dynamics of human infectious diseases including geographical information. For example, 23,25 describe the COVID-19 dynamics in United States counties and in Italy, respectively.

To include the geographical distribution of the population into the model described by Eqs. (1)-(7), we enlarge the transmission matrices. By indexing each site under consideration by l = 1, . . . , m , let β l M , β l H , and β l I be the corresponding transmission matrices. Then, the transmission matrix for mildly infective individuals in the model becomes:

where 1 is the m × n-matrix with all entries equal to one and c l = min i,j [B l ] i,j , where B l is the matrix defining the transmission matrix β l M , for the l-th location. The transmission matrices for the other infective classes are similar. This choice for the matrix that represents the mixture of infective populations from different locations helps to simplify the model, considerably reducing the number of unknowns, facilitating calibration. In addition, the model structure is data-driven in the sense that it depends on the current information and thus reflects more precisely the behavior of the population under different containment measures.

When dealing with large areas, such as states and countries, it is important to consider the distance between locations in the model by using exponential, Gaussian, or power law functions 20 . Due to the interconnectedness of NYC and its boroughs, we prefer to estimate the transmission-matrix components from the reported data.

Estimation procedure. For simplicity, we start by presenting the estimation procedure for the model without sex or geographical dependence. Moreover, the data on new infections, new hospitalizations, and new deaths released by the NYC authorities do not include sex or age ranges. Thus, we use this simpler version of the model, where β M = β(t)B . In addition, to simplify the estimation, the constants a and b, related to the transmission matrices of hospitalized and in ICU individuals, are empirically set as a = 0.1 and b = 0.01 , respectively. To estimate the model parameters from the publicly available curves of daily new infections, we build the so-called posterior distribution relating parameters to data.

We assume that the number of daily new infective cases, denoted by I , is Poisson-distributed with parameter σ n i=1 E i (t) . Thus, denoting the vector of model parameters by θ , the logarithm of the likelihood function is where N is the number of the samples in the data and the log(I!) is approximated by the Stirling's formula

We also assume that the vector of parameters θ is Gaussian-distributed with the mean given by some vector of suitably chosen a priori parameters, denoted by θ 0 , and identity covariance matrix. Thus, the negative of the logarithm of the posterior distribution lP(θ|I, θ 0 ) satisfies The constant α is the so-called regularization parameter in the context of Tikhonov-type regularization methods 1 . The estimated parameters are obtained by minimizing lP(θ|I, θ 0 ).

We estimate the initial proportion of mild infective individual in each age range as follows: where I M,0 is a scalar and p i is the population fraction of infective individuals in the i-th age range. The latter is estimated from census data. Thus, the vector of the parameters to be estimated assumes the following form:

where the values of b j , j = 1, . . . , n − 1 , correspond to the entries of the matrix B in Eq. (9) . The time-dependent transmission coefficient β(t) also appears in the definition of β M . The estimation of θ is carried out as follows:

1. Assume that β(t) is constant and estimate θ from the set of daily reports of new infections; 2. Estimate β(t) for each t j in the dataset by minimizing the following functional:

The estimation of the model with geographical information is carried out as follows. Consider the m-dimensional time series of daily new reported infections from m different locations, where I l , l = 1, . . . , m , denotes the set of reports for the l-th location. For each l, let θ l and β l (t) denote the vector of parameters and the time-dependent transmission coefficient, respectively, of the l-th location. We estimate the set of parameters � = [θ 1 , . . . , θ m ] and the coefficients β(t) = [β 1 (t), . . . , β m (t)] T , by minimizing the log-posterior densities given by:

where lP(θ l |I l , θ l 0 ) and F(β l (t j+1 )|β l (t j ), θ l , I l (t j+1 )) are given by (12) and (14), respectively. The computational cost of estimating this model varies with the number m of the sites considered. Based on the available computational resources, for large m, it may be useful to simplify the model by reducing the number of age ranges or merging the locations into a larger areas, thus reducing the model's dimension.

We implemented the model's solution and the estimation procedures in MATLAB (The MathWorks, Inc., Natick, USA). The code is available upon request. The optimization of the posterior density was performed by the general-purpose gradient-based algorithm LSQNONLIN from MATLAB's Optimization Toolbox.

We start by evaluating the accuracy of the proposed methodology on fitting real data. Then, we perform backtesting of the estimated results with out-of-sample data for periods of 7 and 20 days. Finally, we perform a forecast analysis under different scenarios.

Our results are based on New York City reports of new infections, hospitalizations, and deaths. This dataset is updated daily and contains information about the disease distribution at the five NYC boroughs with agestructured accumulated numbers 27 .

As the number of daily COVID-19 tests have been increasing in NYC and more effective treatments have been tested and implemented 41 , the rates of new hospitalizations and new deaths have been decreasing relative to the rates of new infections. To account for all these features, the rates used in the model were updated.

Estimation results. Our initial example set does not yet consider geographical dependence. This information shall be incorporated in "Including geographical information". The disease dynamics are estimated from the number of daily new infections in the entire NYC area. The data series is from 29-Feb-2020 to 21-Aug-2020 and is obtained from the publicly available data in 27 . This source provides COVID-19 case reports and statistics for NYC and each of its five boroughs. The populations of NYC and of its boroughs are distributed in the 5 age ranges present in the datasets. The population distribution in the age ranges and boroughs is based on the census data publicly available at 33 . The curves of the daily reports of new infections were smoothed-out by a moving average of seven consecutive days. The rates per 100,000 inhabitants were used to define the various model parameters. They include the hospitalization, recovery, and death rates, as well as the vector a in the definition of the transmission matrices.

After preliminary calibration tests, two disease evolution regimes were clearly identified. In the first regime, the number of infective individuals increased exponentially, and in the second regime, the spread was considerably reduced due to the containment measures imposed by the state of emergency declared on 12-Mar-2020. The effect of such intervention was clearly observed by the change in the time-dependent transmission coefficient β(t) on 19-Mar-2020 (Fig. 2) .

To capture possible regime changes, such as the different age-range mixing, we divided the time-series into two parts, namely, the data before and after 19-Mar-2020. For these two time series subsets, we estimate the vector of parameters θ and the time-dependent β . Note that for the second dataset, after 19-Mar-2020, we do not estimate the initial infective population. The estimated parameters and the corresponding 90% confidence 40 . Figure 2 presents the comparison between the reported and model predicted curves for daily new infections. The time evolution of the basic reproduction rate can be found in Fig. 2 . Table 2 depicts the predicted and reported number of infections, hospitalizations, and deaths for 21-Aug-2020 with 90% CIs. In Table 1 , the estimated values of β were far from each other in both periods, indicating the aforementioned change in the transmission regime. For the vector b , after containment, the estimated values decrease consistently, indicating that the interaction between age ranges was also significantly reduced. All these results show that disease spread was contained after 19-Mar-2020. However, as we shall see later on, new infection waves can occur if containment measures are relaxed. Figure 2 shows the adherence of the calibrated model predictions to the 7-day moving average of the reported number of new infections. The hospitalization and death rates were evaluated using the ratios of the reported data defined in Eq. (10) . The proportions of ICU admissions by age were obtained in 14 and were adjusted according to the proportions of deaths available in 27 . The model accuracy is also illustrated in Table 2 , showing that the model predictions for the accumulated numbers of infections, hospitalizations, and deaths, for each age range and sex, are close to the reported values. We also evaluate the relative error of the best fit model infections. It has a median value of 6.5 × 10 −4 and a 90% CI of 8.9 × 10 −6 to 0.09.

It is clearly observed from Fig. 2 that the time-dependent basic reproduction rate R(t) values can be classified into two different categories. Prior to 19-Mar-2020, R(t) has large values, indicating that the disease was spreading uncontrollably. After 19-Mar-2020, the transmission parameter value drops to approximately one, indicating the control of transmission by containment measures imposed from 12-Mar-2020 onwards. Note that the large values for the basic reproduction rate in the first part of the series, i.e., prior to 19-Mar-2020, may be caused by an accumulation of reports in the beginning of the outbreak. Such accumulation can be related, for example, to the difficulties faced by the health authorities in setting up an appropriate diagnosis protocol.

The adherence of the calibrated model predictions to reported data, the accuracy in the number of hospitalizations and deaths, as well as the behavior of the calibrated parameters indicating the effect of disease containment measures for the NYC data show that our proposed model captures the COVID-19 dynamics in NYC well. Therefore, it is useful to track the spread dynamics, allowing to assess the effects of, for example, travel quarantine, social distancing, and reopening strategies. If infection curves for different age ranges are available, it is possible to use the present model to track aspects such as the effects of reopening schools, universities, or parks and public gardens since these spaces are usually frequented by people of well-defined age ranges. Thus, Fig. 3 . Figure 3 shows the model predictions adherence to the curves of reported daily new infections. The high accuracy of the model is also evidenced by the comparison between the reports and predictions of the accumulated number of total infections, hospitalizations, and deaths for 01-July-2020 in Table 3 .

The behavior of the time-dependent transmission coefficient for each borough is similar to the basic reproduction rate R(t) in the previous example. This behavior is expected since transmission containment measures were undertaken since 12-Mar-2020 in the entire NYC.

This example demonstrates the ability of the present model to detect disease transmission patterns in different locations at the same moment. The model also accounts for the interaction between individuals of different age ranges, sexes, and locations. Using such features, it is possible to track the implications of reopening, the necessity of additional containment measures, or modification of the design of vaccination strategies. Such broad applications could not be achieved via simpler models. The median and 90% CI of the relative error of the best fit model infections are 6.7 × 10 −4 and 2.0 × 10 −5 to 0.03, respectively.

To test the short-term forecast capabilities of the model, we consider two different periods of the COVID-19 outbreak in NYC.

Uncontained spread. We calibrate the parameters with data from reports on new infections in the period from 29-Feb-2020 to 19-Mar-2020, and we produce a seven-day forecast starting on 20-Mar-2020. This forecasted period is of particular interest since on 19-Mar-2020, the disease spread pattern changed considerably due to the containment measures undertaken 7 days earlier, as a consequence of the state of emergency declared on 12-Mar-2020. To generate the predictions, we assume that β(t) is constant for dates t after 19-Mar-2020 and that it takes the same value as that estimated on 19-Mar-2020. Figure 4 presents a comparison between the forecasted curves and the reported data. Although the parameters were estimated with information based on uncontrolled disease spread, the 7-day ahead forecast for daily new infections and new hospitalizations starting on 20-Mar-2020 show satisfactory accuracy. The accumulated number of infections, hospitalizations, and deaths for the forecasted period can be found in Table 4 .

Contained spread. We now generate a forecast for the period 24-July-2020 to 12-Aug-2020. For dates t after 23-July-2020, the time-dependent values of the transmission coefficient are given by the mean of values for the period from 13 to 23-July-2020. Table 4 presents the predicted and the reported values of accumulated infections, hospitalizations, and deaths from 24-July-2020 to 12-Aug-2020.

As Fig. 4 and Table 4 show, the model predictions of infections, hospitalizations, and deaths once again show satisfactory accuracy. These results are explained by the model ability to incorporate the disease dynamics through the time-dependent parameters. Note that in these examples, the rates of hospitalizations and deaths are evaluated using appropriate ratios of reported data, as defined previously, until the last day of estimation. In the forecasted period, we repeat the rates for the days 19-Mar-2020 and 23-July-2020, respectively, for the corresponding data ranges studied.

Reopening scenarios. We now apply the calibrated models from the previous sections to a number of plausible scenarios, such as the reopening of the entire NYC region, reopening of schools, and reopening of only one single borough, which we chose to be Staten Island. The predictions for these scenarios will illustrate the predictive capability of our model.

The entire NYC. The aim of this example is to present possible scenarios for the COVID-19 epidemic for long periods without an effective vaccine or appropriate treatment. We consider two different scenarios. In the first scenario, the transmission parameters stay at the level of strict containment, as was observed in the period from 04 to 14-June-2020. Thus, for any date t after 14-June-2020, the transmission coefficient β(t) is set as the mean for the estimated values of β(t) in the period from 04 to 14-June-2020, i.e., β(t) = 1.77 with 90% CI 0.27-4.69. In the second case, we simulate a controlled reopening by allowing the coefficient β(t) to reach twice the values obtained in the previous case, i.e., β(t) = 3.54 (0.53-9.38). However, if we impose the condition that whenever the number of daily new infections reaches 1000 cases, containment measures are undertaken, forcing β(t) to return to lower levels until it reaches the value 1.77 (0. 27-4.69 ). On the other hand, after undertaking containment measures, if the number of daily new infections is below 200 cases, we permit reopening, and β(t) may increase again until it reaches 3.54 (0.53-9.38). Figure 5 (left and center panels) present the curves of daily new infections for the aforementioned scenarios. Before 14-June-2020, we have the estimated curves. The forecasting is carried out from 15-June-2020 to 11-May-2021. According to this example, the reopening without an effective vaccine or the achievement of herd immunity may give rise to new infection waves, even when the number of daily new cases is relatively low. To simulate a controlled reopening of schools from 15-June-2020, we use the set of parameters estimated after 19-Mar-2020. To artificially increase the interaction between school-age individuals, the entries of the transmission matrix β M associated with the mildly infective individuals in the age range 0 to 17 years old are multiplied by 2.5. In addition, the time-dependent transmission coefficient is set to 1.25 times the mean of the estimated values for the period 06 to 15-June-2020. For this specific age range, the transmission parameter values are similar to those obtained for the period before 19-Mar-2020. This result is expected since controlling the mixing in youth population at schools is difficult. Indeed, recent news indicates that the COVID-19 transmission rate among people under 19 years old is similar to that in other age ranges 38, 39 . Reopening Staten island. Two different scenarios are considered in this example: reopening Staten Island with and without restrictive measures. In the first scenario, containment measures are slightly relaxed, without allowing people from different age ranges and boroughs to interact. Quantitatively, we keep the same parameter values estimated in the period 19-Mar-2020 onwards. The only change is in the transmission coefficient for Staten Island, which is set to twice the mean of the corresponding estimated values of β(t) for the period 22-June-2020 to 01-July-2020, for dates t after 01-July-2020. During the forecasted period (02-Jul-2020 to 29-Oct-2020), the transmission coefficients for the other boroughs are kept equal to the mean of the estimated values for the period 22-June-2020 to 01-July-2020.

In the second scenario, people from different boroughs and age ranges are allowed to interact, keeping some social distancing and simple containment measures. However, in Staten Island, only simple containment measures are undertaken. In other words, the time-independent transmission parameters assume the same values estimated in the period 29-Feb-2020 to 19-Mar-2020. In addition, we allow the time-dependent transmission coefficient for Staten Island to reach twice the mean of the estimated values for the period 22-June-2020 to 01-July-2020, for dates t after 01-July-2020. Again, after 01-July-2020, the transmission coefficients for the other boroughs are kept equal to the mean of the estimated values for the period 22-June-2020 to 01-July-2020. Whenever the number of daily new infections reaches 1000 in NYC, containment measures are undertaken again. Therefore, the time-independent transmission parameters are brought back to the same values of the period after 19-Mar-2020, and the transmission coefficient for Staten Island is reset to the mean of the corresponding estimated values for the period 22-June-2020 to 01-July-2020. Reopening reoccurs whenever the number of daily new infections in NYC is below 100 reports. Figures 6 and 7 present the curves of daily new infections for the first and second scenarios, respectively. They show the evolution for each borough and for the entire NYC. In the first case, doubling the mean of the transmission coefficient alone is not sufficient to cause secondary waves of infection. In other words, if restriction of movement between boroughs is maintained and the interaction between age ranges remains contained through strict social distancing measures, infection will not return and the disease outbreak will die out. On the other hand, reopening a borough and allowing people of different ages and from different boroughs to interact, even while keeping some light containment and social distancing measures, can cause new waves of infection for large periods.

We use a 7-day moving average to perform smoothing of the data. After calibration, the model predictions were adherent to the data of daily new infections and accurately predicted the number of daily new hospitalizations and deaths. The adjustment of the rates of hospitalization and death using appropriate ratios of reported data contributed to improving model predictions. While backtesting, forecasts for few days ahead under different contexts were found to be accurate. The model accurately identified the effects of the lockdown undertaken in NYC after 12-Mar-2020. A considerable change in the values of the transmission rates was observed. This change flattened the curve and kept the number of daily new infections low, as shown in 27 . www.nature.com/scientificreports/ The rates of hospitalizations and deaths were lower by the end of August than in earlier periods of the COVID-19 outbreak in NYC. This phenomenon may be caused by the changes in disease virulence or in the protocols used to address COVID-19, such as smaller onset to hospitalization mean time, more precise COVID-19 diagnosis, or the introduction of more effective treatments 41 . Moreover, in NYC, the number of COVID-19 daily tests has been increasing consistently since the beginning of the outbreak in the end of February 2020. To account for such changes in the aforementioned rates, the corresponding model parameters were adjusted using the ratios of reported data, incorporating this feature. This considerably increased the accuracy of the model predictions of hospitalizations and deaths, as shown by the results (Fig. 2 and Table 1) .

With regard to reopening strategies, some simulated scenarios generated with calibrated parameters indicate that there is no completely safe method for reopening schools, boroughs, or the entire city unless people respect strict social distancing protocols, avoiding direct personal contact. As model predictions show, even when only a borough or only schools are reopened, new infection waves may occur, forcing public authorities to establish containment measures again. As long as infective individuals are present in a population, there is the risk of new infection waves in reopening strategies since it is not possible to guarantee that everyone will respect the containment protocols. This phenomenon was recently reported 37, 39 . Thus, reopening must be undertaken with strict control of disease transmission, while applying massive testing and enforcing social distancing measures.

Although the role of children and teenagers in COVID-19 spread is still unknown, some recent events of disease spread among a youth population in an overnight camp in Georgia 38 and in schools in Israel 39 indicate that even though most individuals in this age group present mild symptoms, they can infect other people. Therefore, reopening schools also represents a risk of new infection waves.

Even New Zealand and China, which successfully contained COVID-19 spread and had long periods without registering community transmission, are now facing new cases 35, 36 . All these reports corroborate our model predictions, suggesting that without strict control of COVID-19 through social distancing or after massive and effective vaccination, a completely safe reopening may be impossible.

It is important to mention that the present model is also suited to simulate and analyze vaccination strategies because it addresses the dependence of the disease outbreak on age and spatial distribution. This possibility will be addressed in a forthcoming article. program. NYC entered the fourth phase on 20-July-2020. By 22-Aug-2020, schools and shopping malls were still closed, yet public transportation and a number of economic activities were already operational 42, 43 . Strict social distancing measures were still enforced, and public authorities were closely following their observance. This approach kept the number of daily new infections stable. Figure 2 shows the COVID-19 situation in NYC until 21-Aug-2020. The panel presents the comparison of model predictions and reports of daily new infections, as well as the time-dependent basic reproduction rate R(t). We observe that since 19-Mar-2020, R(t) stays near one, meaning that disease transmission is under control, but not eradicated, and that new infection waves may still occur. Even after reaching the fourth phase of controlled reopening, NYC authorities managed to keep transmission under control through social distancing measures by limiting the operational capacity of numerous services, disinfecting public transportation, and many other practices.

Note that even if COVID-19 is eradicated in NYC, as long as no effective vaccine is ready for worldwide use, containment measures inside NYC must be continued to avoid new outbreaks caused by the reintroduction of the disease from abroad.

SEIR-type models have been proposed by a number of authors to predict qualitative aspects of the dynamics of infectious diseases in general and of the COVID-19 pandemic in particular. See 44 for a recent account of the SIR models and their connections to other models. However, to address the elusive aspects of the complex human interactions within the terrain, we think that it is necessary to forego parsimonious models. Functional and high-dimensional models have been used in a number of areas ranging from financial mathematics to population dynamics 45, 46 . They are directly connected to the mathematical theory of inverse problems 1, 2, 22, 23 .

The novelty of the work includes incorporating several aspects of the COVID-19 outbreak that have received little attention in the recent literature and provide sufficient flexibility for excellent adherence to real data. We consider time-dependent rates of transmission, hospitalization, and deaths, as well as the disease age-and sexdependent severity and transmission, while taking into account the spatial distribution of the population.

The model was extensively tested with real data from NYC and its boroughs. After calibration, the model predictions matched the curves of daily new infections, and the model provided accurate predictions for the number of daily hospitalizations and deaths. The relative error of the best-fit model infections was of the order of less than 1% in all cases, indicating the excellent fitting of the models. The model also correctly detected the change in the transmission pattern on 19-Mar-2020 caused by containment measures undertaken on 12-Mar-2020. Moreover, it illustrated the stabilization of the time-dependent basic reproduction rate around the value 1 in NYC.

Concerning the prediction of new infections, the model was also evaluated while using real data and calibrated parameters. It generated accurate results under controlled and uncontrolled transmission contexts. Moreover, different scenarios such as reopening of schools and of an entire borough of NYC were examined. In both cases, transmission rates increased considerably, demanding new containment measures. In other words, without reaching herd immunity or complete disease eradication, we will always face the risk of new infection waves.

Our proposed model is sufficiently general to track transmission dynamics with dependence on age range, sex, and spatial distribution, while evaluating the disease impact on the population. Thus, it is a powerful tool to evaluate scenarios and build proper vaccination strategies.

The datasets analysed are available in the NYC Health website COVID-19: Data https:// www1. nyc. gov/ site/ doh/ covid/ covid-19-data. page.

The code that was used to obtain the findings of the present study is available from the corresponding author upon reasonable request. 

V.A., R.V. and J.Z. designed research; V.A. performed the research; V.A., J.Z. analyzed data; and V.A., R.V. and J.Z. wrote the paper.

of Mathematics and its Applications

An epidemic prediction from analysis of a combined HIV-COVID-19 co-infection model via ABC-fractional operator

Analysis of caputo fractional-order model for covid-19 with lockdown

Mathematical model to assess the imposition of lockdown during covid-19 pandemic

Stability analysis and optimal control of covid-19 with convex incidence rate in Khyber Pakhtunkhawa (Pakistan)

A mathematical model of coronavirus disease (covid-19) containing asymptomatic and symptomatic classes

The effect of messaging and gender on intentions to wear a face covering to slow down covid-19 transmission

Using social and behavioural science to support COVID-19 pandemic response

Estimates of the severity of coronavirus disease 2019: A model-based analysis

Estimating clinical severity of COVID-19 from the transmission dynamics in Wuhan, China

Characteristics of and important lessons from the coronavirus disease 2019 (COVID-19) outbreak in China: summary of a report of 72 314 cases from the Chinese Center for Disease Control and Prevention

Report of the WHO-China joint mission on coronavirus disease

COVID-19)-United States

Sex differential in COVID-19 mortality varies markedly by age

Gender differences in patients with COVID-19: Focus on severity and mortality

Tracking the economic impact of COVID-19 and mitigation policies in Europe and the United States

New predictions suggest a deeper recession and a slower recovery

A contribution to the mathematical theory of epidemics

Modeling Infectious Diseases in Humans and Animals

Data-based analysis, modelling and forecasting of the COVID-19 outbreak

Bayesian dynamical estimation of the parameters of an SE(A)IR COVID-19 spread model

Metapopulation network models for understanding, predicting, and managing the coronavirus disease COVID-19

Inferring change points in the spread of COVID-19 reveals the effectiveness of interventions

Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study

The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: Estimation and application

Baseline characteristics and outcomes of 1591 patients infected with SARS-CoV-2 admitted to ICUs of the Lombardy region

Clinical characteristics of coronavirus disease 2019 in China

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

The impact of COVID-19 vaccination delay: A modelling study for Chicago and NYC Data

On the definition and the computation of the basic reproduction ratio R 0 in models for infectious-diseases in heterogeneous populations

BBC. New Zealand coronavirus: 14 new Covid-19 cases reported

China records biggest one-day rise in coronavirus cases since march

Europe's second wave? country-by-country breakdown of resurging covid-19 cases

SARS-CoV-2 transmission and infection among attendees of an overnight camp

When covid subsided, israel reopened its schools. it didn't go well

Fitting dynamic models to epidemic outbreaks with quantified uncertainty: A primer for parameter uncertainty, identifiability, and forecasts

What restrictions on reopening remain in New York? The New York Times

The challenges of modeling and forecasting the spread of covid-19

A splitting strategy for the calibration of jump-diffusion models

Transport equations in biology

The authors declare no competing interests.

Correspondence and requests for materials should be addressed to J.P.Z.Reprints and permissions information is available at www.nature.com/reprints.Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.