key: cord-0738889-2o7m89r5
authors: Kevrekidis, P. G.; Cuevas-Maraver, J.; Drossinos, Y.; Rapti, Z.; Kevrekidis, G. A.
title: Reaction-diffusion spatial modeling of COVID-19: Greece and Andalusia as case examples
date: 2020-05-09
journal: Physical review. E
DOI: 10.1103/physreve.104.024412
sha: a593e7ff4acc5d927f4ffa3d17666799fc2e4883
doc_id: 738889
cord_uid: 2o7m89r5

We examine the spatial modeling of the outbreak of COVID-19 in two regions: the autonomous community of Andalusia in Spain and the mainland of Greece. We start with a 0D compartmental epidemiological model consisting of Susceptible, Exposed, Asymptomatic, (symptomatically) Infected, Hospitalized, Recovered, and deceased populations. We emphasize the importance of the viral latent period and the key role of an asymptomatic population. We optimize model parameters for both regions by comparing predictions to the cumulative number of infected and total number of deaths via minimizing the $ell^2$ norm of the difference between predictions and observed data. We consider the sensitivity of model predictions on reasonable variations of model parameters and initial conditions, addressing issues of parameter identifiability. We model both pre-quarantine and post-quarantine evolution of the epidemic by a time-dependent change of the viral transmission rates that arises in response to containment measures. Subsequently, a spatially distributed version of the 0D model in the form of reaction-diffusion equations is developed. We consider that, after an initial localized seeding of the infection, its spread is governed by the diffusion (and 0D model"reactions") of the asymptomatic and symptomatically infected populations, which decrease with the imposed restrictive measures. We inserted the maps of the two regions, and we imported population-density data into COMSOL, which was subsequently used to solve numerically the model PDEs. Upon discussing how to adapt the 0D model to this spatial setting, we show that these models bear significant potential towards capturing both the well-mixed, 0D description and the spatial expansion of the pandemic in the two regions. Veins of potential refinement of the model assumptions towards future work are also explored.

incorporating the ingredients of a generalized SIR formulation. At the same time, a complementary viewpoint that is far less computationally expensive but possibly quite informative in its own right is a metapopulation network approach in the spirit of the work of [41] . While networks and metapopulation studies, see for example [21, 42, 43] , are useful in attempting to examine small numbers of groups of individuals, it is clear that these cannot properly capture the scale of an entire country. The number of nodes would simply be too huge and it would not be possible to capture the tremendous variations of population-density scales. Such models could be well suited to study propagation of waves of an epidemic between metropolitan centers, but they surely are not well suited to capture different and much smaller scales (at least not without significant adaptations). For instance, they cannot capture dynamics that happens within a node (unless they have further structure) and they surely cannot capture dynamics that (spatially) happens in-between nodes. The interest in such models is to consider transport along links if one is interested in long-range transport of a virus by, e.g., airline travel or along highways. Here, we seed viral hotspots and explore how the virus will spread locally thereafter. Our infection initialization by hotspots aims to emulate the initial long-range transport of infectious individuals, and hence it becomes an indirect way to incorporate mobility on a network in the absence of convection. In that sense, our technique too requires significant adaptations but for a different reason: this is in order to properly capture mobility of the population that induces the spreading of a pandemic.

Moreover, one can envision techniques (such as the equation-free modeling framework [44] ) that may enable the cross-linking of the above two approaches, e.g., the use of metapopulation network systems to perform PDE-level tasks. Lastly, there exist isolated examples of models that take into account the structure of different types of networks. A particularly nice example in this direction is the work of [45] , which leverages the availability of Enron, Facebook and social graphs in the form of adjacency matrix patterns that can be used to represent the connectivity within a country's network.

With these considerations in mind, we develop an expanded variant of the classical SIR model (ODE model) and then focus on its PDE spatio-temporal generalization. We incorporate particularities of this virus, such as its latent period, i.e., that individuals exposed to the virus may be infected but not infectious during the latent period, and the significant fraction of infected and infectious individuals that do not develop symptoms. In Section II we present the spatial model, analyze first its ODE variant that will be used to perform appropriate optimization of its parameters in the cases of the country of Greece (but without the islands, i.e., the mainland of Greece), and the Spanish autonomous community of Andalusia. While our study has focused on both regions, for practical purposes, we opted to relegate the presentation and discussion of our results for Greece to an Appendix. This renders the presentation of our model, our approach, and methodology easier to follow, and more focused. The selection of two seemingly unrelated regions may appear a bit disparate, yet we argue these to be particularly interesting examples. Aside from their intrinsic interest to the authors, these roughly equally-sized regions, with similar population densities, exhibit significantly different, i.e., by an order of magnitude, number of deaths, illustrating the potential impact of different policies. Upon optimizing the ODE results, we use their output to formulate the input of the corresponding PDE framework and explain how to set it up within the software package COMSOL. Following the formulation of both the ODE and PDE approaches in section II, the results for Andalusia and the respective interpretations and comparisons with reported date are offered in section III. We provide numerical results for key categories such as cumulative infections and deaths, comparing the PDE results both with the available data for these regions and the associated ODE results. Finally, in section IV we summarize our findings and present our conclusions offering a number of possibilities towards future work. The first Appendix summarizes the calculation of the basic reproduction number by the next-generation matrix approach, whereas the second presents model results for the mainland of Greece for well-mixed and spatially-distributed populations, mimicking our analysis of the spread of the virus in Andalusia.

We first explain the ODE model which is obtained from the full PDE model by removing the convection and diffusion "spatial aspects" in the convection-reaction-diffusion equations of interest.

At the level of an ODE formulation, a relevant extension of the standard SIR model can incorporate some of the key features of this virus, such as, e.g., that a fraction of the exposed population remains asymptomatic [46] . We thus start with a population of susceptibles (S), which may become exposed (E) upon the emergence of the virus within the population. This represents the well-documented [47] feature that the virus is latent within the host for a period of time, before he/she becomes infectious (able to transmit the virus to susceptible hosts). After this latent period, exposed individuals, in turn, may become asymptomatically infectious (A) at rate σ A , or (symptomatically) infected (I) at rate σ I . We assume here that both A and I can interact with the susceptibles S with respective rates β AS and β IS to draw new members of the population into the group E of individuals exposed to the virus. We note that the transmission rates β incorporate the total population size (ODE model) or the total population density (PDE Model). A fraction of hosts in the I class may need hospitalization, thus giving rise to a population of hospitalized (H) at rate γM . Among these, a fraction responds to the treatments, thus leading ultimately to a population of recovered (R) at a rate (1−ω)χ. At the same time, the seriously ill who are hospitalized yield also a number of deceased (D) at a rate ωψ. Asymptomatically infected hosts recover at a rate M AR (i.e., asymptomatic recovered) and move into class AR and the seriously ill recover at a rate (1−γ)M . While AR could, in principle, be merged with R, in our view, it is meaningful to maintain these two populations separate since R is a measurable quantity within available COVID-19 data (and, hence, comparable to model predictions), while AR is not. Notice that the above constants reflect both the population fraction partition (e.g., ω vs 1 − ω) and the (inverse) time scales (e.g., χ vs. ψ) for transition between subgroups.

A weak effect of net change of the population due to other birth or other mortality factors (−µS) can be incorporated in the susceptibles and can be adequately assessed from census data, yet we do not incorporate it in the D population aiming to evaluate purely the deaths stemming from COVID-19. Here, we briefly note two points. The pool of susceptible individuals is not significantly affected (over the time scale of our study) from this term which can be safely neglected for our purposes. Secondly, one can include such a term in the rest of the populations involved in our study. However, the (underlying) health conditions often involved in such mortality often lead to complications in the concurrent presence of COVID-19 and when this leads to mortality, the latter is attributed to COVID-19. Hence, we do not include such a separate term in the rest of the equations.

The above specify the "ODE parameters" within the group; these reflect processes that happen either in an averaged way at the "well mixed" level (when no spatial dependence is assigned) of the ODEs or processes that happen locally at every point in space for the PDEs. We will return to this when we discuss parameter conversions in the next section. The relevant populations and rates of conversion can be seen in a self-contained form in Fig. 1 . The ODE version of the proposed model, described in the figure, is similar to the SEAIHR model used to model the transmission of the MERS coronavirus in the Republic of Korea [48] , but slightly distinct in its treatment of the asymptomatically recovered, the recovered, and the hospitalized who, in the current model, do not transmit the virus, as they are expected to be in isolation. 

We now turn to the PDE properties of the model involving spatial spreading of the pandemic. Initially, we note that we do not anticipate that infected (which should be self-quarantined), hospitalized (or at stages thereafter) will have a diffusivity, i.e., D I = D H = 0 in the initial installment of the model. As regards the R and AR, in principle they can have a diffusivity (although there is a period of recovery), yet since it is fair to assume that these populations have immunity in the immediate interval after their infection, we can assign D R = D AR = 0. However, an interesting possibility within the model is the inclusion of population time-dependent diffusion, possibly anisotropic, and also directed (along the direction of the velocity v) motion, as considered, for example, in [49] where a laboratory case of epidemic propagation along lines of fast diffusion is presented to model the spreading of a virus along a highway. As regards the remaining populations, it may be tempting to examine nonlinear variants where the diffusivity is larger, e.g., where the population is larger, reflecting the existence of a well-established transportation/mobility network. Nevertheless, in the present work, we will initiate relevant considerations by assuming constant diffusivity of the susceptibles, the exposed, and the asymptomatics. The latter are the key, given their mobility and spatial spreading for the corresponding spreading of the pandemic in the context of Eqs. (1)- (8) .

An additional important decision that can be incorporated at the level of the PDE model concerns the functional form of the directional velocity v. In principle, this can be used to capture "daily practices" (e.g., going to work, spending time there, commuting back and resting practices), but also longer temporal or spatial scales (e.g., trips from city to city, or country to country). Motivated partially by the colloquial understanding of some of the case examples considered such as the spreading of the pandemic in Greece [50] , at the present level, we opt not to incorporate these effects but simply allow diffusion to perform the relevant spreading. The idea within a given region then is that arriving infected individuals, e.g., from international travel, form local hotspots within the S population and we examine the diffusional spreading effect of the virus in the presence of the above local viral dynamics. Hence, the initial infectious seeding within the susceptible population is a rough approximation intended to emulate long-range transport in the absence of convection. We will see that this approach is not unreasonable given the results that we obtain for the spreading of the PDE results with both the ODE ones and the data available online for the cumulative infections and the deaths within the regions of interest. Naturally, it is hoped that this will be a seed study towards a further refinement of such considerations on the basis of more accurate spatial data for the spreading of the disease.

In the results given in the next section, we have selected as our illustrative example the autonomous community of Andalusia within Spain. The example of the mainland of Greece is presented in Appendix B. While these examples may seem somewhat disparate, they bear some significant advantages as regards their nature and their comparison. First off, they are regions of similar populations of about 8-10 million inhabitants. Greece has been praised in international media [51] regarding its handling of the first-wave COVID-19 crisis and the effectiveness and promptness of the associated social-distancing measures. Additional relevant features of this region include (a) day 0 of the infection first wave and (b) the origin of the localized events thereof could be successfully identified, as well as (c) strict lockdown effects went into place early on. Another example at the opposite end with very significant numbers of infections and deaths is Spain. However, here there is a significant set of complications. Not only is Spain far larger in spatial and population size, but importantly for the number of reported cases and especially the number of deaths, there is no universally accepted way of reaching the relevant conclusive numbers across the 17 different autonomous communities. For all of the above reasons, and also for reasons of clearer comparison of comparable sizes (and also for ones of intrinsic interest to the authors, admittedly), we selected the autonomous community of Andalusia.

Having selected our target regions, the next complication is to formulate the solution of Eqs. (1)-(8) at the level of the autonomous community/country as a "two-dimensional spatial grid". That is one significant complication toward spatial modeling which we have addressed by utilizing the finite element package COMSOL Multiphysics [52] . We have inserted the regions' map as a geometry within COMSOL and proceeded subsequently to form a triangulated mesh of the computational domain.

The next and also rather complex step is to formulate a population as an initial condition of susceptibles within the relevant grid. Here, we have leveraged tools from the large scale geographic project World Pop [53] . This methodology encompasses census data and enables via random forest models [54] the generation of a gridded prediction of the population density at a resolution of about 90 m. We have imported this type of data within our spatial country grids and via interpolation we are in a position to simulate models of the type of Eqs. (1)-(8) with arbitrary choices of parameters, and, in principle, also initial conditions. This is, in our view, a significant combined asset (the spatial grid of a region combined with an interpolated over this grid realistic representation of population census data) towards modeling spreads of epidemics.

The crucial next step, within this line of modeling the spreading of the epidemic, is to identify suitable parameters, similarly to what has been done in numerous earlier studies [55, 56] at the ODE level. To do so, we utilized a nonlinear optimization algorithm such as the constrained minimization, fmincon function within Matlab. We determined the optimal model parameters by minimizing the Euclidean distance N ( 2 norm) between the time series generated by the model, identified by the subscript "num", and the corresponding "observed" (data) time series, identified by the subscript "obs",

where the index "i" identifies a point in the time series. The parameters were optimized to reproduce the time series of the reported total number of infected cases (C(t) = I(t) + H(t) + R(t) + D(t), the total number of "cases") and total number of deceased (D(t)). We found these two time series to provide the most reliable data. Specifically, for the case of Greece we note some nontrivial lapses in the apparent curation of the data. Particularly noteworthy is the case of the recovered individuals in [50] . The data must evidently be significantly inaccurate, as the number of recovered individuals appears to stay fixed at 53 between March 29 and April 5, only then to jump entirely abruptly to 269 recovered, only to stay there between April 6 and April 29, then to jump on to 1374. Admittedly, the unprecedented circumstances were straining the data collection process, yet it is particularly important to provide accurate data to modelers to calibrate adequately the models towards the future spreading of the pandemic. It is these two data columns (total cases and deaths) from [50] that we thus compare to our 0D model for Greece and the columns that were used in the parameter optimization. As expected for this large parameter space, initial parameter choices and parameter constraints affect the parameters resulting from the optimization algorithms. In the next section, where we present the ODE parameters, we discuss a number of sensitivity studies as well as related issues of parameter identifiability. In addition, we argue that the suggested median values are biologically and socially reasonable, and that they are in line with a number of features known about the SARS-CoV-2 transmission.

At the level of parameters within a certain individual and how the virus acts on it "on average", i.e., as concerns parameters such as (σ A , σ I , M, γ, M AR , ω, χ, ψ), we preserve the same values at the PDE level as at the ODE one. The transmission rates β are more complicated. Keeping in mind that at the PDE level the quantities, S, E, etc. are no longer populations, but rather population densities [57] which integrate over the region's spatial surface (through the respective surface integrals) to the true population of each category, we can immediately infer that the units of such densities are proportional to l −2 where l is a characteristic length-scale of the analysis. In that vein, the β's need to be multiplied by l 2 to adapt dimensionally between the ODE and the corresponding PDE model. Indeed, we found these to be the most complicated parameters to select at the level of the PDE model, as we will explain in the discussion of the results below. It is important to bear in mind that while the results below are given for these two regions our aim is to develop a set of tools that could be in principle used, alongside with data for the pandemic from different countries [3], to perform similar analyses of other regions.

We start the exposition of our results by discussing what we will refer to as the "0D" model (the version of Eqs. (1)-(8) without space dependence) for Andalusia. Data for the evolution of the pandemic in Andalusia were obtained from the Andalusian-Goverment COVID-19 site [58]. The relevant results are given in Fig. 2 . We obtained the optimal (best-fitting) 0D-model parameters for Andalusia (and Greece) from 2,000 optimizations that compared model predictions to jointly the number of cumulative infected and the deceased, as shown in Eq. (9) . For each optimization, the initial guess for each parameter and initial condition was uniformly sampled within a pre-specified range. The upper and lower limits were used as boundaries in the constrained minimization algorithm (implemented in Matlab via the fmincon function). The parameter ranges were determined from epidemiological information. We note that at the initial time of model fitting, the number of exposed and asymptomatic individuals is not known. We thus optimized (and varied) their ratio to the initially infected I(0), a number that was obtained by subtracting the officially reported number of deaths, recovered, and hospitalized from the (reported) number of cases. The sensitivity of the predicted model parameters when the ratio of β AS to β IS is allowed to vary within a specified interval ([0.2, 2]), in steps of 0.02, and parameter-identifiability issues are discussed at the end of this Section.

Upon performing the optimizations we find that the fitting yields the results summarized in Table I . We show the median parameters, as well as the interquartile range, and the range of variation used to sample the parameters (and initial conditions). Model predictions (with median parameter values, solid blue or red lines) are compared graphically to data (black dots) in Fig. 2 . Model output sensitivity to parameter (and initial-condition) variations is represented by the shaded regions.

The optimal parameters were obtained for two scenarios (see Table II ). The first scenario considers that restrictive measures (quarantine, lockdown) in Andalusia were strictly enforced on March 16, 2020. To account for the change in parameters induced by the lockdown, we imposed a time dependence on the transmission rates β

so that the transmission rates β IS and β AS decrease by a factor η IS and η AS (respectively) relatively abruptly at the time t q the lockdown was imposed. The transmission rates effectively incorporate the rate of contact of susceptible individuals with infectious individuals (infected or asymptomatic in our model) that leads to exposure to the virus. In fact, the transmission rates β may be expressed as the product of the average daily contacts (contact rate) a susceptible has with any individual times the probability of infection given a contact: the probability of infection is proportional to the viral load (viral concentration in the respiratory-tract fluid) of expelled respiratory droplets [8] . Hence, government-imposed restrictive measures, e.g., mobility restrictions, social-distance requirement, face-mask wearing, and limitations on the number of persons in a gathering, are expected to decrease the transmission rates (for both infected and asymptomatics), an effect that is reflected in the η's. The two top panels of Fig. 2 illustrate that how well we capture the data for Andalusia depends on the scenario chosen, i.e., when the quarantine is imposed (the events time series is summarized in Table II ). The red solid line (top line at t=150) and shaded region correspond to imposing the quarantine at t q = 3 (scenario one, modeling the beginning of the quarantine almost immediately when it officially occurred), whereas the blue solid line (bottom line at t=150) and shaded region present data (median and range) for t q = 16, scenario two. Both scenarios reproduce reasonably well the number of the fatalities (top right panel), scenario two better, but not so the number of cases. Moreover, the difference between predictions and data increases with simulation time, with the predictions of scenario one becoming progressively very high and dissonant with the trend of the data. There is a characteristic feature in the top left panel (cases) that the first scenario fails to capture: there is an "angle" in the semi-logarithmic plot associated with the curbing of the observed cumulative number of infections C(t) due to containment measures. It is apparent that the attempt to capture the data, due to the relevant mismatch in the associated angle, leads to far more significant deviations. While model prediction seems to minimize the distance to the data by over-predicting C(t) initially, and under-predicting it later, it clearly starts over-predicting the trend of the quantity towards the end of the available (fitted) data time-series. This results in predicted cumulative infections of the order of several (more than 4) tens of thousands. A similar over-prediction seems to develop in D(t) leading to nearly five thousand deaths, while the data seem to clearly tend to values below that.

If we were to shift arbitrarily the time of the application of the quarantine data by 13 days later (t q = 16), we note a nontrivial difference. While we are not missing on D(t) (in fact, the fit is more accurate), we capture accurately the angle in the C(t) data. This suitable shift of the quarantine time clearly does a far better job in capturing the actual trends of both C(t) and D(t) with the C(t) lying between 10 4 and 2 × 10 4 and, correspondingly, D(t) staying below 2 × 10 3 at the end of the fivemonth simulation period. Table I Bottom panels show the other populations. Bottom left panel shows these populations for tq = 3 (distinguishable at t=150 from top to bottom: recovered R(t), asymptomatic recovered AR(t), hospitalized H(t), asymptomatic A(t), and exposed E(t)). The bottom right panel shows them for tq = 16 (distinguishable at t=150 from top to bottom: asymptomatic recovered AR(t), recovered R(t), hospitalized H(t), asymptomatic A(t), and exposed E(t)). In all the panels, shaded regions correspond to the interquartile range for each quantity, whereas the full line corresponds to simulations with the median parameter (and initial-condition) values.

reported in the first and second columns. The most noteworthy difference is that while for scenario one (no reproduction of the angle) β AS > β IS this inequality reverses in scenario two. A concomitant change occurs in the ratio of exposed who turn asymptomatics (σ A /(σ I + σ A )) from 0.41 (scenario one) to 0.56 (scenario two). The relative importance of these changes is discussed at the end of this section where parameter identifiability is addressed. Lastly, note the slight change in the recovery period χ. A potential interpretation of this admittedly somewhat arbitrary shift of the quarantine-parameter imposition may be that at the model level such measures have an immediate, essentially instantaneous effect, while in the realistic country data, there is a time lag before this switch in the number of contacts (due to lockdown) has a perceptible effect, depending on how fast individuals adapt to the imposed restrictions.

Since the second scenario reproduces more satisfactorily the observed data we discuss (with the proviso mentioned at the end of this section) the biological and societal significance of its median parameters in the second column of Table I . For instance, the median latent period is approximately 3 days (2.976), whereas the median incubation period is approximately 4 days (3.77), in reasonable agreement with the values reported in [47] , 3 and 5 days, respectively. The value of M AR suggests a time scale of nearly 7 days (6.85) for the recovery of asymptomatics. On the other hand, M suggests a time scale of about 6 (5.97) days for those with symptoms to potentially need hospitalization. The median timescales associated with leaving the hospitalized compartment imply almost 7 days (1/χ = 6.74) to recovery and almost 8.5 days (1/ψ = 8.33) to fatality. Hence, the approximate recovery period for mild cases is about 6 days, and for severe cases approximately 13 days, while Ref. [47] estimates recovery periods of approximately 2 weeks for mild cases and approximately 6 weeks for the quite severe cases. The value of γ roughly suggests a half-half split between those recovering directly versus those needing some form of hospitalization. The value of ω suggests that among those needing hospitalization nearly 75% recover, while only 25% die. As discussed above, there is no way to evaluate the initial population of exposed E 0 and asymptomatics A 0 . We thus opted to introduce their fraction to the initial population of infected I 0 as two additional parameters (referred to as initial conditions), E 0 /I 0 and A 0 /I 0 in the optimization. Our optimization yields initially (at the beginning of the simulations) approximately equal ratios of exposed and asymptomatics to the infected.

The relative transmission rates of asymptomatics and infected, as well as the ratio of their populations, merit a comment. The reported median values satisfy β IS > β AS , i.e., infected (with symptoms) are predicted to be more infectious than asymptomatics. We surmise that the contact rate would be significantly smaller for the (expected to be) self-isolating infected individuals, than for the asymptomatics who continue their life, not knowing that they are carrying SARS-CoV-2 (and most importantly that they are infectious). Hence, their higher transmission rate would imply a higher emitted viral load. Related to the transmission rates is the ratio of asymptomatics to symptomatically infected. Assuming that the latent time scale of the virus is similar for asymptomatics and infected (as is reasonable to assume), the fraction of turning asymptomatic versus turning infected (σ A /σ I ) is 1.27. Equivalently, the ratio of becoming asymptomatic to the total number of exposed (σ A /(σ A + σ I )) is about 0.56 and that to becoming symptomatically infected 0.44. This ratio reflects the importance of asymptomatics [46] in the transmission of SARS-CoV-2, a particularly important feature that differentiates it from the transmission of other respiratory viruses like influenza and SARS-CoV-1. Lastly, after lockdown measures are imposed, we find that both populations are equally affected, η IS ≈ η AS , possibly because mobility restrictions, and the associated decrease in the average number of daily contacts, apply equally to both populations.

One final observation at the level of data rather than at that of the model is the significance of the early imposition of restrictive measures. In Spain these measures were taken when already the number of cumulative infections and deaths was significantly higher than the corresponding numbers in Greece when the decision was taken. This ultimately appears to have led the smaller of the two regions (Andalusia having 8.4M inhabitants) to have an order of magnitude larger losses of life and infections than the larger of the two regions (Greece having 10.7M inhabitants).

We also calculated the basic reproduction number R 0 reflecting the number of cases expected to be produced by one infectious case in a fully susceptible population. This is used to estimate how the epidemic developed initially. We used the next-generation matrix approach (see Appendix A). The pre-quarantine basic reproduction number (for the second scenario, t q = 16), is R 0 = 1.91(1.86 − 1.95), where we report both the median and the interquartile range. The calculated basic reproduction number is close to the epidemiologically determined range of 2 − 4 [47] , a range that encompasses the variation of the basic reproduction number in space and time. Although our focus is not on the calculation of R 0 , which is of principal interest to a wide range of studies regarding SARS-CoV-2, it is worth noting that a post-quarantine effective reproduction number may be calculated for the lockdown-decreased transmission rates. We find R eff = 0.763(0.72−0.74) reflecting the decline of the epidemic spreading under the lockdown measures. A similar calculation with the scenario-one parameters yields a pre-quarantine R 0 = 2.65 (2.38−2.90), and, interestingly, a post-quarantine effective reproduction number R eff = 1.04 (1.03 − 1.04) a value that suggests the epidemic has not been effectively controlled (as implied by the predicted evolution of the epidemic in Fig. 2, top left) .

We conclude the analyses of the ODE model by commenting on the identifiability of model parameters. We partially addressed it through the previously presented parameter sensitivity analysis with 2000 optimizations. This procedure led to the determination of the median values and their interquartiles. In addition, and as argued in [59] , one approach to specifying confidence intervals is through the Hessian of the variation of the Euclidean norm Eq. (9), the objective function of our optimizations, with respect to model parameters [60] . Specifically, if we denote model parameters by θ the Hessian is H ij = ∂ 2 N /∂θ i ∂θ j , suggesting that if it remains invariant to parameter changes, these parameters would not be identifiable (since their changes would not modify the optimized norm). Alternatively, as discussed in [60] , the inversion of the Hessian leads to the confidence intervals associated with each parameter. When we carried out this programme for a model similar to the 0D model presented here [59] , we found that the Hessian was singular: in fact, it had two zero eigenvalues. Our above line of argumentation (expanded upon in [59] ) suggests that these two "zero-cost" eigendirections are closely connected to the identifiability of the model, and specifically that a number of parameters associated with these eigendirections are not independently identifiable. One eigendirection is easily specified through inspection of the model. The three parameters ω, ψ and χ may be easily combined to two κ 1 = (1−ω)χ and κ 2 = ωψ.

The second combination of parameters that defies identifiability is less immediately transparent. However, we get a hint of the other zero-eigenvalue eigendirection, and the associated not-independently-identifiable parameters by comparing columns one versus two in Table I (and also from the runs of Fig. 3 ; see especially the top panel thereof). When the time the lockdown was imposed is modified (going from scenario one to scenario two), the transmission rates and the time scales shift from β IS < β AS and 1/σ A < 1/σ I (whose ratio, as discussed earlier, determines the fraction of turning asymptomatic to turning symptomatically infected to be 0.70) to the β IS > β AS and 1/σ A > 1/σ I (with the corresponding fraction becoming 1.27). Alternatively, for β AS > β IS the fraction of exposed turning asymptomatic is smaller than when β AS < β IS , i.e., the larger the asymptomatic transmission rate the smaller their fraction. This inverse relation becomes quantitative in the top panel of Fig. 3 where we note that as β AS σ A decreases β IS σ I increases. We chose these two parameter combinations as they appear naturally in the two summands of the basic reproduction number, Eq. (A1). The optimizations, whose results are reported in the figure, were Despite the wide variation of transmission rates (times the associated time scales), we observe a rather small uncertainty in the forward, model predicted evolution of the pandemic, an indication that even though some model parameters are not independently identifiable, they enable an adequate predictor for the quantities of epidemiological interest. performed as previously discussed (i.e., parameters and initial conditions were uniformly sampled within their range of variation) with an additional constraint on the ratio of the two transmission rates β AS /β IS . We chose their ratio to vary between 0.2 to 2, sampled in ten equidistant values. For every ratio of the transmission rates we performed 200 optimizations. The solid blue line denotes the relationship for the median parameters for each choice of β AS /β IS , the red dots correspond to the optimal parameters for each optimization. Monitoring the results obtained in Fig. 3 , it is natural to conclude β AS and β IS are not independently identifiable. Instead, there is effectively a monoparametric freedom (associated with the singular eigendirection) connecting these two parametric combinations.

For this wide range of parameters we also compare the predicted number of cases (bottom left panel Fig. 3 ) and fatalities (bottom right panel) to the reported numbers. Importantly, we note that even though the transmission rates times the associated inverse time scales may vary significantly due to their non-identifiability, the predicted number of cases and fatalities does not (compare also to Fig. 2 ). This is, indeed, a manifestation of the singular eigendirection: a wide range of model parameters provides an equally good predictor of the total number of cases and fatalities. Hence, the uncertainty in the identification of model parameters, and their non-identifiability, has a relatively small effect on the predictions of the model. We believe this provides convincing evidence of the predictive ability of the model and its accuracy. A far more relevant question is not the specification of model parameters and their confidence intervals, but how does the flexibility to specify them, as allowed by the singular eigendirection, modify model predictions. We find that the optimally determined model parameters provide a reasonable, within a given range, estimate of the modeled quantity, even though due to the non-identifiability of the model the model parameters may vary. It is also an interesting direction to explore what additional pieces of data (such as, e.g., on the asymptomatic infections) may render the model identifiable, enabling a more precise identification of the relevant parameters.

We now turn to the PDE simulations. Relevant results for the autonomous region of Andalusia may be found in Fig. 4 for the same diagnostics as for the 0D model. However, now, we complement them with the space-time evolution simulations of Figs. 5-7 that will be compared also in what follows with the data of the map of Fig. 8 . We first explain how we selected the model parameters and how we initialized the PDE model of Eqs. (1)-(8) and then we discuss the numerical results, emphasizing their advantages and deficiencies. At the regional (spatial) level, we must adapt the 0D model parameters. In the PDE model, we retained the same median parameters as the optimized 0D (ODE) model parameters starting with the σ's and beyond in Table I . This is because they involve processes occurring at the level of a single individual, i.e., "locally", and hence we do not expect them to change at the country level in the transition from the ODE to the PDE model. In addition, we kept the same reduction factor of the transmission rates η IS and η AS to model the effect of restrictive measures on them.

On the contrary, we do not expect this to be the case for the transmission rates β. They depend on the interaction between individuals since they may be expressed as the product of the daily average number of contacts times the infectious disease transmission probability [8] . At the ODE level, the presence of S and A or I immediately leads to the conversion of susceptibles to exposed. At spatial (region or county) level, this effect does not occur homogeneously as it does at the ODE level, but rather in a distributed way. As the population is (spatially) distributed over the country in a highly heterogeneous way, the ODE β's have to be modified to obtain their "spatially averaged" variant.

We obtained these spatially averaged transmission rates by first keeping the ratio of the β's the same as that of the ODE, but scaling each one by a scaling factor ξ. The transition from the ODE to the PDE transmission rates involves the introduction of two length scales. The first reflects the transition from the number of individuals (e.g., S, I, E, etc.) in the ODE description to spatial densities of individuals in the PDE description; the other length scale reflects the transition from a spatially homogeneous to a spatially distributed model. We obtained their product by noting that the β's have to be multiplied [57] by an effective inverse density l 2 /N , N being the country population, i.e., by multiplying the ODE transmission rates by the scaling factor ξ ≡ l 2 /N . The product length scale l defines an effective spatial scale over which the ODE transmission rates need to be rescaled to obtain the corresponding PDE transmission rates. The scaling factor was determined by minimizing the 2 norm specified in Eq. (9) . For these optimizations we kept all model parameters at their median values, while diffusivities were varied as subsequently discussed in Eq. (11) .

In addition to the decision regarding the scaling of the β's, an important decision is that of the selection of the diffusivities. Recall that in the present first work we decided to avoid attempting to model convection effects, but rather mostly focus on the role of diffusion. We assume that most of the populations relevant to the infection which have not developed any symptoms, namely the asymptomatics and the exposed, diffuse with a diffusivity of D c = 100 km 2 /day (i.e., associated with a characteristic spatial scale of about 10 km). Our motivation for this choice is that in this small population (for the regions and data considered) associated with the infection, it is relevant to include a wider spatial spread of their motion to enable (through their contacts) the infection to spatially spread. On the other hand, for the far larger population of susceptibles, we assign a smaller diffusivity (0.1D c ) since we consider that the mobility of susceptibles does not significantly change their distribution [49] . In essence the much larger susceptible population provides a background for the relative motion of asymptomatic and exposed carriers of the virus. The rest of the populations (most notably, I, H, and D, since the immunity of R and AR renders their diffusion inconsequential) are assumed to be highly localized/self-isolating and hence bear, for our purposes, a vanishing diffusivity. For all the non-vanishing diffusivities, we assume that the quarantine reduces them to a fraction η D of their original value (see Table I ) in a similar ramped form as before for the transmission rates:

We should also describe the initialization of the model. We selected to populate initially eight key "hotspots" of the infection as they arose in Andalusia. The selected areas (Almería, Córdoba, Huelva, Granada, Jaén, Jerez de la Frontera, Málaga and Sevilla) correspond to the most populated cities of each province of the autonomous community (see Fig. 9 for a map indicating the location of these cities). Initial values for the infections and deaths were provided by the Andalusian Government ("La Junta de Andalucía") [58]. We defined an infection radius of 10 km around the center of each hotspot, within which we placed the source of infection to initialize the epidemic, what we refer to as "blobs" of infection. These epicenters of infection were modeled via Gaussian profiles whose spatial (variance) scale was selected to be the infection radius; their amplitude was chosen such that the total number of infections, deaths, recoveries and hospitalizations, as calculated via the surface integrals of the associated densities through the region, be the same as the one reported in the original data. The population of asymptomatics and exposed was, similarly to the ODE optimization, selected to be proportional to the infected one with the proportionality ratios A 0 /I 0 and E 0 /I 0 maintained as those of the ODE.

With all these choices, the PDE model was run without optimizing at the PDE level the median parameters that are not expected to depend on spatial scales. The quantity that we varied was the diffusion-coefficient reduction factor η D in steps of 0.05 in the interval [0.25, 0.50] (six simulations in total). For each simulation the scaling factor ξ was determined by minimizing the 2 norm, as previously discussed. We used the second scenario parameters (t q = 16) since this choice reproduced better the data and the flattening of the epidemic curves with the imposition of the lockdown. The comparison of the spatially-integrated PDE results to the data for Andalusia and the ODE prediction, shown in Fig. 4 is quite promising. We show the observed data as black dots, the 0D-model predictions as the solid blue line, and the spatially integrated results of the PDE model as the shaded region. The bottom and top of the shaded region are enclosed by the curves corresponding to η D = 0.50 and η D = 0.25 (when the optimal scaling factor ξ is used). We found that both C(t) and D(t) asymptote to lower values for the optimal scaling factor as η D increases. As in the case of the 0D model, the change in the transmission rates and diffusion coefficients as a result of the lockdown leads to a flattening of the epidemic curves. An additional remark on the reported and predicted number of cases, and on the counting of asymptomatics, is in order (Fig. 4 , top left panel). As reported in [61] starting on April 13, 2020 (t = 32) the reported number of cases includes asymptomatic individuals, i.e., susceptibles that have tested positively to the presence of the virus. This might explain the underprediction of cases after t = 32.

Clearly, the spatial model can do an adequate job in capturing both the cumulative infections and the number of deaths (with the caveats to be given in the discussion below). Notice that in the bottom row of the figure, we illustrate the surface integrals of each of the density of E, A, H, R and AR as a function of time, representing the evolution of the pandemic at the "integrated" level of the entire country in an illustration similar to the one that we typically obtain from the ODE models.

In addition, we complemented the spatially-averaged results by the space-time evolution simulations of Figs. 5-7. For the reported cases in the figures, the scaling factor multiplying the β's is ξ = 0.00480, corresponding to a characteristic scale l ≈ 0.200km. A comment on the scaling factor ξ is in order. For the different choices of the diffusion-coefficient reduction factor η D we determined the optimal scaling factor to be in the range ξ [0.00480, 0.00489], for the lower value of which we obtain the reported characteristic length scale (l = √ ξN ≈ 0.200 km). It is tempting to associate the scaling factor, and in particular ξ = 0.00480, with the inverse of the "lived" population density, the population density perceived by a randomly chosen individual [62] . According to [63] the inverse lived density for Andalusia is approximately 0.0046, remarkably close to the smaller value of the scaling-factor range. A similar observation holds for the scaling factor for Greece which, in turn, suggests that this is an important insight (and not a serendipitous occurrence) as concerns the "translation" of the 0D model coefficients into the PDE ones.

Of course, the PDE model has considerable additional information through its spatial resolution. In Figs. 5-7 we can see the spatio-temporal evolution of the infections in Andalusia (i.e., the spatial distribution at a few snapshots over time), the deaths and the cumulative infections C(x, y, t), respectively. We also produced movies of the corresponding evolution that can be found in [64] . We can observe how the biggest fractions of the infections remain in the most populated cities of Sevilla and Málaga, and that the provinces of Huelva and Almería are those with the smallest number of infections, in accordance with the actual status of the pandemic [58]. The predicted spatial distribution of the total confirmed cases C(t) may be compared to the data shown in Fig. 8 , where officially reported data are presented per municipality [65] . The comparison is favorable as both figures show that Málaga is the hardest hit municipality, followed by Sevilla. In addition, the predicted number of cases in Granada, Córdoba, Jaén and Jerez de la Frontera (in decreasing number of cases), which are lower than in the previous two cities, shows the same decreasing trend as that manifested in the reported number of cases shown in Fig. 8 . The predictions for the number of deaths in the different provinces shows a slightly different behavior from the observed data. If the fatalities spatial density is integrated over each province, the predicted final number of deaths is higher in Málaga than in any other province. From the data, the final number of deaths is almost the same in Málaga, Sevilla and Granada. This suggests that modeling human mobility in these provinces solely by diffusion, and specifically by the chosen diffusion coefficients, cannot fully account for this dispersion of the number of fatalities. This may arise, within our reaction-diffusion model, possibly due to the number and intensity of the chosen hotspots of infection or to the requirement that a diffusion coefficient of different value should have been chosen.

In this work we presented a platform for establishing a compartmental epidemiological model both at the level of ODEs (0D, no spatial dependencies) in line with numerous earlier works, as well as at the spatially distributed level of PDEs to study the spatio-temporal spreading of COVID-19. The regions of interest were the mainland of Greece and the Spanish autonomous region of Andalusia for which there has been a small number of studies. As regards Greece, there are some probabilistic [15, 66] , some network-based approaches for time-series analysis [67] , and some based on the SIR variant SEAIR [22] . Studies that focus on Spain also examine Andalusia as a case example using either probabilistic [14] or POD-based decomposition techniques [68] . Our effort has been to explore a model of the SEIR variety that incorporates some of the particular biological features of the SARS-CoV-2 virus [47] , such as its latent period, and the potential to generate a significant fraction of asymptomatic hosts, which, in turn, play a crucial role in spreading the infection. The resulting SEAIHR model involves a number of populations: Susceptible, Exposed, Asymptomatics, symptomatically Infected, Hospizalized, Recovered, and deceased. We found that for the regions of interest the model reproduces the epidemiological data that we determined to be most reliable, namely the data on the cumulative infections and especially the number of deaths. Naturally, more accurate data including also spatially resolved ones (on the spatial scale of our PDE model) would be helpful towards the improved calibration of the results offered herein.

We modeled both the early, pre-quarantine, stage of the epidemic, as well as its development at a later stage when containment measures had been enforced. The effect of quarantine on the spreading of the disease was imposed via a time-dependent (on the time scale of a day) change of the transmission rates and of the diffusivities (the latter in the PDE model). While initiating a quarantine roughly when it was imposed yields more acceptable results in Greece, in Andalusia this is less so. In fact, for both regions, simulations reproduced more accurately the observed data, and in particular they captured the "angle" indicating the curbing of the infection due to government-imposed intervention measures, if a time-lag is imposed on the application of the (instantaneous in the ODE model) quarantine set of parameters. We, thus, considered two scenarios, corresponding to different delays in imposing the quarantine: the second scenarios reproduced the reported data better. We note that this seemingly artificial time shift in imposing the lockdown reflects the fact that model parameters change over a short time scale upon the introduction of lockdown measures, whereas the effect of social intervention measures (self-quarantine, social distancing, face masks, etc) appears to arise later in the epidemiological data.

We determined model parameters via optimizing model predictions with respect to reported total number of (infected) cases and number of deaths. The optimization algorithm minimized the Euclidean distance between model predictions and observed data. For the 0D model, we performed 2000 optimizations for model parameters and (the unknown ones among the) initial conditions uniformly sampled within specified ranges. Median and interquartile ranges for model parameters were determined. Additional sensitivity analyses were performed to conclude that combinations of parameters, specifically the product of the asymptomatic transmission rate times the inverse latent period, is non-identifiable. This parametric combination is intimately connected to the product of the infected transmission rate times the inverse incubation period: in fact, our results suggest that the relation between these two products can be well approximated via a straight line of negative slope. We interpreted both the median time scales of, e.g., the conversion of exposed to asymptomatics and (symptomatically) infected, and the fraction of, e.g., hospitalized that lead to recoveries or deaths. We found them, for both countries and the second scenario, to be in reasonable agreement with current epidemiological estimates. Our median results reinforce the feature prevalent in numerous studies about the importance of asymptomatics in the transmission of the SARS-CoV-2 virus, cf. [46, 69] , a particularity of this coronavirus. The asymptomatic transmission rate β AS was found to be smaller than the symptomatically infected rate β IS , coupled to the fraction of exposed evolving to asymptomatics being larger than those evolving to symptomatically infected. In [70] it was reported that asymptomatic infectious hosts may account for up to 86% of cases, thus further supporting our prediction of their importance in the spread of the disease. We remark that of the four cases studied only one (Andalusia, scenario one, early imposition of lockdown measures) had β AS > β IS (and a lower asymptomatic to infected split), a case that did not reproduce accurately the data: in the more accurate simulations of scenario two the inequality was inverted. Once again, however, we caution the reader that issues of identifiability prevent us from assigning a particular weight to the findings about the relative size of β AS vs. β IS , other than their corroborating the central role of asymptomatics in the transmission of SARS-CoV-2.

We then utilized the median parameters in a spatially distributed, reaction-diffusion model. Here, we overcame the major challenges of formulating a mesh with the boundaries of a region within the software package COMSOL and also leveraged state-of-the-art geographical methods such as the World Pop project (for population mapping based on census data) to set up distributed simulations of the pandemic spreading in the geographical domain. We consider this computational effort a significant and necessary non-trivial step for the eventual inclusion of more realistic long-range human mobility modeling via the inclusion, e.g., of convection or other modeling of directed-motion of individual populations. We pondered on how to adapt the parameters of the ODE model to the PDE framework and argued that "onsite" (i.e., single-individual) parameters can be maintained the same. We also explained the challenge of adapting contact parameters (such as the transmission rates) to the level of the country: this process involves issues of homogeneity at the ODE level vs. substantial heterogeneity at the PDE level. We also made a first series of assumptions at the level of convection (neglected herein) and isotropic diffusion (selected as the primary mechanism for disease spreading herein) to explore the time-resolved dynamics at the country/autonomous community setting.

At the level of our distributed simulations, there exist some promising results. We were able to seed the infection at some of its key epicenters and observe it to produce infections, recoveries, deaths, etc., over the entire region. The "hotspot" seeding at various locations is an indirect attempt to model the movement of infectious individuals as is, e.g., considered by metapopulation or network models. At the cumulative level of the region, surface integrations enabled comparisons with the collected data at the regional level yielding reasonable correspondence between model results and the observed cumulative epidemiological reality. Moreover, the model appears to be promising towards capturing some of the spatial features of the infection progression: for instance, visual comparison of model predictions with reported spatially distributed data for the cumulative infections shows that the model reproduces the persistence of infections in highly populated areas, albeit with a possible time lag. At the spatial level we find that the infection persisted the longest in regions of very high population density. We believe that this effort paves the way for a distributed observation of the relevant spreading, but it also has some weaknesses, challenges, and improvements that are worth considering in future steps. As stated in the Introduction, the spatio-temporal modeling of the epidemic by reactiondiffusion PDEs, and specifically with isotropic diffusion being the dominant mechanism of spatially spreading the virus is an important first step towards developing a continuous description of disease spreading where human mobility is modeled at a fine spatial scale. This approach should be contrasted to discrete network-based metapopulation models and the length scales considered in these models. A distinct advantage of the continuous model over the meta-population one is that the former can model interaction at a finer and more extended inter-nodes scale.

It would be especially useful in the context of the present pandemic of unprecedented information flow [3] to have easily accessible temporally and spatially resolved data for the evolution of the pandemic in different regions. Such "seeding" in a distributed way (rather than the colloquial seeding at hotspots performed herein) would build into the model an accurate spatial distribution of infected population, and, hence, would be far closer to the country's pandemic evolution. Indeed, there is another challenge that is arguably even more significant. Diffusion as a mechanism for spreading a disease is traditionally associated with diseases that have specific transmission characteristics, as for example vector-borne malaria that is transmitted by mosquitoes [71] . Herein, we considered diffusion as a proxy for short-term human mobility (which via its interplay with nonlinear contact interactions provides a mechanism for spreading the disease), relegating long-range transport to the initial seeding of the infection. Yet, admittedly it is not sufficient for expanding the infection at the scale of the country as our results show, at least not via realistic spatial and temporal scales of individual mobility. In particular, it has not escaped our attention that this type of spreading does not account for the directed motion of individuals (possibly infected ones) from the city to the country, or from one city to another for pleasure or business. This is especially important for travel (and hence infection transport) at a longer spatial scale (rather than the shorter one enabled by diffusion). We note that a form of anisotropic diffusion was used to model disease spreading primarily along highways in [49] .

This suggests that some form of a probabilistic element needs to be inserted in the model. One possibility that we are exploring is the spatial distribution of the initial condition of the asymptomatic population. This may generate infections in a more spatially distributed way, leading to the spatial expansion of the pandemic throughout the country in a more consistent way with the observed data [50]. A perhaps even more significant or possibly complementary perspective worth considering is, naturally, a probabilistic one. In addition to deterministic processes like diffusion or convection (which is worth integrating in a subsequent version of the model), it seems relevant to include a probabilistic gain and loss term reminiscent of (a longrange variant of) the conservative Kawasaki dynamics [72] at the level of spins. This type of term would generate infections in a probabilistic way (possibly with a probability weighed upon the region's population density) by allowing individuals to effectively "perform trips" through the country, i.e., disappearing from one location and reappearing (within a short time scale of less than a day for the regions of interest) in another.

As also discussed in the Introduction, there are other ways by which to bypass the practicalities of the application of PDEs at the level of a country. One of the canonical ones involves the application of the theory of networks in the realm of metapopulation models in a way similar to the work of [41] . Such approaches are already being brought to bear, as in the work of [45] or [14] and are certainly also worth expanding upon and refining, as well as comparing with the data available in the context of the SARS-CoV-2 virus. Building such networks for the examples of Greece and Andalusia considered herein (and of course beyond) also constitutes a worthwhile direction of future research. Clearly, further efforts at the level of data collection and curation, at the level of model setup and validation, and then at the level of optimization and utilization for prediction are needed. Our hope, however, is that the approach proposed herein is an initial step towards putting together a number of relevant tools to enable going beyond the 0D approach of ODE models and gradually considering in more detail the expansion of a pandemic at a combined spatial and temporal level.

Appendix A: Next-generation calculation of the basic reproduction number R0

We will use the next generation matrix approach of the system of equations Eqs. (1)-(8) without the spatial term to find R 0 . In particular, we set up the vectors:

The idea is that we rearrange the compartments so that the infectious/infected compartments E, A, I, H, appear first. We then place S, AR, R, D. If we calculate F − V, it should yield a reordered version of the vector field that describes our disease system. We then focus on the 4 infectious/infected compartments and ignore the rest. We find the Jacobians of F, V with respect to E, A, I, H in the order in which they appear. This will yield two 4 × 4 matrices:

The basic reproductive number is the spectral radius of F V −1 which in our case is

This result is in accordance with epidemiological intuition: the first contribution to R 0 is proportional to β SA and S * , namely, the transmission rate and total susceptible population S * . It is also proportional to σ A /(σ A + σ I ), namely the fraction of exposed hosts becoming infectious, yet asymptomatic, A. Finally, it is inversely proportional to the loss rate M AR of the infectious asymptomatic class A. The second contribution to R 0 is analogous to the first one and stems from the second mode of transmission, i.e., through contact with I. We present the "0D" model predictions for Greece in Fig. 10 . The data we consider [50] start on March 12, 2020 when losses of life started to occur and the cumulative number of infected (total number of cases) was already a bit over 100 individuals. We follow the evolution of the pandemic till May 11, 2020 using official data up to this point to optimize model parameters and initial conditions. As in the case of Andalusia, we minimized the combined Euclidean distance of model results and observed data for the total number of cases (C(t) = I(t) + H(t) + R(t) + D(t)) and the total number of deceased (D(t)), see Eq. (9) . A few remarks on the quality of the data and our choice of the most reliable time series (total number of reported cases and fatalities) are presented in Section II, following Eq. (9). We performed 2,000 optimizations to obtain the parameters shown in Table III . We present the medians of all model parameters and initial conditions, as well as the interquartile range and the range of variation of the initial values of parameters.

The model-predicted evolution of the pandemic shown in Fig. 10 was calculated for the median parameters (solid blue line) and for the cloud of the parameter variations represented as the shaded region in the figure. Officially reported data are denoted by black dots. As in the case of Andalusia, optimal parameters are shown for two scenarios, see Table II for the events time sequence. The first scenario considers that quarantine was strictly enforced at March 24, 2020 (t q = 13). In reality, the lockdown in Greece started at the end of March 22, 2020;  it is reasonable to assume that it was strictly enforced 1-2 days later. In the second scenario the lockdown was considered to have been imposed on April 3, 2020 (t q = 23). To account for the change in parameters due to the lockdown, we imposed the time dependence of the transmission rates β as shown in Eq. (10) in the main text. This In all panels, shaded regions correspond to the interquartile range for each quantity, whereas the full line corresponds to simulations with the median parameter (and initial-condition) values. time dependence forces the transmission rates β IS and β AS to decrease by a factor η IS and η AS , respectively, relatively abruptly at the time the lockdown was imposed, t q .

The top panels of Fig. 10 compare model predictions for the two scenarios for the total number of cases (left) and fatalities (right). Scenario-one median parameters, i.e., imposing the quarantine practically at the time when it was officially announced, capture the data for fatalities in Greece fairly well. However, the number of reported cases, top left panel, is not that accurately reproduced. We attribute this discrepancy to the previously mentioned characteristic feature of the data (top left panel) that model predictions fail to capture adequately: the "angle" in the semi-logarithmic plot associated with the curbing of the cumulative number of infections C(t) due to containment measures. As the optimization algorithm attempts to minimize the distance from the observed data, initially it slightly over-predicts and then under-predicts the data and eventually the long-term predictions seem to over-predict the flattening of the cases curve. Nevertheless the overall differences are relatively small: the model prediction flattening out (over 5 months) around 200 deceased and slight over 3K infected individuals seem reasonable, were the lockdown measures potentially extendable to such a long time interval. If, as in the case of Andalusia, were we to shift the time of the application of the quarantine date by about 10 days later (scenario two), then we note in the top panel of the figure a nontrivial difference. Most notably, without significantly missing on D(t) we capture accurately the angle in the C(t) data. The relevant parameters (medians, interquartile range) are presented in the second column of Table III . As before, we justify our decision to consider a second scenario in that lockdown measures have an almost immediate effect in model predictions, whereas in reality there is a time lag before restrictive measures have a measurable effect. Lastly, we note that the effect of the shift of t q is far less TABLE III. ODE parameters for Greece: optimal (best-fitting), median and interquartile range, and variation range used in the optimization algorithm. Initial parameters and initial-condition guesses were uniformly sampled within these ranges.

Initial value (tq = 13) (tq = 23) severe than in the case of Andalusia (Fig. 2) . It is worthwhile to compare the scenario-two median parameters (Table III, second column) to those we found for Andalusia. There are no particularly noteworthy differences, although some do exist. For example, as in the case of Andalusia, the median incubation period is approximately 4 days (3.72), the latent period approximately 3 days (2.89), the asymptomatic infectious period about 7 days (6.87) and that of the infected 6 days (6.13). A slight difference is noted in the fraction of infected that need to be hospitalized (γ ≈ 0.44 instead of 0.55), and the fraction of hospitalized that become fatalities (ω ≈ 0.21 instead of 0.25). As in the case of Andalusia, we find β IS > β AS , but note the proviso related to parameter identifiability reported later on, and that the fraction of exposed who turn asymptomatic is approximately 0.56 while the ratio of asymptomatics to symptomatically infected is 1.29. After lockdown measures are imposed, the two transmission rates decrease by the same amount (η IS = η AS ), again as we found for Andalusia.

The calculated basic reproduction number R 0 (see Appendix A) for the scenario-one pre-quarantine period (t q = 13), i.e. calculated with data from the first column of Table III with the η set to unity, is R 0 = 1.64 (1.60−1.68), median and interquartile range. The effective reproduction number, i.e., the reproduction number at the beginning of the quarantine with the associated change of the transmission rates (η's as reported in column one) was calculated to be R eff = 0.849(0.842 − 0.854), reflecting the curbing of the epidemic curves, as shown in Fig. 10 . Similarly, the calculated pre-quarantine basic reproduction number for scenario two (t q = 23) is R 0 = 1.32 (1.31 − 1.34). This post-quarantine effective reproduction number decreases to R eff = 0.64 (0.63 − 0.65), an indication that intervention measures lead to a curbing of the epidemic.

We conclude this section by a brief discussion of parameter identifiability and the zero eigenvalues of the Hessian of the variation of the Euclidean norm with respect to model parameters, in the spirit of the corresponding discussion for the case of Andalusia. Figure 11 presents our calculations relevant to parameter identifiability and sensitivity to parameter variations (for the more reliable simulations of the second scenario). As for Andalusia, we performed 200 optimization with parameters uniformly sampled within their variation range with the ratio β AS /β IS fixed at one of ten equidistant values chosen within [0. 2, 2] . The top panel in the figure shows the inverse (apparently nearly linear) relationship between properties of asymptomatics (β AS σ A ) and the corresponding properties of infected, a relationship that we argued may be interpreted as the manifestation of the singular eigendirection of the Hessian. As the fraction of asymptomatics to infected increases the asymptomatic transmission rate decreases. Lastly, we note that the cloud of points for the total number of cases slightly overpredicts the data initially en route to eventually (slightly) underpredicting their long-term evolution. The cloud of points follows rather closely the data for total fatalities, with the data eventually lying near the bottom edge of the prediction interval. . While the transmission rates (times the associated time scales) vary significantly, the forward, model-predicted evolution of the pandemic does not. This is an indication that even though the model parameters are not identifiable, they provide an adequate predictor of the evolution of the pandemic if chosen within a suitable range.

Results relevant to the scenario-two PDE simulations for the mainland of Greece may be found in Fig. 12 for the same parameters as for the 0D model, and in Figs. 13-15 for the spatio-temporal simulations of the pandemic.

The initialization of the model for Greece consisted of selecting five of the key "hotspots" of the infection, as they arose in Greece, to populate initially. We defined an infection radius of 10 km around the center of Athens (largest city and capital), Thessaloniki (second largest city and source of the first infection), Patras (third largest city and the location where a key imported group of infected individuals was transferred), as well as Kozani and Xanthi. These cities where the epidemic was initialized are identified in the 2000 population density map shown in the right panel of Fig. 16 . The latter two are two significant peripheral centers where infections were seeded early on. In Athens, we placed the largest (by a factor of two) source of infection, while similar "blobs" of infection were initialized in the remaining four cities. As in the case of Andalusia, these epicenters of infection were initiated via Gaussian profiles whose spatial (variance) scale was selected as the infection radius; their amplitude was chosen so that the total number of infections, deaths, recoveries and hospitalizations, as calculated via the surface integrals of the associated densities through the country, be the same as the one reported in the original data. The population of asymptomatics and exposed was, similarly to the ODE optimization, selected to be proportional to the infected one with the proportionality ratios A 0 /I 0 and E 0 /I 0 maintained as those of the ODE.

As for Andalusia, having initialized the PDE model, it was run without optimizing the median parameters that are not expected . The plots displayed on the left panels correspond to tq = 13 (March 22, 2020) and those on the right panels hold for tq = 23 (April 3, 2020). The solid blue line (top two rows) reproduces the 0D simulation, cf. Fig. 10 . The median parameters shown in Table III are used, except for the transmission rates that were scaled by ξ [0.00216, 0.00234] for tq = 13 and ξ [0.00245, 0.00248] for tq = 23. The optimal scaling factor ξ increases as the diffusion-coefficient reduction factor ηD increases. The latter was varied in the interval [0.25, 0.50] in steps of 0.05. Shaded regions are delimited by the optimal plots for ηD = 0.5 and ηD = 0.25 ; their order at t = 150 is the same as that of the bottom panels of Fig. 10. to depend on spatial scales at the PDE level. Instead, we varied the diffusion-coefficient reduction factor η D in steps of 0.05 in the interval [0.25, 0.50] (six simulations in total). For each simulation the scaling factor was determined by minimizing the 2 norm, as previously discussed. Similarly to the 0D simulations, the comparison of the spatially integrated PDE results to the data for Greece (see Fig. 12 ) is not particularly good for the first scenario (t q = 13, left column) that considers only a minor shift of the quarantine time. In fact, the 0D results do not fall within the range of the spatially averaged PDE simulations. Further modification of the quarantine time (scenario two) can also help capture once again the "angle" in the relevant data (right column). Clearly, the spatial model, via the surface integral of the population densities, does a very adequate job at capturing both the cumulative infections and the number of deaths. Notice that in the bottom row of the figure, we illustrate the surface integrals of each of the densities of E, A, H, R and AR as a function of time, representing the evolution of the pandemic at the "integrated" level of the entire country in an illustration similar to the one that we typically obtain from the ODE models. That being said, of course, the PDE model provides considerable additional information through its spatial resolution. In Figs. 13-15, we can see the spatio-temporal evolution of infections (i.e., the spatial distribution at a few snapshots over time), fatalities and cumulative infections C(x, y, t), respectively. These figures were generated using the median parameters of the second scenario, and with the transmission rates multiplied by the scaling factor ξ = 0.00246. The corresponding length scale is l = √ ξN ≈ 0.163 km, comparable to what we found for Andalusia l ≈ 0.200km. In addition and as in the case of Andalusia, if we associate the scaling factor with the inverse of the "lived" population density we find that in Greece [73] it is around 1/379 ≈ 0.0026, remarkably close to the numerically-determined ξ. We believe that the identification of the scaling factor as very closely matching the inverse of the lived population density in two entirely independent (and quite distinct in their number of infections) cases suggests this scaling as a nontrivial insight stemming from these studies about the connection of ODE and PDE models.

We also produced movies of the corresponding evolution that can be found in [64] . It is important to re-iterate here that we have not included the islands of Greece in this effort (i.e., we are looking at the mainland of Greece). Obviously if one were to model the disease spread in each of these islands it would be relevant to seed the infection in each island individually and study the spreading there rather than together with the spatially disconnected from the islands mainland of Greece.

We can clearly see how the infection spreads throughout the country, affecting most significantly the regions of higher pop- ulation density. Indeed, it is clear that over time the infection is extinguished in most regions and it finally persists chiefly in Attica, the region of highest population density (and where the main metropolitan center of the country, Athens, lies); see, in particular, the bottom panels of Fig. 13 . Nevertheless, it is evident from Fig. 14 that a number of deaths develops in each of the 5 regions where the infection was initially seeded, in line with the corresponding expectation from the country's data. Indeed, also, each region features a discernible fraction of cumulative infections in Fig. 15 , although clearly once again the lion's share of infections pertains to Attica. The second biggest fraction of infections pertains to Thessaloniki (the second biggest metropolitan center) and so on. It is clear from these figures that the model yields a reasonable prediction of the spatial evolution of the disease spread, in line with the cumulative totals of deaths and infections throughout the country. Nevertheless, a comparison with the spatial distribution of the pandemic throughout the country [50], see Fig. 16 , suggests also some limitations. The spatial snapshot of the number of confirmed cases per prefecture shown in that figure was created on March 29, 2020, corresponding to t = 18. Hence the data are intermediate between the two bottom panels of Fig. 15 , at t = 11 and 21 days. The persistence of the epidemic in Attica (dark red) and to a lesser degree Thessaloniki (red) is evident in both figures. The other areas seem to be slightly underpredicted. The comparison suggests a partial time mismatch between the simulations and the data. One can also observe that since the initial spots for Andalusia were spatially close, the infection spreads more easily in that region than in the mainland of Greece. 

SARS outbreaks in Ontario, Hong Kong and Singapore: the role of diagnosis and isolation as a control mechanism

Interhuman transmissibility of Middle East respiratory syndrome coronavirus: estimation of pandemic risk

Novel Vaccine Technologies: Essential Components of an Adequate Response to Emerging Viral Diseases

World Health Organization Writing Group

Impact of non-pharmaceutical interventions (NPIs) to reduceCOVID-19 mortality and healthcare demand

What aerosol physics tells us about airborne pathogen trannsmission

Estimating the asymptomatic proportion of coronavirus disease 2019 (COVID-19) cases on board the Diamond Princess cruise ship

COVID-19 Outbreak Associated with Air Conditioning in Restaurant

Feasibility study of mitigation and suppression strategies for controlling COVID19 outbreaks in London and Wuhan

Projected Development of COVID-19 in Louisiana

A spatial model of CoVID-19 transmission in England and Wales: early spread and peak timing

Modeling the Spatiotemporal Epidemic Spreading of COVID-19 and the Impact of Mobility and Social Distancing Interventions

Estimating the infection horizon of COVID-19 in eight countries with a data-driven approach

A spatial-temporal model for the evolution of the COVID-19 pandemic in Spain including mobility

Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe

Modeling and forecasting the early evolution of the Covid-19 pandemic in Brazil

Model studies on the COVID-19 pandemic in Sweden

Modeling Covid-19 dynamics for real-time estimates and projections: an application to Albanian data

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

On the transmission dynamics of SARS-CoV-2 in a temperate climate

A mathematical model for the novel coronavirus epidemic in Wuhan, China

Estimation of the transmission risk of the 2019-nCoV and its implication for public health interventions

Tracing day-zero and forecasting the COVID-19 outbreak in Lombardy, Italy: A compartmental modelling and numerical optimization approach

Mathematical modeling and the transmission dynamics in predicting the Covid-19-What next in combating the pandemic

Modelling insights into the COVID-19 pandemic

Wrong but useful-What COVID-19 epidemiological models can and cannot tell us

Understanding variants of SARS-CoV-2

The mathematics of infectious diseases

The mathematical theory of epidemics

Mathematical structures of epidemic systems

Localized outbreaks in an S-I-R model with diffusion

A two-phase epidemic driven by diffusion

Mathematical biology II: Spatial models and biomedical applications

Modeling infectious diseases in humans and animals

Demography and Diffusion in Epidemics: Malaria and Black Death Spread

Dynamics of infectious disease transmission by inhalable respiratory droplets

Geographic and temporal development of plagues

Universal model of individual and populatrion mobility on diverse spatial scales

A reaction-diffusion system to better comprehend the unlockdown: Application of SEIR-type model with diffusion to the spatial spread of COVID-19 in France

Simulating the spread of COVID-19 via spatially-resolved susceptible -exposed -infected -recovered -deceased (SEIRD) model with heterogeneous diffusion

Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: Theory and simulations

Estimating spatial coupling in epidemiological systems: A mechanistic approach

Modeling human mobility responses to the large-scale spreading of infectious diseases

Equation-Free, Coarse-Grained Multiscale Computation: Enabling Mocroscopic Simulators to Perform System-Level Analysis Comm

Parallel algorithm for simulating the spatial transmission of influenza in EpiGraph

Presymptomatic SARS-CoV-2 infections and transmission in a skilled nursing facility

SARS-CoV-2 (COVID-19) by the numbers

Modeling the transmission of Middle East Respiratory Syndrome corona virus in the Republic of Korea

Propagation of epidemics along lines of fast diffusion

Disaggregating census data for population mapping using random forests with remotely-sensed and ancillary data

Neural network aided quarantine control model estimation of global Covid-19 spread

Accounting for Symptomatic and Asymptomatic in a SEIR-type model of COVID-19

Lockdown Measures and their Impact on Single-and Two-age-structured Epidemic Model for the COVID-19 Outbreak in Mexico

Parameter set uniqueness and confidence limits in model identification of insulin transport models from simulation data

The COVID-19 pandemic as experienced by the individual

URL will be inserted by publisher] for animations on the time evolution of the predictions of infections, deaths and total cases for Andalusia and Greece

Estimating the number of infections and the impact of non-pharmaceutical interventions on COVID-19 in European countries: technical description update

The effect of anti-COVID-19 policies on the evolution of the disease: A complex network analysis of the successful case of Greece

Predictive data assimilation through reduced order modeling for epidemics with data uncertainty

Temporal dynamics in viral shedding and transmissibility of COVID-19

Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2)

A reaction-diffusion malaria model with incubation period in the vector population

Acknowledgments. The authors are indebted to Dr. Maksym Bondarenko for his substantial help with the WorldPop maps and setup, and also thank Dr. Jinlan Huang for assistance in setting the relevant computation up in COMSOL. In addition, we thank Francisco Rodríguez Sánchez for providing the data and R code used to generate the map with the spatial distribution of the epidemic in Andalusia, Fig. 8 . PGK gratefully acknowledges discussions and input from Andy Ludu including regarding the layout of Fig. 1 . This material is based upon work supported by the US National Science Foundation under Grants No. DMS-1815764 (ZR), PHY-1602994, and DMS-1809074 (PGK). J.C-M. also thanks the Regional Government of Andalusia under the project P18-RT-3480 and MICINN, AEI and EU (FEDER program) under the project PID2019-110430GB-C21. ZR and PGK also acknowledge support through the C3.ai Digital Transformation Institute. PGK also acknowledges support from the Leverhulme Trust via a Visiting Fellowship and thanks the Mathematical Institute of the University of Oxford for its hospitality during this work.Disclaimer The views expressed in this manuscript are purely those of the authors and may not, under any circumstances, be regarded as an official position of the European Commission.