Challenges in the mathematical modeling of the spatial diffusion of SARS-CoV-2 in Chile

The COVID-19 pandemic, which started in 2019, has become one of the most severe public health crises in recorded history. As of June 2024, according to the World Health Organization (WHO), there have been a total of 776 million reported cases and more than 7 million deaths around the world ^[1]. Specifically, in Chile, there have been more than five million cases and 57,918 deaths ^[2]. During the pandemic, extended shelter-in-place timelines by the Centers for Disease Control and Prevention (CDC) fostered public distrust in science and government, worsening patient outcomes and increasing fatalities as citizens disregarded CDC regulations. Epidemic models are crucial for a better understanding of the infection rate predictions and improving prevention strategies.

A variety of works have proposed mathematical models for the dynamics of the COVID-19 pandemic ^{[3,4,5,6,7,8,9,10]}, and these are often based on ordinary differential equations (ODEs). However, there are other works that have used difference equation to model the COVID-19 pandemic ^[11]. For instance, in ^[11] a difference equations model was used to study the COVID-19 pandemic in Santiago de Cali from March 2020 to June 2022.

The most common and simplest approach is to compartmentalize the population into susceptible, infected, and recovered individuals (i.e., SIR framework) ^[12,13]. However, the population can be further subdivided according to the progression of the COVID-19 disease, the status of vaccination, or the age ^[14,15]. The increased number of compartments introduces added complexity to the model, often enhancing its capacity to replicate real-world phenomena, albeit at the cost of incorporating additional state variables and parameters. From a practical point of view, this can raise identifiability issues related to the parameters, thereby raising uncertainty ^[4,16,17]. Fred Brauer was an eminent mathematician who made important contributions to mathematical epidemiology, which is very related to real data. He recognized that calibrating mathematical models to data is oftentimes necessary in order to use infectious disease transmission models in real public health situations ^[18].

It is important to note that the proportion of the previously mentioned models that deal with real-world data is relatively small ^[15,19,20]. When working with real-world data, many difficulties arise, in part due to the necessity to give numerical values to the model parameters ^[21]. This often encourages researchers to avoid working with data when proposing mathematical models for disease dynamics. Moreover, the computational challenges increase when dealing with spatial-temporal models, which are often based on partial differential equations (PDEs) ^{[22,23,24,25,26,27]}. For example, ^[22] addresses the introduction of a new strain with a structural change at time $t_1$ and proposes a spatial-temporal model for population mobility. Similarly, ^[25] presents a SIR framework with intra-city and inter-city mobility to study COVID-19. In ^[28] a mathematical model to investigate the spatial spread of an infectious disease with a moving boundary was presented.

The first COVID-19 reported case in Chile occurred on March 3, 2020. The peak of the first epidemic wave was reached on June 14, 2020, with $6938$ new cases. The peak daily death rate was reached on June 13, 2020. The first epidemic wave was primarily focused on Chile's capital. By July 30, 2020, $9457$ people had died in Chile ^[2]. In this paper, we work with real data on reported infected cases. It should be noted that earlier research has examined COVID-19 data from Chile through different methodologies and found engaging results ^{[29,30,31,32,33]}. For example, ^[32] estimated trend changes, showing that non-pharmaceutical interventions (NPIs) before March 2021 impacted case numbers in most regions. In ^[34] ARIMA models were used to predict confirmed cases from March 2, 2020 to July 14, 2020, assessing accuracy with training and test sets. ^[30] proposed a multi-group SEIRA model, studying NPI effectiveness by estimating the effective reproduction number's evolution. In ^[35] a multi-parameter logistic growth model was fitted to early phase data, predicting pandemic growth until October 2020. In ^[33] the authors estimated the reproduction number at $1.8$, noting a shift to linear growth post-lockdowns on July 7, 2020. In ^[31] the authors used a generalized SEIRQ model to predict COVID-19's course, showing differences across three scenarios. In ^[36] the authors fitted a generalized non-autonomous SIR model to early-phase data from Antofagasta, Metropolitan, and Nuble regions, using constant delays and varying transmission rates. In ^[37] the authors studied case fatality risk by age and gender, finding the highest mortality rates among male seniors ($\geq70$ years).

Notable works based on PDEs have also been introduced to model different infectious diseases, but few studies have dealt with the COVID-19 pandemic ^{[38,39,40,41,42]}. Even fewer of these works included real-world data and the calibration of the mathematical models used ^[8,42,43]. In ^[38], a spatial SIR model based on a system of PDEs that can be used for many infectious diseases was introduced. Sufficient conditions for local and global stability were found in terms of the basic reproduction number $\mathcal{R}_0$. The effect of the diffusion and transmission rates was studied under Neumann boundary conditions. In ^[43], a spatio-temporal model based on PDEs was used to quantify compliance with the US COVID-19 mitigation policies at a regional level and was validated with short-term predictions. In ^[10], a COVID-19 model based on a system of PDEs is proposed to study how different vaccination-isolation strategies affect the COVID-19 pandemic. The authors studied the system's asymptotic distributions of endemic equilibrium under different conditions, and the formula for the basic reproduction number was found. The proposed model includes Neumann boundary conditions and more than 15 parameters. Although the study is interesting, having a large number of parameters makes it very difficult to calibrate the model with real data. In ^[44] a spatial-temporal model to study the dynamics of COVID-19 pandemic was presented. Theoretical results regarding positivity, boundedness, and stability were obtained. Numerical simulations that show the effect of isolation and social distancing rate were presented. The authors used only the temporal COVID-19 model to fit the real data. In ^[42], a spatial model for the spread of the coronavirus in the Canadian province of Nova Scotia was presented. Several generalizations of the model were proposed. In ^[45] graph theory was used to analyze spatial COVID-19 transmission dynamics between municipalities in Aveiro district, Portugal, and between different age groups. There are other works that have used spatial-temporal models to study different phases of the COVID-19 pandemic or even within-host dynamics of SARS-CoV-2 ^{[46,47,48,49,50]}. In ^[51] a review of spatio-temporal models for COVID-19 pandemic is presented and it was found that oftentimes the modeling results are not accurate due to the unavailability of detailed data.

We propose a mathematical framework to analyze the early dynamics of COVID-19 in Chile. This framework uses spatial-temporal models based on systems of PDEs to incorporate both the spread of SARS-CoV-2 and population mobility. The model includes diffusion effects to represent population movement and categorizes the population into susceptible, infected, and recovered groups. By employing a reaction-diffusion model of the SIR epidemiological type, we can account for various transmission rates across different regions and times. The proposed model aims to describe the early situation where individuals exhibited significant fear of getting infected with SARS-CoV-2 infection and consequently avoided high-risk regions, such as densely populated locations ^[52,53,54].

In this study, we approach the modelization of the spatial-temporal COVID-19 dynamics in Chile using different methodological approaches for the spread of SARS-CoV-2. First, we employ a constant SARS-CoV-2 transmission rate across all of Chile. Second, we proceed similarly but allow the SARS-CoV-2 transmission rate to vary over the space. This latter approach seems more realistic since this rate might differ regionally. In our work, we account for a different diffusion rate for infected individuals since their mobility is likely limited by the disease itself and by self- and government-imposed measures to prevent the spread of SARS-CoV-2.

Performing a calibration of the spatial-temporal model with real data from Chile is a very challenging task ^[55] even if we do not consider spatial effects ^{[25,31,34,36]}. An alternative approach that can be used to approximate and characterize the dynamics of the COVID-19 pandemic in Chile is to calibrate a mathematical model based on ODEs (no spatial effects) to the data of each region ^{[25,31,34,36]}. For the particular case of Chile, there are 16 regions, which would require a total of 16 calibrations. Each region would have its own SARS-CoV-2 transmission rate without taking into account the mobility of individuals in Chile. One might assume that calibrating one spatial-temporal model would be easier than 16 calibrations. However, it is a trade-off since the spatial-temporal model has more parameters and more variability that could give rise to more uncertainty ^{[10,20,23,25]}. In summary, in this work we construct models and present the challenges related to modeling the COVID-19 pandemic with spatial-temporal real data with a mathematical model based on systems of PDEs. These are significant contributions that can help with future pandemics and also provide a deeper understanding of what happened during the COVID-19 pandemic. In ^[56] an interesting theoretical work to study the COVID-19 epidemic was presented. The authors used a classical SIR model with diffusion, and an additional compartment composed by the infected individuals traveling on a line of fast diffusion (roads).

This paper is organized as follows: In Section 1, we present the spatial-temporal model, the main assumptions related to this model, and the procedure for the calibration process. In Section 2, the main aspects regarding the mathematical modeling approaches, calibration process, initial and boundary conditions are presented. In Section 3, we present the results for different modeling approaches and the calibration of the mathematical model to the real spatial-temporal data. In Sections 4–6, we present summary, discussion, and the main conclusions, respectively.

1.   Mathematical modeling approaches for the early phase of COVID-19 pandemic in Chile

In this section, we present the spatial epidemiological mathematical model for the dynamics of the early COVID-19 pandemics in Chile. The model is designed to include the diffusion effects of SARS-CoV-2 in Chile. The compartmental model divides the total population based on the COVID-19 disease progression status. In particular, the model considers susceptible, infected, and recovered individuals. The model describes the dynamics of these subpopulations with regard to space and time. A system of first-order partial differential equations is used for the mathematical model in order to include spatial and temporal effects.

1.1.   Spatial-temporal data of the COVID-19 pandemic in Chile

Chile's territory shape resembles a narrow strip between the Pacific Ocean (west) and the Andes mountains (east). Chile's length from north to south is 4270 km, and from west to east, it measures an average of only 177 km. Moreover, most of the east side has a very low population density due to the Andes mountains ^[57,58]. Chile's narrowest section is just approximately 66 km wide. It has been considered one of the most unconventional territory shapes in the world ^[57,58]. The population of Chile has been estimated as approximately 18.5 million for the year 2020 ^[57]. Based on these previous facts in this study, we assume a spatial mathematical model that considers one spatial variable that relates to the latitude position within Chile.

Daily COVID-19 data was collected from March 10, 2020, through September 15, 2020, for the country of Chile. The data is separated for each region in Chile. We sort these data into new infections per week in order to smooth it, which is a standard procedure in epidemiology.

Figure 1 illustrates the COVID-19 infected cases data in Chile separated by regions, from March 10, 2020, to September 15, 2020. Note that here, for simplicity, all regions are represented with the same length. In the Appendix, we present Table 2 which shows the inhabitants of each of the regions in Chile and their respective identification numbers. The inhabitants and density of each region were collected from ^[59]. The lengths of each region are ad hoc approximations. Notice that the Metropolitan region (region 7) is the most populated. It also has the highest density of all of Chile's regions.

1.2.   Spatial-temporal mathematical model for COVID-19 pandemic in Chile

In this subsection, we present a spatial-temporal mathematical model that describes the dynamics of the early spread of SARS-CoV-2 in Chile. In particular, we use a reaction-diffusion model of SIR epidemiological type. The model attempts to describe the early situation in which people were very afraid of getting infected with SARS-CoV-2 and moved away from high-risk regions, such as crowded places ^[52,53,54]. In many countries it was observed that some people moved from urban areas to rural areas in order to be in a lower risk place with regard to SARS-CoV-2. In ^[53] it was found a higher than average mobility from high density population regions to low density regions. The factors that people moved from urban to rural areas are various. As mentioned above, one main factor was to decrease the likelihood of getting infected, but also economic factors since many workers living in urban areas were not able to continue working or their revenues decreased ^[4,52,54,60]. Let $\Omega$ be a bounded domain in $\mathbb{R}$ with smooth boundary $\partial \Omega$. Then, based on the previous description and the mechanisms of the spread of SARS-CoV-2 in the population, we present a spatial-temporal mathematical model as follows:

The initial condition can be written as

The boundary condition is given by

This condition implies that there was no mobility between Chile and Peru due to the lockdown restrictions during the early COVID-19 pandemic ^[32]. The most southern region of Chile is delimited in the south part by the South Pacific Ocean, so people cannot move south of that region. Thus, the selected boundary condition (1.3) is plausible.

The uniqueness and continuity of the solutions are guaranteed if the initial values are positive and the functions for the diffusion and transmission in the system (2) are assumed to be sufficiently smooth in the domain $\Omega$ ^[38].

The state variables, parameters, and their descriptions are listed in Table 1. In the mathematical model the density is given in terms of population/km, but can be changed without loss of generality. Notice that the mathematical model (1.1) includes self-diffusion terms that represent the movements of individuals, that are modulated by the population densities over the space $x$. Also, the PDE model assumes constant population size, which is well justified for short dynamics ^[12]. The importance of the parameters and initial conditions of the model (1.1) is much more relevant.

1.3.   Diffusion coefficients

In this section, we present some modeling approaches for the variation of diffusion coefficients. There are many different forms that can be used for the diffusion coefficients $d_i$ of the mathematical model (1.1). Some of the options that we can set for the diffusion coefficients are the following:

● Constant and equal: In this case, all diffusion coefficients have the same value ^[38]. This form seems unrealistic, since we can expect that some part of the infected population (symptomatic) should not travel or would reduce mobility ^[10].

● Constant and different: In this case, all diffusion coefficients are constant, but some of them have different values to take into account the low mobility of symptomatic people ^[10].

● Space-varying diffusion coefficients $d_i(x)$. This form can be suitable for scenarios where travel restrictions in each region are different or when due to social behavior the mobility is different.

● Time-varying diffusion coefficients $d_i(t)$. This form can be suitable for scenarios where travel restrictions vary over time depending on awareness of the pandemic or NPIs ^{[32,61,62,63]}.

● Space- and time-varying diffusion coefficients $d_i(x, t)$. This form is a combination of the two previous ones, where the diffusion coefficients vary depending on the space and time. This form is more complex, since it has more variety in the potential forms of the functions $d_i(x, t)$. This form has more flexibility and might be more realistic.

It has been mentioned that when introducing spatial effects, the virus often propagates from its origin and if the mobility is low, the propagation occurs in a wavelike manner ^[42].

1.4.   Numerical solution of the spatial-temporal model

The numerical solution is obtained using the pdepe() built-in Matlab function ^[64]. In ^[65] the authors compared the reliability of the built-in Matlab function pdepe() by solving one-dimensional steady-state advection dispersion reaction equations. The pdepe() solver was shown to be reliable and can be used for teaching or in modeling water quality. The pdepe() built-in Matlab function has been used to solve mathematical models of brown stock washing problems where a system of governing partial differential equations arises ^[66].

The pdepe() function is designed to solve systems of parabolic and elliptic PDEs in one space dimension and time. This built-in Matlab function belongs to the Partial Differential Equation Toolbox from Matlab and it is based on the algorithm developed in ^[67]. This method uses the Galerkin and Petrov-Galerkin methods depending on the singularity of the PDEs problem. The algorithm discretizes the space and then solves an explicit system of ODEs that results from the discretization. The algorithm can use a variety of integrators in time, but the pdepe() function uses the ode15s built-in Matlab function. This ode15s function is suitable for numerically solving differential-algebraic equations that emerge from elliptic PDEs ^[67]. The ode15s() function uses a variable-step and it is a variable-order solver based on the numerical differentiation schemes of orders 1 to 5 ^[68,69]. The pdepe(), performs well with initial conditions that correspond to parabolic PDEs. It is important to note that in all the numerical simulations performed in this study, we deal with parabolic PDEs.

1.5.   An initial numerical simulation of the spatial-temporal model

First, we use the pdede solver to obtain an initial numerical solution of the spatial-temporal mathematical model (1.1). The diffusion parameters used for this first numerical simulation were chosen small in order to approximate a low mobility, $d_i = 0.91\times10^{-12}$ for $i = 1, 2, 3$. Computing a good estimate of the diffusion coefficient requires a great amount of data that is not available. However, in ^[70] the authors used $d\in (10^{-1}, 10^{-4})$ for mosquitoes which we assume that have higher mobility than people during COVID-19 pandemic. Nevertheless, the two scenarios (low and high mobility) considered in this work regarding diffusion allows us to show the effect of the diffusion. We assume an average length of infection of 1 week, which is the same as the time step of our model making $\gamma = 1$. Since we are interested in the early phase of the COVID-19 pandemic, we assume initial conditions such that most people are susceptible and the infected population is small. To begin, we use a Gaussian form to distribute infected people through the regions of Chile, assuming that there are more infected people cases in central Chile due to higher population density ^[57,71]. The initial condition for this first numerical simulation is $S(x, 0) = 1,250,000, I(x, 0) = 1/(\sigma \sqrt{2 \pi})\, exp\left(-\frac{1}{2}\left({\frac{x-8}{\sigma}}\right)^2\right)$. The first density allows to approximate the total population of Chile. The boundary conditions are of Neumann type in order to assume zero mobility between the north and south boundaries of Chile.

Figures 2 and 3 show two different preliminary numerical solutions of the spatial-temporal mathematical model (1.1) for the infected population density $I(x, t)$. In Figure 2, it can be seen that the solution is negative in some regions. However, model (1.1) cannot have negative values since it is a combination of the heat equation and the SIR ODE model. From the equations of the SIR diffusion model, it can be proven that $I(t, x)$ cannot have negative values ^[38]. Therefore, the numerical solution shown in Figure 2 is inaccurate. We have found many scenarios where negatives values are obtained for the state variable $I(x, t)$. This is especially more common when the PDE system becomes stiff, which, for example, occurs when the transmission rate $\beta$ is large. However, this can be fixed by decreasing the space and time step sizes ^[72,73]. The Courant–Friedrichs–Lewy (CFL) condition provides a necessary condition for convergence of the numerical solution ^[74]. This condition involves choosing the length of the time step size $\Delta t$ depending on the space step size $\Delta x$ and the PDE under consideration. Thus, for the reaction-diffusion model (1.1) this CFL condition is $\Delta t\le (\Delta x)^2/2d$. Using a time step size $\Delta t = 0.0001$ and a spatial step size $\Delta x = 0.001$, one obtains a computation time of approximately half an hour in a computer with the 12th Gen Intel(R) Core(TM) i7-12700 2.10 GHz and 64-bit Windows 10 Enterprise LTSC operating system. This aspect creates a computational issue related to the calibration process, which is related to solving an inverse problem related to the model (1.1). Figure 3 shows another numerical solution with smaller step sizes. Note that the solution for the component $I(x, t)$ becomes positive.

2.   Mathematical modeling approaches, resembling dynamics, initial and boundary conditions

In this section, we describe the mathematical modeling approaches and the resembling process to the real-world data related to the early phase of COVID-19 in Chile. There are many paths that can be taken to perform the resembling process. This process is very challenging if we aim to approximate the real-world scenario for each region in Chile, since NPIs were implemented regionally at different times. In addition, social behavior and population density differ regionally as well.

2.1.   Initial conditions

Let us first proceed with some options that can be used for the initial conditions for the susceptible population density $S(x, 0)$, the infected population density $I(x, 0)$ and the recovered population density $R(x, 0)$. The form of the initial condition is crucial for the calibration process, since this affects the transient dynamics and therefore the early phase of the COVID-19 pandemic ^[75,76,77]. There are many options for the initial conditions that can be used as approximations of the real-world scenario. We can choose a function $N(x)$ that represents the initial population density of Chile. Thus, one gets

Now, since at the beginning of the early phase of the COVID-19 pandemic we have that $R(x, 0) = 0$ and $I(x, 0)\approx0$. Then, one obtains that $S(x, 0)\approx N(x)$. Through this study we will use these approximations for the numerical simulations and results.

2.1.1.   Uniform initial conditions for each region

One plausible approach for the initial condition is a uniform distribution of the population over each region in Chile. In this way, each region has an accurate approximation of its respective population. The integral of the population density $N(x)$ in each region provides the approximated population in each of these regions. However, this particular approach creates a computational problem, since the initial conditions for the susceptible $S(x, t)$, infected $I(x, t)$ and recovered $R(x, t)$ populations would be discontinuous functions. From a numerical computational viewpoint the problem becomes challenging due to the jumps in the boundaries between the regions. The initial condition is given by a discontinuous piecewise function. This can be written as

where $\hat{N_{i}}$ is such that the integral of $N(x, 0)$ over the region $i$ gives the population from region $i$. Then, using Eq (2.1) and assuming $R(x, 0) = 0$, one can obtain the initial conditions $S(x, 0)$ and $I(x, 0)$, which are also discontinuous piecewise functions. Moreover, this initial condition has issues with the compatibility condition of the the boundary condition since at the boundaries the derivatives should be approximately zero.

2.1.2.   Spline interpolating polynomial for continuous initial conditions for each region

One first choice that is common when a continuous and differentiable function is required to approximate a set of points is to use spline interpolating polynomials. As we have mentioned previously, choosing the initial condition by using Eq (2.1) has the discontinuity issue. However, one can attempt to solve this aspect by using spline interpolating polynomials. Nevertheless, we have found that this modeling approach creates unrealistic oscillations in the initial conditions due to the jumps. Moreover, depending on the number of interpolating points used to approximate the initial conditions as in Eq (2.1) this approach generates negative values. Thus, this modeling approach is not suitable to approximate the initial conditions $S(x, 0)$, and $I(x, 0)$.

2.1.3.   Piecewise cubic Hermite interpolating polynomial for continuous initial conditions

Another modeling approach for the initial condition is to approximate the uniform distribution over each region in Chile by a piecewise cubic Hermite interpolating polynomial (pchip) that is continuous and differentiable. In this way, the uniqueness of the solution of the PDE system (2) is guaranteed ^[38]. We can implement this approach by relying on the Matlab built-in function pchip that generates a polynomial $P(x)$. This polynomial has fewer oscillations than the spline interpolating polynomial when the data is not smooth ^[78]. This is relevant for the mathematical modeling approach since the polynomial $P(x)$ is used to approximate the initial population density of susceptible $S(x, 0)$ and infected $I(x, 0)$. Both populations should be positive for all times, and oscillations are not expected due to the nature of these populations. The pchip() function is less expensive to set up than the spline, so it has some additional advantages ^[78,79]. However, from a numerical computational point of view, there is still a problem since, despite the continuity and differentiability of the initial conditions, we have a large gradient in some of the boundaries between the regions. This generates a numerical challenge for the numerical solution and the calibration process, since piecewise cubic Hermite polynomials with a few interpolating points fail to accurately approximate the different initial conditions for the susceptible $S(x, t)$ and infected $I(x, t)$ subpopulations. However, if more interpolation points are used to gain accuracy, then the piecewise cubic Hermite polynomials increase their oscillations. This creates an issue for the numerical solution, and the initial conditions also become unrealistic due to the oscillations. In addition, unrealistic internal migrations within a region occur due to the oscillations.

The left-hand side of Figure 4 illustrates an example of a piecewise cubic Hermite interpolating polynomial for the initial condition $I(x, 0)$. Note the large gradients between regions and particularly between the Metropolitan and O'Higgins regions (see Table 2). This creates a computational issue since the pdepe() solver requires a finer grid. Moreover, the values of the integrals of the function $N(x)$ that represent the population in each region become inaccurate as time increases due to the large gradients and diffusion between regions. On the right-hand side of Figure 4 we show a piecewise cubic Hermite interpolating polynomial that approximates the initial condition of the regions 7 and 8, where this aspect can be seen.

2.1.4.   Gaussian initial condition for Chile

Another potential modeling approach for the initial condition is to use a Gaussian distribution function to distribute the population throughout Chile. In this way, the population density would be higher in the middle and lower at the northern and southern boundaries of Chile. Moreover, the function that approximates the initial condition is continuous and differentiable. This approach seems plausible, particularly for the southern region of Chile since the population density is low there. However, this particular approach can be misleading since each region of Chile has its own population density. In this approach, the initial condition can be written as

where $\sigma$ is a parameter that allows to control the dispersion or variance of the population over the space and $s$ is a scale factor that can be chosen in order to obtain a specific population when integrating. In this case we have normalized the length of Chile to $L = 16$, then the midpoint of Chile is represented at $x_0 = 8$. This approach has one drawback since the density is symmetric, which might not be a realistic approximation but is one of the challenges that real-world modeling brings and that we are pointing out in this research. Note that this initial condition satisfies approximately the compatibility condition with the boundary condition since at the boundaries the derivatives should be approximately zero.

2.1.5.   Sum of Gaussian functions for the initial conditions

We can use a more realistic and accurate modeling approach for the initial condition by using a Gaussian function in each region, where the integral of the Gaussian function gives the initial subpopulations of each region ^[80,81]. However, this has a main drawback since this approach generates several discontinuities in the boundaries between the regions in Chile, and then numerical issues arise. A better modeling approach is to approximate the initial condition of the subpopulations by a sum of Gaussian functions where each function is centered at the midpoint or center of each region ^[82,83,84]. This modeling approach solves the discontinuity aspect. However, careful attention is needed since all the Gaussian functions overlap over the whole space that represents Chile. This computational issue can be partially addressed by selecting a variance such that the overlap between Gaussian functions is very small ^[85]. The smaller the variance, the less overlap occurs. Nevertheless, a very small variance might create numerical issues due to the high mobility within each region. Thus, there is a trade-off between the overlapping of the Gaussian functions and the mobility of the subpopulations. Therefore, in this study, one of the modeling approaches that we utilize for the initial condition is the following

where $n$ is the number of regions in Chile, $x_i$ are the midpoints of each region $i$, and $N_i$ are the initial population of each region $i$. The main advantage of this approach is that the function $N(x)$ is continuous and differentiable. This allows us to guarantee the existence and uniqueness of the solution of the PDE system (1.1). In addition, it allows to have larger population densities at the middle of the regions in Chile and therefore to have more specific locations for the subpopulations. Variations on this modeling approach can be implemented and especially scale factors can be used. Notice that this particular initial condition satisfies approximately the compatibility condition with the boundary condition since at the boundaries the derivatives should be approximately zero. Also, note that using a sum of Gaussian functions is similar to using a mollifier for the initial conditions. Indeed, there are many choices for the mollifier and therefore for the ICs ^[86,87].

The graph on the left hand side of Figure 5 shows a function composed by a sum of Gaussian functions that approximates the initial condition of the susceptible population density $S(x, 0)$ of Chile. As it can be seen there is a relatively small overlap between the Gaussian functions. Using a variance $\sigma = 1/6$ for each Gaussian function we can guarantee that 97.3% of the population in each region is included in their respective region. The graph on the right hand side of Figure 5 shows the initial condition of the infected population density $I(x, 0)$ of Chile. It can be seen that the values are much smaller since at the early phase of the COVID-19 pandemic there were few infected people. For the initial recovered people we assume that $R(x, 0) = 0$.

2.2.   Boundary conditions

The boundary condition is given by Eq (1.3). This condition implies that there was not traveling of people from Chile to the northern neighbors countries, Bolivia and Peru, due to the lockdown restrictions during the early phase of the COVID-19 pandemic ^[32,88,89]. This is just an approximation of the real-world scenario where the traveling between the borders was relatively small. On the other hand, the most southern region of Chile is delimited by the South Pacific Ocean, so the people cannot move southern of that region. Thus, the selected boundary condition (1.3) is plausible. Other boundary conditions could be taken in order to consider a small flux of people during the early phase of the COVID-19 pandemic ^[88,89,90].

2.3.   Transmission rates

In this section we present some options regarding the transmission rates and the form of the force of infection ^[12]. There are many different forms that can be used for the transmission rates and the form of the force of infection of the mathematical model (1.1). The most common forms for the force of infections are the mass and standard bilinear incidence ^[12]. We can set the transmission rate using the following forms:

● Constant: In this case there is only one transmission rate $\beta$. This form seems unrealistic since it could be possible that the transmission rate may be different in each region due to social behavior or population density ^[30,89,90].

● Space varying: In this case transmission rate $\beta(x)$ is varying over the space $x$. This modeling approach has many options. For instance, it can be assumed that it varies depending on the total population density $N(x)$. Thus, assuming standard incidence one gets that $\beta(x) = \beta/N(x)$. Note that the units of $beta$ here are different from the one of the previous case. However, this option still requires to choose the initial condition $N(x, 0)$ for the distribution of the population over the space $x$. Note, that the distribution of the population $N(x)$ is affected by the diffusion of the mathematical model (1.1). This form can be suitable when the transmission rates are relatively different over the regions.

● Space and time varying transmission rates $\beta(x, t)$. This form implies that the transmission rate varies over regions and also over time. This form includes the possibility that the NPIs or social behavior affect the transmission rate over the early phase of the COVID-19 pandemic. This form is more complex and flexible and might be more realistic.

2.4.   Model's calibration process

As we have mentioned above, we consider Neumann boundary conditions and specifically with no flux due to the assumption of COVID-19 pandemic border traveling restrictions. This is a modeling approximation of the real-world scenario where the mobility between the borders was relatively small ^[89,90]. As mentioned previously, a particular attention is required to setup the initial condition for each of the populations of susceptible, infected and recovered. It is reasonable to use the total population of each region to set the initial condition. We have mentioned some options to achieve this. In the next section we will present some numerical results by using different initial conditions.

Let us focus on the calibration process of the spatial-temporal mathematical model (1.1). The first step of the calibration process consists in defining an objective function. For spatial-temporal models, the calibration process requires careful attention since the state variables are defined over space and time ^[38,91]. In the spatial-temporal model (1.1) the state variables $S(x, t)$, $I(x, t)$, and $R(x, t)$, represent population densities over the one-dimensional space variable $x$. Therefore, for a particular value of $t$ and over a specific region $\Delta x$, we have that the integral $\int I(x, t) dx$ represents the number of infected people over that region and at some particular time $t$. The calibration process can be done in many ways depending on the objective function and the number of parameters that need to be estimated ^[92]. Regarding the objective function, we can use the classical least square error between the infected data and the simulated infected data. With regard to the parameters, the calibration process and identifiability depends on the assumptions regarding the parameters. The most simple case is when the transmission rate is constant and the diffusion is fixed beforehand. In this particular case, the calibration process can be done by finding the optimal $\beta$ that minimizes the SSR (sum of squared residuals). Thus, the objective function is given by

where $\hat{I}_{i, j}$ is the infected cases corresponding for the $i^{th}$ region and $j^{th}$ week; $I(x, t, \beta)$ is the infected population density, which depends on the value of $\beta$; and $\Delta x_i$ corresponds to the $i^{th}$ one-dimensional region. Notice, that in this case, the only parameter that needs to be estimated is the transmission rate since the parameter $\gamma$ can be fixed based on biological data. This approach allows the identifiability of the parameter $\beta$, but the mathematical model would have almost no freedom to adjust to the real data $\hat{I}_{i, j}$. More complex calibration processes can be implemented by for instance leaving the diffusion coefficients as parameters to be estimated, but the identifiability of the model would not be guaranteed ^[92]. Moreover, assuming that the transmission rate and the diffusion depend on the space variable $x$ and time $t$ would make the calibration process unmanageable due to the great freedom of the functions $\beta(x, t)$ and $d_i(x, t)$. Another addtional aspect that can be modified when solving the inverse problem is to choose a different objective function. For instance, using logarithmic values of the data and the model or choosing different weights for each of the terms in Eq (2.5) ^[93,94].

3.   Numerical results for different mathematical modeling approaches

In this section we present some numerical results when different mathematical modeling approaches are implemented. A variety of approaches are applied in order to get insight into the challenge of the mathematical modeling of the early phase of the COVID-19 pandemic in Chile. From a practical point of view the modeling aspects are quite different from the theoretical ones where oftentimes we are interested in qualitative results such as existence and uniqueness ^[38]. In this study we focus on the early phase of the COVID-19 pandemic in Chile, which is the transient dynamics.

In this study, the mathematical modeling approaches vary depending on the following:

● The form used to approximate the initial conditions

● The transmission rate and force of infection

● The level of diffusion or mobility of people

Other approaches that can be implemented, such as different boundary conditions, but in this study, we focus on the previous ones which provide enough versatility to the mathematical modeling approach. In addition, we initially implemented the models with regions measuring 1 unit size and Chile 16 units. However, we also implemented other models with different lengths in order to approximate more accurately the real-world scenario. In order to numerically solve the spatial-temporal mathematical model (1.1) we rely on the pdepe() built-in Matlab function. Although, the solution of the mathematical model (1.1) depends on $(x, t)$, and the graph is a surface, in this section we present the numerical results using curves that represent the subpopulations in each region of Chile. Thus, we are able to compare with the epidemic data that is available by regions and not given by densities. An alternative approach would be to approximate the densities in each region and compare with the numerical solution (compare two surfaces). However, this will create additional challenges to approximate the subpopulations densities.

Mathematical model (1.1) has five parameters related to the COVID-19 pandemic. The average infectious period $\gamma$ is a biological parameter that can be assumed fixed in the model ^[42,44]. However, the diffusion coefficients $d_{i}$ and transmission rate $\beta$ are unknown. Approximating or estimating these unknowns is very challenging due to the lack of specific data and also due to identifiability aspects ^{[4,92,95,96,97,98,99]}. Thus, in this study we focus on presenting the mathematical modeling approaches, challenges and the calibration of the model to the data. We perform some soft calibration processes in order to illustrate the applicability and challenges of the model (1.1) to real-world scenarios. In this way, we can provide insight into the construction of mathematical modeling approaches.

3.1.   Modeling approach with constant SARS-CoV-2 transmission rate $\beta$

For the first mathematical modeling approach that consists in the use of the spatial-temporal model with a space-invariant transmission rate we approximate the initial condition by using the piecewise cubic Hermite interpolating polynomials. We can try to calibrate the model and in particular the transmission rate $\beta$ to resemble the spatial-temporal data of Chile. The first model's calibration has difficulties approximating the dynamics of the early COVID-19 pandemic in Chile due to the constant transmission rates over the regions of Chile.

Figure 6 shows the dynamics of the early phase of the spread of SARS-CoV-2 in Chile by using a constant transmission rate $\beta$ and low diffusion rates. It can be seen that the model fits relatively well to the infected cases for the growing phase in the Metropolitan region. Nevertheless, in the other regions the model is not able to approximate well to the infected cases. Increasing the transmission rate $\beta$ can improve the fit for the regions different than the region 7, but the the fit of the model to the 7 region would be inferior since there would be a large increase of infected cases in the region 7. We can infer that performing a calibration process for a constant transmission rate would try to adjust the model to the data from the Metropolitan region since the number of infected cases is larger and therefore it has more impact on the SSR. We conducted additional tests using the initial condition in the form of a sum of Gaussian functions as in Eq (2.4). While this initial condition addresses the numerical issues of the Hermite polynomials at the boundary zones between regions, the modeling approach fails to achieve a good calibration. Thus, we can suggest that the modeling approach with a constant transmission rate is not suitable to characterize the early phase of the COVID-19 pandemic in Chile. Later, we will introduce a time varying transmission rate that reduces the transmission rate due to NPIs in order to resemble the decay phase of the infected cases.

3.2.   Modeling approach with a space-variable transmission rate $\beta(x)$

In this section we use a mathematical modeling approach with a space-variable SARS-CoV-2 transmission rate $\beta(x)$. There are many options for accomplish this approach. Notice that the transmission rate depends on the space variable $x$. Thus, the transmission rate in the PDE system (1.1) can be replaced by a function $\beta(x)$ that varies in terms of the space variable. Therefore, one gets the following PDE system:

The initial and boundary conditions can be written as in Eqs (1.2) and (1.3). Note that we have many options to choose for the functions that represent the initial conditions and this can affect the dynamics of the spread of SARS-CoV-2 since we are studying the transient dynamics (early phase of COVID-19 pandemic). Therefore, careful attention is required to select the mathematical forms of the initial conditions.

3.2.1.   Space-variable transmission rate $\beta(x)$

Let us assume that the transmission rate depends on the population or population density of each region ^{[100,101,102,103]}. Thus, we assume an approximated standard bilinear incidence where the transmission rate in the PDE system (3.1) is replaced by

where $N$ is the total population density as presented in Eq (2.4). This form gives rise to a standard incidence, which implicitly assumes that an individual has a fixed number of contacts per unit time ^[12]. The factor $N_i$ scales the Gaussian functions of each region and therefore for the regions with larger populations the transmission rates are decreased. As a consequence, the transmission rate would increase such that the models fits well the region 7 and therefore opens the possibility of having higher peaks in the other regions.

Before we proceed with the numerical simulations, first we need to distribute the initial population by using the function $N(x)$ given by a sum of Gaussian functions (see Eq [2.4]). Since the numerical simulation provides the density of the populations we need to integrate over the space of each region in order to obtain the total population in each region at a specific time.

Figure 7 shows the dynamics of the infected cases that comes from the numerical simulation of mathematical model (3.1) and the real data from Chile. It can be seen that this mathematical modeling approach is able to improve the approximation of the real dynamics by means of the space-variable transmission rate $\beta(x)$ across the regions of Chile. Notice that the epidemic wave occurs in all the regions, but overestimates the real data. We have used a scale factor in order to adjust the peak of the model in the region 7 to the real data. This scale takes into account the fact that there are asymptomatic cases which are not reported and therefore are not included in the real data ^[104,105]. It also shows that in the regions located beside the high dense region 7 there are large epidemic waves due to smaller population densities at the neighborhoods of the boundaries between regions. In Figure 8, it can be seen that in the last step of the simulation there is a large number of infected cases in those neighborhoods. For instance, at the first time steps we can see that the infected cases are distributed as Gaussian functions with mean values located at the middle of the regions and then during the numerical simulation the infected cases shift to the boundaries where the transmission rate is higher. In the last time step of the simulation we can see that the infected cases are located at the boundaries between regions.

3.2.2.   Space-variable transmission rate $\beta(x)$ in terms of sum of scaled Gaussian functions

In this section we assume that the transmission rate function $\beta(x)$ is given by a sum of scaled Gaussian functions. In this way, we can make the transmission rate larger in the middle of the regions and smaller in the boundaries. From a real-world viewpoint this translates in that there are more contacts in the cities and that the cities are closer to the middle of the region. The transmission rate function $\beta(x)$ is given by

where $n$ is the number of regions in Chile, $x_i$ are the midpoints of each region $i$, $\sigma_i$ is the density variance of each region $i$, and $N_i$ is the initial population of each region $i$. Furthermore, we define $\sigma_i = \frac{L_i}{6 L_C }$ where $L_i$ is the length in spacial units of region $i$ and $L_C$ is the length in the same spacial units of the entire country of Chile. In the first numerical models we set all regions to measure 1 unit and Chile 16 units. After Subsection 3.3, we will use the lengths in km indicated in the third column of Table 2. Defining $\sigma$ in this manner guarantees that $99.7\%$ of the region's density is located within its spacial limits. This transmission rate function $\beta(x)$ is similar to the one presented in Eq (3.2), but uses a particular scale factor $s_i$ for each region and guarantees that the regional centers of mass are located in each regions midpoints $x_i$. Thus, we can avoid very large transmission rates at the boundaries between regions. Note that the parameter $\beta$ modulates the transmission rate over the whole space and therefore interconnects in some way all the regions. Moreover, it allows to model lower mobility of people at the boundary of the regions. Figure 9 shows the transmission rate function $\beta(x)$. Note that the transmission rate is larger in the middle of the regions and smaller in the boundaries, which is a crucial aspect of this particular modeling approach.

3.2.3.   Results with high diffusion

In this subsection, we show the impact of diffusion on the transient dynamics of the COVID-19 pandemic. We set the diffusion coefficients as time and space invariant. Figure 10 shows the dynamics of infected cases resulting from the numerical simulation of mathematical model (3.1) with high diffusion. It can be seen that the infected population is larger at the last step of the simulation and the population density is larger in the middle of the regions. Figure 11 partially shows the transient dynamics of the total population density $N(x)$ resulting from the numerical simulation of mathematical model (3.1) with a high diffusion and using the transmission rate defined in Eq (3.3). It can be seen that the initial population profile is partially lost at the end of the simulation. This is due to the high diffusion and therefore the population density becomes more homogeneous regionally thereby reducing the distinctiveness of the boundaries. This result is expected and can be interpreted as people moving from regions of high population density to ones with low density. Notice that with accurate data on mobility, the diffusion coefficients could be estimated. However, due to quarantine interventions and the complexity of data gathering, most likely there are no records of mobility between regions during the early phase of the COVID-19 pandemic in Chile.

Figure 12 shows the dynamics of infected cases derived from the numerical simulation of mathematical model (3.1) compared to the real data from Chile. It can be seen that this mathematical modeling approach improves the previous results by scaling the transmission rate, but still overestimates the real infected data in all regions except the Metropolitan region.

3.2.4.   Results with low diffusion

Here we show the transient dynamics for a low diffusion. Figure 13 shows the dynamics of infected cases and susceptible when diffusion is low. It can be seen that the density of the infected population initially grows. A very interesting behavior is observed here: in each region, the middle has a lower infected population density due to the fact that the transmission rate is higher in the centroids of each region. This higher transmission causes the density of the susceptible population $S(x, t)$ to decrease in the centroids due to the infection process. Therefore, due to the reduced number of susceptible people individuals in the central areas of the regions, a delayed decline in the number of infected individuals these central zones occurs. This epidemiological phenomenon is related to the herd immunity effect ^[106,107]. Figure 14 shows the dynamics of infected cases in Chile. It can be seen that the model overestimates in all regions, since now people reduce their mobility from the spots with high transmission rate. A different modeling approach is needed to better resemble the real dynamics.

3.3.   Modeling with regions' length and with a spatial-temporal transmission rate $\beta(x, t)$

In this section, we consider spatial-temporal mathematical model (3.1) with the approximated length of the regions, since this might improve the approximation of the real dynamics. In this way, we are taking into account the space of each region. This affects the population density $N(x)$ and therefore the other state densities $S(x, t)$, $I(x, t)$, and $R(x, t)$. This may have some effect on the dynamics of the COVID-19 pandemic and could differ from the dynamics obtained using other previous modeling approaches. This particular modeling approach seems more realistic than the previous ones, and more accurate results may be expected despite the great freedom of the mathematical model (3.1). In addition, we introduce a time varying transmission rate $\beta(x, t)$ of the following form:

Figure 15 shows the numerical solution of the spatial-temporal mathematical model (3.1) with a time-varying transmission rate $\beta(x, t)$ in order to include a NPI (lockdown) around week 12.2. In particular, for $t > 12.2$ we reduce the function $\beta(t)$ by half. The numerical values of $\beta$ before and after the time peak were found by doing a manual calibration to avoid a very large computational time since the main aim is not to solve exactly an inverse problem. The model is able to relatively approximate the dynamics of the infected cases in the region 7 of Chile. However, improvements are needed for the other regions. Note that the regions closer to the Metropolitan region are smaller. We assumed low diffusion rates and in particular reduced by $10^{-2}$ for the infected people due to their symptoms and lower mobility. In addition, we considered low mobility during the early phase of the COVID-19 pandemic due to the NPIs in Chile. Interestingly, in ^[108], the diffusion term was only include in the infected people and disregarded for the susceptibles.

3.4.   Model's calibration for each region with a spatial-temporal transmission rate $\beta(x, t)$

We have shown that calibrating the mathematical model (3.1) to real infected data provides better results than the spatial-temporal model (1.1). Nevertheless, the calibration process is still complex due to the variety of regions and dynamics of each region. One alternative modeling approach is to calibrate model (3.1) for each region (independently) in order to obtain particular spatial-variable transmission rates in each region. In this approach, there is no mobility between regions, since the model is used for each region individually. In this section, we provide some results regarding this alternative modeling approach.

We calibrated model (3.1) for each region of Chile and used the Neumann boundary conditions in each region. Note that here, the calibration of each region is performed independently. The numerical result for each region improves the previous results. In Appendix, Table 3 we present a summary of the parameters values obtained from the calibration for each region. Figure 16 depicts the calibration of model (3.1) by region. The model adjusts relatively well to the infected population. However, it is important to mention that to calibrate the model, we considered a spatial-temporal transmission rate $\beta(x, t)$ as the one presented in Eq (3.4). We introduced only one change per region for the parameter $\beta$ at a time $t_{pi}$ in order to mimic the effect of NPIs on the transmission rate ^[30,33]. This approach is plausible since in each region the social behavior might be different and the NPIs were implemented at different times ^[32,37]. The data for the southern regions show epidemiological behaviors without any recognizable pattern. Only the the region of Los Lagos (region 14) shows a typical epidemiological behavior. The irregular behavior of the COVID-19 pandemic in the southern regions could be due to spatial effects or social behavior related to NPIs ^{[109,110,111]}. Also, the most southern regions have a much lower population density, and thus the homogeneous mixing underlying assumption might not be a good approximation of the reality. Some of the southern regions show multiple waves that might be related to spatial effects or data collection ^{[112,113,114]}.

3.5.   Global model's calibration for Chile with a spatial-temporal transmission rate $\beta(x, t)$

In this section we improve all the previous model's calibrations by introducing a more realistic transmission rate. In particular, we utilize a spatial-temporal transmission rate $\beta(x, t)$ and a global calibration is performed by simultaneously including all the regions of Chile. Therefore, there is mobility of people between regions. In this section, we scale the global transmission rate $\beta(x, t)$ in order to calibrate the model (3.1) to the data. The spatial-temporal transmission rate is given by a sum of scaled Gaussian functions with time-varying scales $s_i(t)$ as follows:

The scale factor where $s_i(t)$ is a piecewise function that allows to modulate the transmission rate of each region during the growth and decay phases ^[6,115]. This gives the model the flexibility to replicate the impact of government imposed restrictions on mobility after the peaks of infected individuals are attained.

Figure 17 shows the infected population resulting from a numerical simulation of spatial-temporal model (3.1). The model adjusts relatively well to the data despite the high irregularity of the data. This results show the potential use of the spatial-temporal model (3.1). For some regions, the real data show irregularity and several waves that might be due to a combination of spatial factors and the implementation of NPIs. In addition, as we have mentioned before, the social behavior of the population could cause these temporal waves ^{[19,20,113,116]}. The regions that are more challenging to approximate by the mathematical model are regions 9-16. All these regions except region 14 feature a second epidemic wave. These regions are from the southern part of Chile, and it might be possible that spatial factors play a role in the dynamics of the COVID-19 pandemic ^{[19,20,113,116]}.

Figure 18 shows a normalization of the infected population resulting from a numerical simulation of spatial-temporal model (3.1). The computation time is approximately 11 minutes in a computer with 12th Gen Intel(R) Core(TM) i7-12700 2.10 GHz and 64-bit Windows 10 Enterprise LTSC operating system. The plot shows all the numerical results for each region by using an ad-hoc normalization that enables to observe the results for all the regions.

4.   Summary of the spatial-temporal mathematical modeling

● In the scenario of a constant transmission rate $\beta$, we obtained that the region 7 has the highest increase in cases due to a large population density. However, in the other regions, the epidemic waves generated by the mathematical model underestimated the real data.

● For the simulation results we took into account that some proportion of the infected people is not reported.

● In the scenario of a space-variable transmission rate $\beta(x)$, the transmission rate is reduced in the most populated regions to obtain higher epidemic waves to match real data. The total population density plays a crucial role. We analyzed a few options that can be used for $N(x)$.

● For the standard incidence, choosing a sum of standard Gaussian forms for the total population density $N(x)$ creates a computational issue since it generates large values for the transmission rate $\beta$ at the boundaries between the regions.

● For the standard incidence, the scaled Gaussian forms for the total population density $N(x)$ avoids the computational issue at the boundaries of the regions, but still the modeling approach is not still able to resemble the real dynamics of all the regions.

● Taking into account the length of the regions creates a more realistic scenario. Nevertheless, in this scenario, the model approximates the region 7 but underestimates the real data for the other regions.

● We also independently calibrated the spatial-temporal model to each region. In this case, the calibration improves, but at the expense of the lack of spatial connection between regions. This local calibration helped with the global calibration of the model for Chile.

● In the scenario of a spatial-temporal variable transmission rate $\beta(x, t)$, where we varied the scale piecewise functions, we obtained relatively good calibrations due to its versatility.

5.   Discussion

We presented and analyzed several mathematical modeling approaches using variants of a spatial-temporal epidemiological model. We performed numerical simulations that attempted to resemble the spatial-temporal data of Chile and provide insight into the suitability of the constructed mathematical modeling approaches. One first approach considered a space-invariant transmission rate, where it is implicitly assumed that the SARS-CoV-2 transmission rate does not vary between regions. This later modeling approach is very common in many mathematical models that rely on ODEs. We found that the spatial-temporal model with constant transmission rate was not suitable to resemble the dynamics of the early phase of the COVID-19 pandemic in Chile. In particular, in this approach, only the region 7 is resembled approximately when a least-square calibration process is performed. This is due to the fact that the large population of the metropolitan region allows to use a relatively small transmission rate that is not large enough to have epidemic waves that resemble the dynamics of other regions.

Another crucial aspect of the mathematical modeling approach is the mathematical form used for the initial conditions. The best option was to assume that the initial population density is concentrated in the middle of the regions. Other options such as using a uniform distribution can generate a large diffusion or mobility at the boundaries between the regions. Moreover, a uniform distribution creates a computational challenge since an approximating function for the initial conditions is necessary. We explored the use of the classical splines and the piecewise Hermite interpolating polynomials. We found that the latter one was more suitable for when a uniform distribution of the population is used.

In the second general mathematical modeling approach, we used different space-variant transmission rates. First, we utilized a space-variable rate form $\beta(x)$, which is composed of the sum of sixteen Gaussian functions in order to have higher transmission rates in the middle of the regions. This mathematical form also includes a factor related to the number of people in each region in order to reduce the transmission rate in the more dense regions. This form is similar to the well-known standard incidence but uses spatial effects ^[12]. With this approach, we were able to improve the previous results of a constant transmission rate. Nevertheless, the model was not able to resemble the real dynamics since using a transmission rate that only depends on the population density of the region is very restrictive. Therefore, we designed a new transmission rate form $\beta(x)$ that has a scale factor $s_i$ in each of the sixteen Gaussian functions. Thus, the scale factor of each region can modulate the spread of the SARS-CoV-2. Using this modeling approach the mathematical model was able to approximately resemble the growth of the COVID-19 pandemic, but not very well the decay phase. Therefore, we designed a new modeling approach for the transmission rate.

In the third general modeling approach we considered a spatial-temporal varying transmission rate. This approach has more flexibility since it allows to include NPIs in each of the regions. We implemented only one NPI per region, but further NPIs can be included. Furthermore, we assembled the NPIs at different times in order to approximate the real dynamics for each region. In addition, in this approach we took into account the length of each region of Chile for the initial conditions, the transmission rate, and the equations of the model. This provides a more realistic scenario for the population density, which can affect the real dynamics. All these nice features come with a drawback from a mathematical viewpoint. Adding more freedom to the mathematical model implies that more parameters need to be estimated. However, there is not enough available real data to allow us to uniquely identify the values of some of the parameters, such as the distribution of the population over the space ^{[16,17,55,117]}.

It is important to mention that the mathematical modeling approach presented in this study faced some important challenges due to the lack of data and also due to the heterogeneous effects of NPIs. For instance, in ^[32] it was mentioned that some regions (e.g., regions 11–13) presented completely different trends regarding COVID-19 cases. Another aspect that complicates the mathematical modeling of the spatial-temporal dynamics is the fact that some NPIs were imposed by communes (administered by municipalities) and not by regions. Moreover, individual communes in Chile implemented different quarantine protocols ^[29,34,118]. In addition, the numerical results presented in this study suggest that the spread of SARS-CoV-2 in the early phase of the COVID-19 pandemic in Chile was heterogeneous by region. These results agree with the statistical analysis presented in ^[32]. Our results highlighted the fact that the spatial-temporal mathematical modeling approach needs careful attention regarding calibration since there were heterogeneous NPIs during the early phase of the COVID-19 pandemic in Chile. One important feature was the introduction of quarantine interventions especially during June and July of 2020 in the regions 6–8. These regions have a great percentage of cases concerning to the total cases in Chile. We observed that the mathematical modeling approach that does not consider explicit temporal changes in the transmission rate has difficulties resembling the real data during the downward trend of the June and July months of 2020. We introduced time-varying transmission rates in order to deal with the introduction of quarantine interventions around the peaks of the first wave in each region ^[36].

One important limitation of many continuous spatial-temporal models such as the ones presented in this work, resides in the imposition of uniform population diffusion. This approach fails to accurately represent the relocation of family units that transition from densely populated areas to sparsely populated ones through discrete and isolated means of transportation, such as private vehicles. Such dynamics could be captured by agent-based models which include spatial and temporal effects ^[119]. Nevertheless, we think that the spatial-temporal models provide plausible insights into the dynamics of COVID-19 pandemic, especially for regions where the population is approximately distributed over all the region and therefore the density can be considered continuous. In addition, during the COVID-19 pandemic, people living in large cities left those cities and traveled to regions with a lower population density ^[120]. Thus, the reaction-diffusion models used in this work are suitable.

6.   Conclusions

In this paper, we constructed several spatial-temporal mathematical models to characterize the dynamics of the early phase of the COVID-19 pandemic in Chile. The model divides the total population into susceptible, infected, and recovered individuals. A system of first-order PDEs is used to build the mathematical models in order to include spatial and temporal effects. Thus, the model was designed to include the diffusion effects of SARS-CoV-2 in Chile and also to take into account different social behaviors and NPIs in the regions. First, we implemented and studied the spatial-temporal model with a constant SARS-CoV-2 transmission rate. Secondly, we implemented and studied the spatial-temporal model with different modeling approaches using a variety of space-variant SARS-CoV-2 transmission rates. In a third modeling approach, we constructed the model with a space-time-variant SARS-CoV-2 transmission rate and considered the length of each region of Chile in order to make the model more realistic.

We found that the first modeling approach has difficulties approximating the dynamics of the early COVID-19 pandemic in Chile. The second modeling approach was able to improve the previous result, but was still not able to approximately resemble the real dynamics. The third modeling approach with a space-time-variant SARS-CoV-2 transmission rate was able to better approximate the real dynamics due to its flexibility. In summary, we have shown that calibrating the mathematical model (3.1) to the real infected data provides better results than the counterpart model (1.1). Nevertheless, the calibration process is still very complex due to the variety of regions and dynamics of each region.

The results presented in this work show the advantages and drawbacks of the mathematical modeling approaches proposed to describe the early phase of the COVID-19 pandemic in Chile. The whole mathematical modeling process is explained in detail from the first step in order to show the challenges of the modeling process. Thus, this work can help other studies that deal with spatial-temporal epidemics. In addition, the results suggest that the SARS-CoV-2 transmissions in the Chilean regions were different.

We think that the construction of the models presented in this work and the challenges related to modeling the COVID-19 pandemic with spatial-temporal real data with a mathematical model based on PDEs are significant contributions. These contributions can help with future pandemics and also provide a deeper understanding of what happened during the COVID-19 pandemic. The results provide additional insights into the study of COVID-19 pandemics, since very few studies have explored similar approaches and even less with real-world data ^[56,108]. We show that the calibration of the mathematical model under different scenarios is very challenging due to the features of the model and the parameters that affect the dynamics. Calibrating or estimating the parameters of a spatial-temporal model, such as the ones presented in this work, is a difficult task, and future work can explore this avenue to improve knowledge in the field ^[121]. Thus, this paper also provides a didactic component for mathematical modeling with systems of PDEs. Future works can also study the mathematical model by considering space-time-varying diffusion rates. This approach is even more complex, since there are more parameters in the model.

Future studies can incorporate additional features of the spatial-temporal mathematical model in order to get more information about the characterization and modeling of the early phase of the COVID-19 pandemic in Chile and other countries. Extending the spatial-temporal model to other countries requires adding one more independent spatial variable since in reality SARS-CoV-2 spreads in a two-dimensional spatial coordinate system. However, the main ideas proposed in this work can be extrapolated when considering an additional spatial variable. For instance, in order to do the calibration of a two-dimensional spatial model, performing double integrals over the regions is required. Also, the transmission rate and diffusion coefficients can vary depending on the two spatial variables. This increases the variability of the spatial-temporal model, but the underlying mathematical approach is similar to the one presented in this work.

As in any mathematical modeling study, there are limitations that we need to be aware of. The main limitation of this study is the non-identifiability of the parameters of the spatial-temporal model due to the lack of some specific spatial-temporal data. It is important to point out this since this can encourage other researchers to gather data that can be used to estimate the parameters of the model. For instance, having the specific initial population density for the susceptible, infected, and recovered populations would allow for a more accurate model and in that way better characterize the dynamics of the early phase of the COVID-19 pandemic. In addition, having the exact dates on which the NPIs were implemented in each region would improve the estimation of the temporal change in transmission rate. In reality, the transmission rate changes continuously, so having only one change is a coarse approximation. Finally, the SARS-CoV-2 transmission rate can also vary depending on the social behaviors of people from different regions, but estimating these variations is a very challenging task and future work might attempt to address this key aspect.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

First author has been supported by grant María Zambrano (UPV, funding from the Spain Ministry of Universities funded by the European Union-Next Generation EU) and an Institutional Development Award (IDeA) from the National Institute of General Medical Sciences of the National Institutes of Health under grant number P20GM103451. Cristina-Luisovna Pérez has been supported by an FPU21/02343 grant (Formación de Profesorado Universitario 2021) from the Ministerio de Universidades, Spain. The authors are grateful to the reviewers for their careful reading of this manuscript and their useful comments and suggestions to improve the content of this paper.

Conflict of interest

Gilberto González-Parra is a guest editor for Mathematical Biosciences and Engineering and was not involved in the editorial review or the decision to publish this article. All authors declare that there are no competing interests.

Appendix

In this section, we provide the data from Chile, results regarding model's parameters, and some additional numerical aspects of the solutions.

Data of Chile

Table 2 shows the population, population density and approximated length (north to south) for each region of Chile ^[57,58,59]. The regions are identified by name and number. The data is used for the initial conditions and spatial dimension of the different mathematical models constructed in this work.

Model's parameters

In Table 3 it is shown the scales that were found by calibrating the mathematical model (3.1) to the infected data of each region in Chile.

Numerical aspects of the solution

The calibration requires utilizing the Matlab built-in function pdepe() requires very small spatial and time step sizes in order to obtain a reliable numerical solution. The accuracy of the numerical solution can be improved by decreasing the spatial and time step sizes. Special attention is required whenever the density of the infected population $I(x, t)$ becomes almost zero. This might generate numerical issues and negative solutions may arise. There are numerical schemes that can guarantee positivity solutions for biological systems ^[122,123].In this work we can obtain positive solutions by choosing a relatively small spatial and time step sizes for the Matlab built-in function pdepe(). Figure 19 shows two numerical solutions at different times and different spatial and time step sizes. As it can be seen the accuracy of the numerical solution improves as the spatial and time step sizes become smaller. The numerical issue arise whenever the density of the infected population $I(x, t)$ becomes close to zero. This numerical aspect is faced by many numerical schemes that cannot guarantee positivity solutions for biological systems where the positivity of the system is guaranteed ^{[122,123,124]}.