Estimating Cross-Border Mobility from the Difference in Peak Timing: A Case Study of Poland–Germany Border Regions

Abstract

Human mobility contributes to the fast spatiotemporal propagation of infectious diseases. During an outbreak, monitoring the infection on either side of an international border is crucial as cross-border migration increases the risk of disease importation. Due to the unavailability of cross-border mobility data, mainly during pandemics, it becomes difficult to propose reliable, model-based strategies. In this study, we propose a method for estimating commuting-type cross-border mobility flux between any pair of regions that share an international border from the observed difference in their infection peak timings. Assuming the underlying disease dynamics are governed by a Susceptible–Infected–Recovered (SIR) model, we employ stochastic simulations to obtain the maximum likelihood cross-border mobility estimate for any pair of regions. We then investigate how the estimate of cross-border mobility flux varies depending on the transmission rate. We further show that the uncertainty in the estimates decreases for higher transmission rates and larger observed differences in peak timing. Finally, as a case study, we apply the method to some selected regions along the Poland–Germany border that are directly connected through multiple modes of transportation and quantify the cross-border fluxes from the COVID-19 cases data from 20 February to 20 June 2021.

1. Introduction

In recent times, sporadic infectious disease outbreaks have created massive disruptions not only in public health functioning but also in socioeconomic developments. Incomplete knowledge of disease transmission mechanisms and the unavailability of effective medical treatment cause difficulties in controlling disease spread. Advances in modern transportation technology further accelerate the spatial spread of infection and complicate the management of disease spread in affected regions [1,2]. Mobility can occur on different spatial scales within a country. The complex interconnections among multiple scales limit the predictability of spatiotemporal patterns of disease spread [3,4,5]. Alternatively, intercountry migration often acts to introduce new pathogens in nonaffected countries, and within-country mobility can then quickly spread the disease throughout the country [6]. Intercountry migration occurs primarily in two ways: international air travel and daily commuting between countries sharing a border. The availability of digital technology makes it possible to keep track of the fluxes of global air-travel passengers. Moreover, during an outbreak, it helps to implement measures like airport screening, quarantining upon arrival, and contact tracing, which reduce the risk of disease import due to air travel. On the other hand, for daily cross-country mobility, due to a lack of convenient real-time database management [7,8], it is often difficult to estimate the mobility flow, and, in turn, to quantify the component of infection risk because of human mobility.

Mathematical models are important tools for assessing the risk of disease import from a neighboring country and investigating the consequences of different control strategies [9,10]. Reliable information on cross-border mobility flow not only hints at the effectiveness of intervention strategies like lockdowns or border-control measures but is also a key input parameter for multinational epidemiological models in evaluating the performance of spatially explicit intervention strategies [11,12]. In very few cases, mobility flux information can be obtained from technology companies, public transit data, census data, and survey data. However, these data are very expensive to obtain and often not publicly available due to privacy laws. Therefore, the estimation of cross-country mobility flux remains a challenging task in the context of infectious disease modeling.

In spatial epidemiological analyses, mobility flux is usually accounted for via the gravity model, the radiation model, or related approaches [13,14,15,16,17,18,19]. However, despite their wide use in the literature, these models have certain limitations as they rely on tuneable parameters which can vary with spatial location [20]. Similarly, several theoretical frameworks have been proposed to reconstruct the underlying connection topology from the observed dynamics of the focal system [21,22,23]. These methods have been mainly applied to those systems that exhibit oscillatory dynamics, such as coupled chaotic oscillators [24,25,26], EEG time series data [27], or neuronal systems [28]. Recently, [29] proposed a method to estimate the coupling strength between two regions from invasion time, which is the time taken to reach the focal region after an infection is seeded in a neighboring region. In practice, however, the invasion time will not likely be a reliable enough basis for extracting mobility flow. This is because disease monitoring systems typically face the greatest difficulties in the early phase of a disease outbreak due to incompleteness or delays in reporting. Consequently, all existing methods for estimating spatial connectivity are either data-hungry or depend heavily on precise information related to data reporting. A need therefore exists for a practical method that overcomes these limitations and can estimate cross-border mobility from imperfect monitoring data.

To this end, we propose an intuitive approach based on maximum-likelihood estimation to quantify the cross-border mobility flux for a pair of regions situated in two different countries. Our method assumes the underlying disease dynamics are governed by a simple Susceptible–Infected–Recovered (SIR) model and estimates the mobility flux from the observed difference in the timing of the peak of infection between two regions. We first theoretically investigate, under realistic epidemiological assumptions, how the observed difference in peak timing can effectively recover the underlying mobility between two regions and how the estimate of mobility varies depending on the disease transmission rate. Next, as a case study, we use our method to estimate the cross-country mobility flux from COVID-19 incidence data for the pair of regions that are connected directly through multiple modes of transportation and are located on either side of the Poland–Germany border.

2. Methods

2.1. Deterministic Model

We consider a standard deterministic two-patch SIR model incorporating short-term or commuting-type migration. Each spatial unit in our study, such as a district or even a state that shares a border, can be considered a patch. In each patch i, the total population is divided into three classes based on health status: susceptible ($S_i$), infected ($I_i$), and recovered ($R_i$). Since the patches are connected through migration, the disease dynamics in each patch have two components: (i) disease dynamics within the patch and (ii) disease dynamics between the patches (see [30,31,32,33,34] and references therein). The mathematical equations describing the mechanism of disease dynamics in this two-patch scenario are given as follows:

(1)$\begin{array}{l}\frac{d{S}_i}{dt}=-{\displaystyle \sum _{j=1}^2}{\beta }_j{m}_{ij}{S}_i\left(\frac{{\displaystyle \sum _{k=1}^2{m}_{kj}{I}_k}}{{\displaystyle \sum _{k=1}^2{m}_{kj}{N}_k}}\right),\\ \frac{d{I}_i}{dt}={\displaystyle \sum _{j=1}^2}{\beta }_j{m}_{ij}{S}_i\left(\frac{{\displaystyle \sum _{k=1}^2{m}_{kj}{I}_k}}{{\displaystyle \sum _{k=1}^2{m}_{kj}{N}_k}}\right)-{\gamma }_i{I}_i,\\ \frac{d{R}_i}{dt}={\gamma }_i{I}_i,i=1,2.\end{array}$

The rate at which the susceptible population in patch i gets infected in patch j is given by the product of three terms: the risk of infection in patch j (${\beta }_j$), the number of susceptibles from patch i who are currently in patch j ($m_{ij}{S}_i$), and the proportion of the population that is infected in patch j ($\frac{{\sum }_{k=1}^2{m}_{kj}{I}_k}{{\sum }_{k=1}^2{m}_{kj}{N}_k}$). Here, $N_j$ denotes the total population in patch j and ${\sum }_{k=1}^2{m}_{kj}{N}_k$ gives the total population currently present in patch j, called the effective population of patch j [30,31]. The model, (1), has two epidemiological parameters: transmission rate (${\beta }_i$) and infectious period ($\frac{1}{{\gamma }_i}$). The matrix $m_{ij}$ is called the residence-time matrix, where the element $m_{ij}$ can be interpreted as the fraction of people from patch i that visited patch j [30,31]. Since $m_{ij}$ denotes the fraction of mobility flux, we have the constraint ${\sum }_{j=1}^2{m}_{ij}=1$, for $i=1,2$.

For a fixed set of model parameters and initial conditions, we define $t_d$ as the difference between the peak timing of these two regions, which can be expressed as follows:

$t_d=|t_1^{max}-t_2^{max}|,$

where, $t_1^{max}$ and $t_2^{max}$ denote the time points at which the trajectories of the infected populations $I_1\left(t\right)$ and $I_2\left(t\right)$ reach their respective maxima.

2.2. Stochastic Model

Here, we are interested in estimating the cross-border mobility flux from the observed difference in peak timing, $t_d$. It should be noted that in a deterministic setup, for given model parameters and initial conditions, we always end up with exactly one $t_d$. Therefore, we introduce randomness into the model as a form of demographic stochasticity. Here, demographic stochasticity refers to the fluctuations in the disease propagation process arising from the random nature of disease transmission and recovery events at the individual level. We consider a standard event-driven approach called the tau-leaping method [35]. This method is a modified version of the Gillespie algorithm, where the interevent duration is kept fixed instead of following an exponential distribution, which increases computational efficiency [35,36].

Let us denote the number of transmission and recovery events at time t for each patch i (i = 1, 2) by ${\mathcal{E}}_{\mathrm{Trans}}^i\left(t\right)$ and ${\mathcal{E}}_{\mathrm{Rec}}^i\left(t\right)$, respectively. The transition probabilities for the events ${\mathcal{E}}_{\mathrm{Trans}}^i\left(t\right)$ and ${\mathcal{E}}_{\mathrm{Rec}}^i\left(t\right)$ that occur during the small but fixed time interval $\delta t$ can be written as

(2)$\begin{array}{l}{\displaystyle \mathbb{P}(\delta {\mathcal{E}}_{\mathrm{Trans}}^i=1|S_i,I_i,R_i)=\sum _{j=1}^2{\beta }_j{m}_{ij}{S}_i\left(\frac{{\displaystyle \sum _{k=1}^2{m}_{kj}{I}_k}}{{\displaystyle \sum _{k=1}^2{m}_{kj}{N}_k}}\right)\delta t+o\left(\delta t\right),}\\ \mathbb{P}(\delta {\mathcal{E}}_{\mathrm{Rec}}^i=1|I_i)={\gamma }_i{I}_i\delta t+o\left(\delta t\right),\end{array}$

where, $\delta {\mathcal{E}}_{\mathrm{Trans}}^i={\mathcal{E}}_{\mathrm{Trans}}^i(t+\delta t)-{\mathcal{E}}_{\mathrm{Trans}}^i\left(t\right),$ and $\delta {\mathcal{E}}_{\mathrm{Rec}}^i={\mathcal{E}}_{\mathrm{Rec}}^i(t+\delta t)-{\mathcal{E}}_{\mathrm{Rec}}^i\left(t\right)$ for $i=1,2$.

For small $\delta t$, the increments $\delta {\mathcal{E}}_{\mathrm{Trans}}^i$ and $\delta {\mathcal{E}}_{\mathrm{Rec}}^i$ are approximated as Poisson distributed random variables:

(3)$\begin{array}{l}{\displaystyle \delta {\mathcal{E}}_{\mathrm{Trans}}^i\approx Poisson(\sum _{j=1}^2{\beta }_j{m}_{ij}{S}_i\left(\frac{{\displaystyle \sum _{k=1}^2{m}_{kj}{I}_k}}{{\displaystyle \sum _{k=1}^2{m}_{kj}{N}_k}}\right)\delta t),}\\ \delta {\mathcal{E}}_{\mathrm{Rec}}^i\approx Poisson\left({\gamma }_i{I}_i\delta t\right).\end{array}$

The population sizes in different compartments are then updated as

(4)$\begin{array}{c}\hfill \begin{array}{c}{S}_i(t+\delta t)=S_i\left(t\right)-\delta {\mathcal{E}}_{\mathrm{Trans}}^i,\hfill \\ I_i(t+\delta t)=I_i\left(t\right)+\delta {\mathcal{E}}_{\mathrm{Trans}}^i-\delta {\mathcal{E}}_{\mathrm{Rec}}^i,\hfill \\ R_i(t+\delta t)=R_i\left(t\right)+\delta {\mathcal{E}}_{\mathrm{Rec}}^i.\hfill \end{array}\end{array}$

Given model parameters ${\beta }_i,{\gamma }_i,N_i,I_i\left(0\right)$, for a particular residence-time matrix $m_{ij}$, we generate multiple stochastic realizations. In each realization, we calculate the difference in peak timing ($t_d$) and consequently the distribution of it. We show an example in Figure 1A,B, where we consider the parameters as ${\beta }_1=0.30$, ${\beta }_2=0.20$, ${\gamma }_1={\gamma }_2=0.07$, $N_1$ = 10,000, $N_2$ = 20,000, $I_1\left(0\right)=10$, $I_2\left(0\right)=1$, m = 0.05. For this parameter combination, the corresponding basic reproduction number is 3.79, which aligns with previous estimates [37]. Next, we generate multiple stochastic realizations, and for each of the realizations, we calculate $t_d$ (see Figure 1A) and plot its distribution as a histogram (see Figure 1B).

To produce the likelihoods, we first consider a range of values for the mobility parameter and divide it into 100 linearly-spaced values. Then, for each value of the mobility parameter, we produce 10,000 stochastic simulations, and for each simulation, we calculate $t_d$. The full range of $t_d$ values is then divided into 100 bins. Then, we calculate the number of simulations for which the corresponding bin contains the observed $t_d$ and, consequently, this gives the likelihood of mobility parameter given $t_d$. Finally, we take the logarithm of these likelihoods (i.e., log-likelihoods (LL)) and calculate the maximum likelihood estimate (MLE) of the mobility parameter given $t_d$.

3. Results

We first investigate the relationship between the difference in peak timing ($t_d$) and cross-country mobility ($m_{ij}$). For the sake of simplicity, we assume the mobility in both directions is the same, i.e., $m_{12}$ = $m_{21}$ = m. We also take the recovery rates to be the same for both the regions, ${\gamma }_1={\gamma }_2=\gamma$. We fix $\gamma =\frac{1}{14}da{y}^{-1}$ throughout our study unless specified otherwise. For the transmission rate, we consider two cases: (i) ${\beta }_1$ and ${\beta }_2$ are equal, i.e., ${\beta }_1={\beta }_2=\beta$, and (ii) ${\beta }_1$ and ${\beta }_2$ are not equal. We consider the total population $N_1$ = 10,000 and $N_2$ = 20,000, which is fixed over time, and we assume initially there is no recovered individual in either region, i.e., $R_1\left(0\right)$ = 0, and $R_2$(0) = 0.

Following the method described in Section 2.2, for each value of the mobility parameter (m), we take 10,000 realizations and calculate $t_d$ for different combinations of the initial infected population ($I_1\left(0\right)$ and $I_2\left(0\right)$) and disease transmission rate ($\beta$). We see that the difference in peak timing ($t_d$) gradually decreases with increasing mobility flux m (see Figure 2A–F).

We see that for a particular value of m, the difference in peak timing is higher for low disease transmission rate, and if we increase the transmission rate it starts to decrease. For instance, if we fix initial infecteds as $I_1\left(0\right)=1$, and $I_2\left(0\right)=10$, we see that for lower transmission rate $\beta =0.15$, and mobility parameter $m=0.025$, the difference in peak timing varies in the range [0, 18] days (see Figure 3A). If we slightly increase the transmission rate, i.e., $\beta =0.2$, we observe that the range of $t_d$ shrinks to [0, 10] days (see Figure 3B). For a higher transmission rate ($\beta =0.3$), the range of $t_d$, in this case, reduces to [0, 6] days (see Figure 3C). Similar patterns are observed even if we change the number of initially infected individuals. For $I_1\left(0\right)=10$, and $I_2\left(0\right)=1$, we see, for $\beta =0.15$, the range of $t_d$ is [0, 20] (see Figure 3D). It shrinks to [0, 14] for $\beta =0.2$ and finally to [1, 9] for $\beta =0.2$ (see Figure 3E,F). We also note that for low transmission rate values, the distribution of $t_d$ is flatter than that for higher transmission rates (see Figure 3). This indicates that estimation of the mobility parameter from observed differences in peak timing ($t_d$) would be associated with greater uncertainty in the scenario where the transmissibility of the infection is comparatively lower.

Now, we explore the feasibility of estimating the mobility flux (m) from observed differences in peak timing ($t_d$). We calculate the maximum likelihood estimate (MLE) of the parameter m, given the observed value of $t_d$, and other model parameters ($\beta ,\gamma ,N_1,N_2,I_1\left(0\right),and{I}_2\left(0\right)$). Here, the observed value of $t_d$ in each case is obtained by integrating the deterministic model (1) with the same set of fixed parameters. Since it is evident from the above analysis that the transmission rate ($\beta$) is influential in calculating $t_d$, for an observed value $t_d$, we consider two different values of $\beta$ and estimate the MLE of the parameter m.

We calculate the log-likelihood (LL) of the parameter m for a range of possible values using the stochastic simulation described in Section 2.2. Then, using cubic splines, we smooth the likelihood profile to obtain the maximum likelihood estimate and the corresponding confidence interval using a likelihood ratio test. For example, we consider a case where the observed value of $t_d$ is 5 days. If we fix $\beta =0.2$, then the MLE of m is $0.031$, and if $\beta$ is taken as $0.3$, the MLE of m is reduced to $0.021$ (see Figure 4A,C). In another example, we consider a higher observed value of $t_d$, i.e., 10 days. In this case, for $\beta =0.2$ the MLE of m is $0.018$ and for $\beta =0.3$ it becomes $0.01$ (see Figure 4B,D). From Figure 4A–D, in all four cases, we can see that MLEs of m (vertical solid lines) are quite close to the true value of m (dashed vertical line), the value of mobility parameter obtained from (1) for the corresponding value of $t_d$.

Note also that for a given $t_d$, the uncertainty in estimating the mobility parameter decreases with increasing disease transmission rate $\beta$. For both the lower and higher observed $t_d$ scenarios, the confidence intervals on m shrink when we increase $\beta$ (see Figure 4A–D). For a fixed value of $\beta$, if the observed value of $t_d$ is higher, then we also observe reduced uncertainty. For example, if $\beta$ is fixed as $0.2$, then the confidence interval on m is wider for $t_d=5$ days than $t_d=10$ (see Figure 4A,B). A similar trend is also observed for higher transmission rates, e.g., $\beta =0.3$ (see Figure 4C,D).

Now, we consider the case when ${\beta }_1\ne {\beta }_2$. The scenario essentially captures the realistic situation in which the intensity of disease transmissibility varies due to different levels of prevention measures being implemented. We explore how the estimate of the mobility parameter changes with the difference in transmission rates. To illustrate this, we consider a case where the observed value of $t_d$ is 5 days. If we fix ${\beta }_1=0.25$, and ${\beta }_2=0.20$ then the MLE of m is $0.076$, and if ${\beta }_1$, and ${\beta }_2$ are taken as $0.30$ and $0.20$, respectively, the MLE of m is increased to $0.108$ (see Figure 5A,C). In another example, we consider a larger observed value of $t_d$, i.e., 10 days. In this case, for ${\beta }_1=0.25$ and ${\beta }_2=0.20$, the MLE of m is $0.038$, and for ${\beta }_1=0.30$ and ${\beta }_2=0.20$, it becomes $0.051$ (see Figure 5B,D). Similar to the case when ${\beta }_1={\beta }_2$, from Figure 5A–D, we can see that the MLEs of m (vertical solid lines) are quite close to the true value of m (dashed vertical line), the value of the mobility parameter obtained from (1) for the corresponding value of $t_d$.

Case Study: COVID-19 Incidence in Poland–Germany Border Region

We apply our proposed method to COVID-19 incidence data during the period from 20 February 2021 to 20 June 2021 for the pair of regions located on either side of the Poland–Germany border directly connected through different modes of transportation, e.g., bus transportation and train transportation (see Figure 6A). For each pair of regions, we obtain the corresponding information for the total population and initial number of infections (see Table 1). All these datasets are obtained from the database pipeline by [38]. We fix the value of the average infectious period ($\frac{1}{\gamma }$) at 14, following [39].

We have previously demonstrated that the disease transmission rate plays an influential role in the estimation of mobility flux given the difference in peak timing. Therefore, in the case of real data, instead of fixing the disease transmission rate arbitrarily, we follow a simpler approach to estimate these values from the initial phase of the epidemic. Since early in a epidemic the daily number of new cases grows exponentially over time, for a SIR-type model the number of infections at time t can be approximated as [40]:

$I\left(t\right)\approx I\left(0\right)exp\left(\lambda t\right).$

Here, $\lambda$ is the growth rate and is given by $\lambda =\beta -\gamma$. Since the cumulative number of new cases per day ($C\left(t\right)$) and new cases per day ($I\left(t\right)$) are linearly related [41,42], i.e., $I\left(t\right)$∼$\lambda$$C\left(t\right)$, we estimate the growth rate $\lambda$ using a linear fit. Therefore, once the parameter $\gamma$ is kept fixed, the disease transmission rate $\beta$ can be approximated as $\beta \approx \lambda +\gamma$. Following this approach, for each region, we estimate the disease transmission rate from the time series data of daily new cases (see Table 1).

We first consider the regions: Gryfiński, and Uckermark, which are located at the upper part of the Poland–Germany border. The observed difference in peak timing obtained from the real data is 20 days (see Figure 6B). Given all the other model parameters, we estimate the mobility flux as 0.038 [0.022, 0.064] (see Figure 6C). Next, we move to the middle part of the border and choose the regions: Słubicki and Frankfurt(Oder). In this case, the observed difference in peak timing is 10 days (see Figure 6D). The maximum likelihood estimate of m is 0.032 [0.004, 0.052] (see Figure 6E). Now, for the regions Zgorzelecki and Görlitz, the observed difference in peak timing is 20 days (see Figure 6F). The estimated value of the parameter m in this case is 0.012 [0.002, 0.018] (see Figure 6G).

Comparing the first case and third case, we see that there is a substantial difference between the estimates of mobility parameters even though the observed differences in peak time ($t_d$) for both cases are the same. Instead, the difference in estimated transmission rates (i.e., ${\beta }_1$ and ${\beta }_2$) is higher in the first case (Gryfiński and Uckermark) than that of the third case (Zgorzelecki and Görlitz) (see Table 1). This marked difference in transmission rates likely explains the different estimates of the coupling parameter.

4. Discussion

In general, understanding the interplay between human mobility and spatiotemporal propagation of infectious diseases is a long-standing problem in the field of infectious disease modeling. The unavailability of reliable information on human mobility poses challenges to designing effective public health strategies. In this study, we proposed a method to address the issue of estimating cross-country mobility flux during a disease outbreak, which is important yet less explored in the literature. We showed that from the observed difference in peak timing between any two regions that share a border, the underlying mobility flow between them can be retrieved effectively. Since the information on peak timing is less affected than the invasion time by the delay or incompleteness in reporting infection, the difference in peak timing can be considered as a more robust quantity than invasion time. Consequently, our approach would be more broadly applicable for estimating mobility flux in real-world case studies than existing alternatives that rely on the observed disease invasion time [29]. It is noteworthy that our approach, relying on minimal data requirements, can be utilized to understand the degree and impact of human mobility on disease spread. This is particularly valuable for rapid assessments during the initial phase of an outbreak when fine-grained data are unavailable.

We investigated how different model parameters influence the difference in peak timing in a simple scenario, where the disease transmission rates in both regions are the same. We see that the disease transmission rate plays a key role in determining the distribution of differences in peak timing. The difference in peak timing decreases with the disease transmission rate (Figure 2). This is because of the fact that a higher disease transmission rate causes the earlier occurrence of the peak in infection trajectory and thus the difference in peak timing between the two regions also decreases. It has also been revealed that, for lower transmission rates, the distribution of differences in peak timing becomes flatter and can have a wider range of possible values (Figure 3A,D). On the other hand, with increasing transmission rates the distribution progressively narrows (Figure 3B,C,E,F).

Based on a simple SIR-type model, with the help of stochastic simulations, we provided a maximum likelihood estimate of the cross-border mobility under different epidemiologically relevant scenarios. From our numerical investigation, it is evident that the estimation method based on the observed difference in peak timing can effectively recover the underlying mobility rate for both scenarios: (i) when the disease transmission rates in both regions are assumed to be the same (Figure 4A–D) and (ii) when the transmission rates are assumed to be different (Figure 5A–D). In the first scenario, for a given difference in peak-timing, the uncertainty of the estimate decreases with the transmission rate (Figure 4A–D). Also for a fixed transmission rate, the uncertainty in estimation decreases with increasing values of the observed difference in peak timing (Figure 4A–D). In the second scenario, for a given observed difference in peak timing, the estimate of the mobility parameter increases with the difference between the disease transmission rates (Figure 5A–D). Similar to the first scenario, for fixed disease transmission rates, the uncertainty in the estimate decreases with the observed difference in peak timing (Figure 5A–D).

To demonstrate its real-world utility, we applied our proposed method to COVID-19 incidence data obtained for the regions located along the border of Poland and Germany. We chose only three pairs of regions located in different parts of the Polish–German border which are well connected through different modes of transportation. However, the method can also be applied to any other pair of regions that share a border. Using least squares techniques, we first estimated the disease transmission rate for each of the regions from the early phase of epidemics. We then calculated the maximum likelihood estimate of the mobility parameter for three pairs of regions located in different areas along the border. It is noteworthy that if we compare the pair of regions Gryfiński, Uckermark, Zgorzelecki, and Görlitz, we see that, even though the differences in peak timing are the same, the estimated values for the cross-border mobility are very different. Since the difference in transmission rates between the regions Gryfiński and Uckermark is higher than that of Zgorzelecki and Görlitz, we see higher mobility flow for the regions Gryfiński and Uckermark. This result is in good agreement with our theoretical findings (Figure 5).

Even though our study addresses the important issue of estimating cross-border mobility and our proposed method worked well in the real case study, it has some limitations. We assumed that the mobility between two regions is symmetric, which is not generally true, particularly when considering mobility between a large city and a small village. This assumption simplifies the model but may not accurately reflect the asymmetrical nature of real-world human movement patterns.

Additionally, our model considered the regions on either side of the border as isolated, ignoring the contributions from adjacent regions. In reality, disease transmission is influenced by a network of interactions among multiple regions, and ignoring these interactions can lead to an incomplete understanding of the mobility flux and its impact on disease spread. Addressing this issue requires a more comprehensive model that includes a fixed number of regions connected through explicit, not necessarily symmetric, migration patterns.

Another limitation of our approach is the reliance on the Susceptible–Infected–Recovered (SIR) model, which, while useful for many infectious diseases, may not capture the complexity of diseases with more intricate transmission dynamics, such as those with significant latent periods or multiple stages of infection. Extending our method to incorporate more complex epidemiological models could enhance its applicability to a broader range of diseases.

Despite these limitations, our study provides a foundational approach to estimating cross-border mobility using readily available data, which is a significant step forward in epidemiological modeling. Future work could focus on developing and testing more general models that account for asymmetric mobility and interactions among multiple regions. Such models could provide more accurate estimates of mobility flux and improve our ability to predict and control the spread of infectious diseases across borders. Investigating alternative statistical techniques or incorporating machine learning approaches to handle such data imperfections could be a promising direction for future research.

We leave these modifications and extensions for future study, acknowledging that addressing these challenges is crucial for advancing the field of spatial epidemiology and improving our preparedness for managing cross-border infectious disease outbreaks.

Summing up, this work provides a simple and intuitive method for estimating mobility flow between the regions connected through human migration. The methodology and findings can serve as a basis for designing intervention strategies in situations where direct information on human mobility is inaccessible.

Author Contributions

Conceptualization, A.S. and J.M.C.; methodology, A.S.; software, A.S.; validation, A.S. and J.M.C.; formal analysis, A.S.; investigation, J.M.C.; data curation, A.M. and W.S.-W.; writing—original draft preparation, A.S.; writing—review and editing, A.S., J.M.C., A.M. and W.S.-W.; visualization, A.S.; supervision, J.M.C.; funding acquisition, J.M.C. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The Python codes for generating the figures presented in the paper are available in the GitHub repository https://github.com/abhisheksena/cross_border_mobility (accessed on 1 May 2024).

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. (A) Example of stochastic realizations of the model for a fixed combination of model parameters. The red and blue trajectories refer to [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.], respectively. The dashed red and blue line indicate the points at which [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] attain their maxima, respectively. The difference between these time points is shown as [Forumla omitted. See PDF.]. The unit for time is in days. (B) Histogram of [Forumla omitted. See PDF.] for 10,000 stochastic realizations.

Figure 2. (A–F) Behavior of the difference in peak timing ([Forumla omitted. See PDF.]) with respect to coupling strength (m) for different combinations of transmission rate and initial infection. Here, the solid line represents the mean of the stochastic realizations and the shaded region represents a [Forumla omitted. See PDF.] confidence interval around the mean. The subfigures (B,C,E,F) share the same vertical axis label as (A,D).

Figure 3. (A–F) Distribution of [Forumla omitted. See PDF.] corresponding to the mobility parameter [Forumla omitted. See PDF.], for different parameter combinations as mentioned in Figure 2. The subfigures (B,C,E,F) share the same vertical axis label as (A,D).

Figure 4. (A–D) Maximum likelihood estimates of the mobility parameter (m) given the difference in peak timing ([Forumla omitted. See PDF.]). The solid circles are the log-likelihood obtained from stochastic simulations and the solid curve is the smoothed likelihood profile. The solid vertical line is the maximum likelihood estimate of m and the dashed vertical line is the true value of m, for which the deterministic model (1) gives the corresponding observed value of [Forumla omitted. See PDF.]. The shaded region depicts the [Forumla omitted. See PDF.] confidence interval. The remaining parameters are fixed as [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] = 10,000, [Forumla omitted. See PDF.] = 20,000, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.].

Figure 5. (A–D) Maximum likelihood estimate of the mobility parameter (m) given the difference in peak timing ([Forumla omitted. See PDF.]), when [Forumla omitted. See PDF.]. The solid circles are the log-likelihood obtained from stochastic simulation and the solid curve is the smoothed likelihood profile. Solid vertical line is the maximum likelihood estimate of m and the dashed vertical line is the true value of m, for which the deterministic model (1) gives the corresponding observed value of [Forumla omitted. See PDF.]. The shaded region depicts the [Forumla omitted. See PDF.] confidence interval. The remaining parameters are fixed as [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] = 10,000, [Forumla omitted. See PDF.] = 20,000, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.].

Figure 6. (A) The map shows the three pairs of regions that share a border between Poland and Germany. The Polish regions (called Powiat) are filled with dark blue color, and the German regions (called Kreise) are shaded dark green. The name of the regions are as follows: top—(Gryfiński, Uckermark); middle—(Słubicki, Frankfurt (Oder)); and bottom—(Zgorzelecki, Görlitz). (B,D,F) The time series of the number of infecteds in the selected regions. The blue lines represent the time series for the Polish regions and green lines represent the time series for the German regions. Here, time zero refers to the date 20 February 2021 and time 120 refers to the date 20 June 2021. (C,E,G) The likelihood profiles of the cross-border mobility m. The gray circles are the log-likelihoods obtained from the stochastic simulation and the solid curve is the smoothed likelihood profile. The dashed vertical line is the maximum likelihood estimate of m and the shaded region denotes the 75% confidence interval.

Table 1

Parameters and observed difference in peak timing ($t_d$) for each pair of region used for estimating the mobility flux m. The estimated value of m is presented in the last column of the table.

Pair (Powiat, Kreise)	Transmission Rate (${\mathit{\beta }}_1$, ${\mathit{\beta }}_2$)	Total Population (${\mathit{N}}_1$, ${\mathit{N}}_2$)	Initial Infection $({\mathit{I}}_1\left(0\right),{\mathit{I}}_2\left(0\right))$	Observed Difference in Peak Timing (${\mathit{t}}_{\mathit{d}}$)	Estimated Mobility Flux ($\mathit{m}[75\%\mathit{CI}]$)
(Gryfiński, Uckermark)	(0.10, 0.075)	(82,951, 119,552)	(11, 8)	20 days	0.038 [0.022, 0.064]
(Słubicki, Frankfurt(Oder))	(0.13, 0.14)	(47,068, 57,873)	(8, 1)	10 days	0.032 [0.004, 0.052]
(Zgorzelecki, Görlitz)	(0.10, 0.11)	(90,584, 254,894)	(9, 23)	20 days	0.012 [0.002, 0.018]

References

1. Brockmann, D. Human mobility and spatial disease dynamics. Rev. Nonlinear Dyn. Complex.; 2009; 2, pp. 1-24.

2. Findlater, A.; Bogoch, I.I. Human mobility and the global spread of infectious diseases: A focus on air travel. Trends Parasitol.; 2018; 34, pp. 772-783. [DOI: https://dx.doi.org/10.1016/j.pt.2018.07.004] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30049602]

3. Alessandretti, L.; Aslak, U.; Lehmann, S. The scales of human mobility. Nature; 2020; 587, pp. 402-407. [DOI: https://dx.doi.org/10.1038/s41586-020-2909-1] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33208961]

4. Barbosa, H.; Barthelemy, M.; Ghoshal, G.; James, C.R.; Lenormand, M.; Louail, T.; Menezes, R.; Ramasco, J.J.; Simini, F.; Tomasini, M. Human mobility: Models and applications. Phys. Rep.; 2018; 734, pp. 1-74. [DOI: https://dx.doi.org/10.1016/j.physrep.2018.01.001]

5. Schläpfer, M.; Dong, L.; O’Keeffe, K.; Santi, P.; Szell, M.; Salat, H.; Anklesaria, S.; Vazifeh, M.; Ratti, C.; West, G.B. The universal visitation law of human mobility. Nature; 2021; 593, pp. 522-527. [DOI: https://dx.doi.org/10.1038/s41586-021-03480-9]

6. Wilder-Smith, A.; Gubler, D.J. Geographic expansion of dengue: The impact of international travel. Med. Clin. N. Am.; 2008; 92, pp. 1377-1390. [DOI: https://dx.doi.org/10.1016/j.mcna.2008.07.002]

7. Cavallaro, F.; Dianin, A. Cross-border commuting in Central Europe: Features, trends and policies. Transp. Policy; 2019; 78, pp. 86-104. [DOI: https://dx.doi.org/10.1016/j.tranpol.2019.04.008]

8. Docquier, F.; Golenvaux, N.; Nijssen, S.; Schaus, P.; Stips, F. Cross-border mobility responses to COVID-19 in Europe: New evidence from facebook data. Glob. Health; 2022; 18, 41. [DOI: https://dx.doi.org/10.1186/s12992-022-00832-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35436927]

9. Shyamsunder,; Bhatter, S.; Jangid, K.; Abidemi, A.; Owolabi, K.M.; Purohit, S.D. A new fractional mathematical model to study the impact of vaccination on COVID-19 outbreaks. Decis. Anal. J.; 2023; 6, 100156.

10. Kumawat, S.; Bhatter, S.; Suthar, D.; Purohit, S.; Jangid, K. Numerical modeling on age-based study of coronavirus transmission. Appl. Math. Sci. Eng.; 2022; 30, pp. 609-634. [DOI: https://dx.doi.org/10.1080/27690911.2022.2116435]

11. Shubin, M.; Brustad, H.K.; Midtbo, J.E.; Gunther, F.; Alessandretti, L.; Ala-Nissila, T.; Scalia Tomba, G.; Kivela, M.; Yat Hin Chan, L.; Leskela, L. The influence of cross-border mobility on the COVID-19 epidemic in Nordic countries. medRxiv; 2023; 2023–11. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1012182] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38865414]

12. Kühn, M.J.; Abele, D.; Binder, S.; Rack, K.; Klitz, M.; Kleinert, J.; Gilg, J.; Spataro, L.; Koslow, W.; Siggel, M. et al. Regional opening strategies with commuter testing and containment of new SARS-CoV-2 variants in Germany. BMC Infect. Dis.; 2022; 22, 333. [DOI: https://dx.doi.org/10.1186/s12879-022-07302-9] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35379190]

13. Balcan, D.; Colizza, V.; Gonçalves, B.; Hu, H.; Ramasco, J.J.; Vespignani, A. Multiscale mobility networks and the spatial spreading of infectious diseases. Proc. Natl. Acad. Sci. USA; 2009; 106, pp. 21484-21489. [DOI: https://dx.doi.org/10.1073/pnas.0906910106] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20018697]

14. Jandarov, R.; Haran, M.; Bjørnstad, O.; Grenfell, B. Emulating a gravity model to infer the spatiotemporal dynamics of an infectious disease. J. R. Stat. Soc. Ser. (Appl. Stat.); 2014; 63, pp. 423-444. [DOI: https://dx.doi.org/10.1111/rssc.12042]

15. Xia, Y.; Bjørnstad, O.N.; Grenfell, B.T. Measles metapopulation dynamics: A gravity model for epidemiological coupling and dynamics. Am. Nat.; 2004; 164, pp. 267-281. [DOI: https://dx.doi.org/10.1086/422341] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15278849]

16. Senapati, A.; Sardar, T.; Ganguly, K.S.; Ganguly, K.S.; Chattopadhyay, A.K.; Chattopadhyay, J. Impact of adult mosquito control on dengue prevalence in a multi-patch setting: A case study in Kolkata (2014–2015). J. Theor. Biol.; 2019; 478, pp. 139-152. [DOI: https://dx.doi.org/10.1016/j.jtbi.2019.06.021] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31229456]

17. Mari, L.; Bertuzzo, E.; Righetto, L.; Casagrandi, R.; Gatto, M.; Rodriguez-Iturbe, I.; Rinaldo, A. Modelling cholera epidemics: The role of waterways, human mobility and sanitation. J. R. Soc. Interface; 2012; 9, pp. 376-388. [DOI: https://dx.doi.org/10.1098/rsif.2011.0304]

18. Yang, W.; Zhang, W.; Kargbo, D.; Yang, R.; Chen, Y.; Chen, Z.; Kamara, A.; Kargbo, B.; Kandula, S.; Karspeck, A. et al. Transmission network of the 2014–2015 Ebola epidemic in Sierra Leone. J. R. Soc. Interface; 2015; 12, 20150536. [DOI: https://dx.doi.org/10.1098/rsif.2015.0536] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26559683]

19. Mertel, A.; Vyskočil, J.; Schüler, L.; Schlechte-Wełnicz, W.; Calabrese, J.M. Fine-scale variation in the effect of national border on COVID-19 spread: A case study of the Saxon-Czech border region. Spat. Spatio-Temporal Epidemiol.; 2023; 44, 100560. [DOI: https://dx.doi.org/10.1016/j.sste.2022.100560]

20. Simini, F.; González, M.C.; Maritan, A.; Barabási, A.L. A universal model for mobility and migration patterns. Nature; 2012; 484, pp. 96-100. [DOI: https://dx.doi.org/10.1038/nature10856]

21. Timme, M. Revealing network connectivity from response dynamics. Phys. Rev. Lett.; 2007; 98, 224101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.98.224101] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17677845]

22. Casadiego, J.; Nitzan, M.; Hallerberg, S.; Timme, M. Model-free inference of direct network interactions from nonlinear collective dynamics. Nat. Commun.; 2017; 8, 2192. [DOI: https://dx.doi.org/10.1038/s41467-017-02288-4] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29259167]

23. Shandilya, S.G.; Timme, M. Inferring network topology from complex dynamics. New J. Phys.; 2011; 13, 013004. [DOI: https://dx.doi.org/10.1088/1367-2630/13/1/013004]

24. Gao, T.T.; Yan, G. Autonomous inference of complex network dynamics from incomplete and noisy data. Nat. Comput. Sci.; 2022; 2, pp. 160-168. [DOI: https://dx.doi.org/10.1038/s43588-022-00217-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38177441]

25. Ching, E.S.; Lai, P.Y.; Leung, C. Reconstructing weighted networks from dynamics. Phys. Rev.; 2015; 91, 030801. [DOI: https://dx.doi.org/10.1103/PhysRevE.91.030801] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25871037]

26. Shi, L.; Shen, C.; Jin, L.; Shi, Q.; Wang, Z.; Boccaletti, S. Inferring network structures via signal Lasso. Phys. Rev. Res.; 2021; 3, 043210. [DOI: https://dx.doi.org/10.1103/PhysRevResearch.3.043210]

27. Alves, C.L.; Pineda, A.M.; Roster, K.; Thielemann, C.; Rodrigues, F.A. EEG functional connectivity and deep learning for automatic diagnosis of brain disorders: Alzheimer’s disease and schizophrenia. J. Phys. Complex.; 2022; 3, 025001. [DOI: https://dx.doi.org/10.1088/2632-072X/ac5f8d]

28. Sysoeva, M.V.; Sysoev, I.V.; Prokhorov, M.D.; Ponomarenko, V.I.; Bezruchko, B.P. Reconstruction of coupling structure in network of neuron-like oscillators based on a phase-locked loop. Chaos Solitons Fractals; 2021; 142, 110513. [DOI: https://dx.doi.org/10.1016/j.chaos.2020.110513]

29. Hempel, K.; Earn, D.J. Estimating epidemic coupling between populations from the time to invasion. J. R. Soc. Interface; 2020; 17, 20200523. [DOI: https://dx.doi.org/10.1098/rsif.2020.0523]

30. Bichara, D.; Kang, Y.; Castillo-Chavez, C.; Horan, R.; Perrings, C. SIS and SIR epidemic models under virtual dispersal. Bull. Math. Biol.; 2015; 77, pp. 2004-2034. [DOI: https://dx.doi.org/10.1007/s11538-015-0113-5]

31. Bichara, D.; Iggidr, A. Multi-patch and multi-group epidemic models: A new framework. J. Math. Biol.; 2018; 77, pp. 107-134. [DOI: https://dx.doi.org/10.1007/s00285-017-1191-9]

32. Espinoza, B.; Castillo-Chavez, C.; Perrings, C. Mobility restrictions for the control of epidemics: When do they work?. PLoS ONE; 2020; 15, e0235731. [DOI: https://dx.doi.org/10.1371/journal.pone.0235731] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32628716]

33. Mhlanga, A.; Mupedza, T. A patchy theoretical model for the transmission dynamics of SARS-cov-2 with optimal control. Sci. Rep.; 2022; 12, 17840. [DOI: https://dx.doi.org/10.1038/s41598-022-21553-1] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36284219]

34. Agusto, F.; Goldberg, A.; Ortega, O.; Ponce, J.; Zaytseva, S.; Sindi, S.; Blower, S. How do interventions impact malaria dynamics between neighboring countries? A case study with Botswana and Zimbabwe. Using Mathematics to Understand Biological Complexity: From Cells to Populations; Association for Women in Mathematics Series Springer: Cham, Switerland, 2021; Volume 22, pp. 83-109. [DOI: https://dx.doi.org/10.1007/978-3-030-57129-0_5]

35. Gillespie, D.T. Approximate accelerated stochastic simulation of chemically reacting systems. J. Chem. Phys.; 2001; 115, pp. 1716-1733. [DOI: https://dx.doi.org/10.1063/1.1378322]

36. Gillespie, D.T.; Petzold, L.R. Improved leap-size selection for accelerated stochastic simulation. J. Chem. Phys.; 2003; 119, pp. 8229-8234. [DOI: https://dx.doi.org/10.1063/1.1613254]

37. Linka, K.; Peirlinck, M.; Kuhl, E. The reproduction number of COVID-19 and its correlation with public health interventions. Comput. Mech.; 2020; 66, pp. 1035-1050. [DOI: https://dx.doi.org/10.1007/s00466-020-01880-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32836597]

38. Abdussalam, W.; Mertel, A.; Fan, K.; Schüler, L.; Schlechte-Wełnicz, W.; Calabrese, J.M. A scalable pipeline for COVID-19: The case study of Germany, Czechia and Poland. arXiv; 2022; arXiv: 2208.12928

39. Byrne, A.W.; McEvoy, D.; Collins, A.B.; Hunt, K.; Casey, M.; Barber, A.; Butler, F.; Griffin, J.; Lane, E.A.; McAloon, C. et al. Inferred duration of infectious period of SARS-CoV-2: Rapid scoping review and analysis of available evidence for asymptomatic and symptomatic COVID-19 cases. BMJ Open; 2020; 10, e039856. [DOI: https://dx.doi.org/10.1136/bmjopen-2020-039856] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32759252]

40. Ma, J. Estimating epidemic exponential growth rate and basic reproduction number. Infect. Dis. Model.; 2020; 5, pp. 129-141. [DOI: https://dx.doi.org/10.1016/j.idm.2019.12.009]

41. Favier, C.; Dégallier, N.; Rosa-Freitas, M.G.; Boulanger, J.P.; Costa Lima, J.; Luitgards-Moura, J.F.; Menkes, C.E.; Mondet, B.; Oliveira, C.; Weimann, E. et al. Early determination of the reproductive number for vector-borne diseases: The case of dengue in Brazil. Trop. Med. Int. Health; 2006; 11, pp. 332-340. [DOI: https://dx.doi.org/10.1111/j.1365-3156.2006.01560.x]

42. Pinho, S.T.R.d.; Ferreira, C.P.; Esteva, L.; Barreto, F.R.; Morato e Silva, V.; Teixeira, M. Modelling the dynamics of dengue real epidemics. Philos. Trans. R. Soc. Math. Phys. Eng. Sci.; 2010; 368, pp. 5679-5693. [DOI: https://dx.doi.org/10.1098/rsta.2010.0278] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21078642]