Using hybrid automata to model mitigation of global disease spread via travel restriction

1.   Introduction

Mathematical modeling of disease spread became more important than ever during the COVID-19 pandemic. Various models were used in an attempt to forecast possible outcomes of the pandemic, often in order to help political leaders make informed decisions about implementing policies intended to mitigate the potential negative impacts of the spreading disease ^[1]. In that regard, most models were solely focused on making local predictions in whatever particular geographical location was of primary concern: a state, country, or larger region, for example. Notably lacking is a global model designed to better understand how many interacting regions, all simultaneously implementing different mitigation policies, could possibly effect each other and impact disease spread on a larger scale. The Federal Aviation Administration alone accounts for $2.9$ million airline passengers each day ^[2]. Other estimates list $4.5$ billion flights in 2019 prior to the pandemic ^[3]. In this manner, a virus could quickly cover vast distances, and thereby couple the local disease dynamics of two regions between which air travel exists. Thus, the world wide network of air traffic has the potential to control or effect the global impact of a pandemic.

In this paper, we focus on the mitigation of global disease spread via air travel restrictions. Geographic locations have their own local dynamics for disease spread, which vary continuously. Travel restriction policies, even if updated and implemented regularly, naturally vary in a discrete manner and potentially alter the dynamics of the spread of the disease across multiple regions. The general formulation of the problem we are interested in is of dynamical systems that exist within a larger network, i.e., dynamical systems for which we have a set of many (potentially different) kinds of local dynamics, coupled together via a number of pairwise connections. The system can be written in a compact form as

where $x = [x_1\quad x_2\cdots x_n]^T$ is in $\mathbb{R}^{mn}$ with each $x_i\in \mathbb{R}^m$. The function $g(x) = [g_1(x_1), \dots, g_n(x_n)]^T$ distinguishes the local dynamics from the coupling between different local dynamics obtained by multiplication by $\mathcal{L}$, which is a matrix representation of a graph Laplacian. Here, $n$ is the number of different local dynamics and $m$ is the number of states for each node representing local dynamics. The local dynamics can be associated with vertices in a graph representing the network, with edges (which may potentially be weighted and/or directed) modeling connections between them, as shown in Figure 1.

For our purposes, the coupling is determined by the graph Laplacian, whose matrix representation has particularly desirable properties ^[4,5]. In addition, we are interested in systems that can evolve continuously in time, but further allow for discontinuities in the dynamics when some part of the system has reached a threshold that fundamentally changes the way the system will evolve subsequently ^[6,7]. Of primary interest are discontinuities in the dynamics that effect the coupling from $\mathcal{L}$ but leave the local dynamics $g(x)$ unaffected. The superstructure itself can also then be viewed as a graph, with a vertex representing a particular dynamical system, and outgoing edges encoding the threshold values that may potentially be reached, thereby leading to an entirely different instance of a dynamical system. We illustrate an example in Figure 2.

For our application, referring to Eq (1.1), we have that $g(x)$ describes the local dynamics for disease spread, for example, representing a compartmental Susceptible-Exposed-Infected-Recovered (SEIR) model and $\mathcal{L}$ models the coupling between regions due to potentially infected individuals traveling between them. As mentioned before, modeling during the COVID-19 pandemic was most often done locally, albeit out of the necessity of immediacy, but also due to problems of scalability. Indeed, a compartmental model like a SEIR model intended to accommodate the entire world population would not automatically account, for example, for variants of a virus that are active in different regions, or for other details that are specific locally. Additionally, a compartmental model would not inherently allow for the mixing of populations via travel, while factoring in travel restrictions that vary from region to region. On the other hand, an agent-based model at the scale of the world is computationally unfeasible. In response, we aim to provide a proof of concept for a hybrid automaton model that functions at the scale of the entire world and regulates travel in a deterministic way via so-called "guard conditions".

The outline of this paper is as follows. In Section 2, we develop carefully our epidemiological model. Specifically in Section 2.1, we provide the general formalism and definitions for hybrid automata. Next, in the five parts of Section 2.2, we describe in detail each of the individual components in the definition relative to the particular hybrid automaton, which we will use to model the mitigated spread of disease. In Section 3, we implement the model and run simulations. Finally, in Section 4 we summarize our results and note goals for future work.

2.   Hybrid automaton for disease spread

A typical approach in epidemiological modeling is compartmental models written in terms of differential or difference equations ^{[8,9,10,11,12]}. Another is agent-based models, whose complexity is both a drawback, analytically and computationally, and a strength, as they can capture more nuanced behaviors by tracking individuals separately ^[13]. Our modeling efforts for the state of Hawai'i included work using both approaches and provided our political leaders with insights, while the spreading disease required them to make difficult decisions about mitigation policies ^[14]. However, neither modeling approach is perfectly suited for developing a worldwide model. Indeed, a single compartmental model does not automatically include travelers moving between regions, a potentially significant effect to neglect, as was made clear in modeling the islands of Hawai'i. On the other hand, using an agent-based model to track individual agents and their movements is computationally too expensive to handle a global population of more than $8$ billion people, and, most importantly, would require the construction of a synthetic global population.

Of all of the mitigation policies implemented in the state of Hawai'i throughout the COVID-19 pandemic, perhaps the most drastic were travel restrictions. These included a quarantine period required for all Spring 2020 travelers, which temporarily stopped the important economy of tourism, or later requiring only proof of vaccination, which helped restore the state's tourism rate prior to the pandemic ^[14]. Inspired by the aforementioned work and the potentiality of various types of travel restrictions, we present an innovative global model that accounts for coordination between regions, each using travel restrictions in order to mitigate the worldwide spread of disease. The model is designed via a family of compartmental SEIR models and using the formalism of hybrid automata.

2.1.   Hybrid automaton

In this section, we define formally the notion of a hybrid automaton. Equivalent definitions are given in ^[6,7], but we prefer the treatment found in ^[15].

Definition 2.1. A hybrid automaton $H$ is given by the following collection of data: $H = (Q, X, f, {I}nit, {D}om, E, {G}rd, {R}es)$, where:

● $Q = \{q_1, q_2, \dots\}$ is a set of discrete states (a finite or countable number of them),

● $X\subseteq\mathbb{R}^d$ is a set of continuous states,

● $f(\cdot, \cdot):Q\times X\to\mathbb{R}^d$ is a vector field,

● ${I}nit\subseteq Q\times X$ is a set of initial states,

● ${D}om(\cdot):Q\to 2^X$ is a domain (assigns to each discrete state $q\in Q$ a set of continuous states ${D}om(q)\subseteq\mathbb{R}^d$ that make sense for that particular $q$),

● $E\subseteq Q\times Q$ is a set of (directed) edges (encoding the eligible changes the discrete state may make),

● ${G}rd(\cdot):E\to 2^X$ is a guard condition, (explained subsequently)

● ${R}es(\cdot, \cdot):E\times X\to 2^X$ is a reset map (explained subsequently).

Definition 2.2. We refer to $(q, x)\in Q\times X$ as the state of the hybrid automaton $H$.

Hybrid automata necessarily involve both continuous flow determined by the vector field $f(q, \cdot)$ and also discrete changes to said flow, determined by the edges of a digraph. To describe this in more detail, we desire a set of times that include both continuous intervals over which the continuous evolution occurs and particular distinguished moments in time when the transitions along an edge occur. As such, we introduce a notion of a hybrid time set.

Definition 2.3. A hybrid time set is a sequence of intervals $\tau = \{I_0, I_1, \dots, I_K\}$ (either a finite number $K < \infty$ or an infinite number $K = \infty$), such that

● $I_k = [\tau_k, \tau_k']$ for each $k < K$,

● if $K < \infty$, then either $I_K = [\tau_K, \tau_K']$ or $I_K = [\tau_K, \tau_K')$,

● $\tau_k\le\tau_k' = \tau_{k+1}$ for all $k\ge 0$.

The idea is that our primed times are the times at which discrete transitions of the hybrid system take place. In other words, the rightmost point of the $I_k$ interval, $\tau_k'$, corresponds to the time instant just before a discrete transition, the leftmost point of the $I_{k+1}$ interval, $\tau_{k+1}$, corresponds to the time instant just after the discrete transition, and that $\tau_k' = \tau_{k+1}$ corresponds to the fact that we assume discrete transitions happen instantaneously. One advantage of this convention is that it easily affords multiple discrete transitions taking place one after the other. This is because we allow specifically in the definition degenerate intervals $\tau_k\le\tau_k'$ instead of imposing $\tau_k < \tau_k'$, so situations like $\tau_{k-1}' = \tau_k = \tau_k' = \tau_{k+1}$ can very well occur.

We can now precisely define trajectories of our hybrid system and, thus, elaborate on the meaning of the guard condition and reset map in Definition 2.1.

Definition 2.4. A hybrid trajectory is a triple $(\tau, \hat{q}, \hat{x})$ consisting of a hybrid time set $\tau = \{I_k\}_{k = 0}^{K}$ and two sequences of functions $\hat{q} = \{q_k(\cdot)\}_{k = 0}^{K}$, with each $q_k(\cdot):I_k\to Q$ and $\hat{x} = \{x_k(\cdot)\}_{k = 0}^{K}$, with each $x_k(\cdot):I_k\to X$.

Definition 2.5. An execution of a hybrid automaton $H$ is a hybrid trajectory $(\tau, \hat{q}, \hat{x})$, which satisfies the following conditions:

● Initial Condition: $(q_0(0)), x_0(0))\in{I}nit$

● Discrete Evolution: For all $k$, we have that

- $(q_k(\tau_k'), q_{k+1}(\tau_{k+1}))\in E$

- $x_k(\tau_k')\in{G}rd(q_k(\tau_k'), q_{k+1}(\tau_{k+1}))$

- $x_{k+1}(\tau_{k+1})\in{R}es(q_k(\tau_k'), q_{k+1}(\tau_{k+1}), x_k(\tau_k'))$

● Continuous Evolution: For all $k$, we have that

- $q_k(\cdot):I_k\to Q$ is constant over $t\in I_k$

- $x_k(\cdot):I_k\to X$ is the solution to the differential equation: $\frac{dx_k(t)}{dt} = f(q_k(t), x_k(t))$ over $t\in I_k$ starting at $x_k(\tau_k)$

- $x_k(t)\in{D}om(q_k(t))$ for all $t\in[\tau_k, \tau_k')$

The definition above specifies which hybrid trajectories are executions of $H$ and which are not. The first restriction simply states that the execution should start at an acceptable initial state in our set $\text{I}nit\subseteq Q\times X$. We typically write $(q_0, x_0): = (q_0(\tau_0), x_0(\tau_0))$ for simplicity. Moreover, we will usually make the assumption that $\tau_0 = 0$, without loss of generality. The second restriction states when discrete transitions can occur and what states are possible after discrete transitions. It relates the state before the discrete transition, $(q_k(\tau_k'), x_k(\tau_k'))$, to the state after the discrete transition, $(q_{k+1}(\tau_{k+1}), x_{k+1}(\tau_{k+1}))$, by noting that $(q_k(\tau_k'), q_{k+1}(\tau_{k+1}))$ should be an existing directed edge, $x_k(\tau_k')$ should belong to the guard condition image of this edge, and $x_{k+1}(\tau_{k+1})$ should belong to the reset image of this edge. Importantly, we think of the guard image as enabling the discrete transition that potentially occurs, i.e., the execution may take a discrete transition so long as the continuous value is in the corresponding guard image (or it may not, but the physical meaning is that it is a physically valid potential transition). The third restriction determines what happens during continuous evolution, and when continuous evolution must give way to a discrete transition. Furthermore, we require that along continuous evolution, the continuous state remains in the domain assigned to the discrete state, $\text{D}om(q)\subseteq\mathbb{R}^d$. More elaborate definitions, constructions, and examples can be found in ^[6,7,15].

2.2.   Hybrid automaton for SEIR model

Our goal is to construct a global model capable of accounting for various travel restrictions between different regions of the world. The spread of disease in each region is modeled as a continuous process, while the coupling between regions is incorporated via a directed graph whose weighted edges depict travel between them. We assume that travel between different regions is done primarily by plane, for instance, regions could be islands or continents. Conceptually, the model would still be valid in the case of regions which are adjacent by land, for example, the $48$ states in the contiguous United States, but it would be assuming control along entire borders in order to enforce the travel restrictions, which is clearly unrealistic. During a pandemic, as it was seen with COVID-19 in Hawai'i, travel restrictions are subject to change based on government policy. They may be more severe at times and less at others, but the change that is implemented is a discrete phenomenon. In this section, we define precisely both the continuous dynamical system and the space on which it is defined, as well as the discrete state of the hybrid system, including the necessary guard conditions.

2.2.1.   SEIR model's continuous states

Suppose we are in a world with a total population of $N$ individuals separated into $n$ disjoint regions that have been labeled uniquely from $1$ to $n$. The continuous portion of the hybrid automata is a SEIR compartmental model, which is the epidemiological model within each region used to model the spread of disease. This choice is motivated by the fact that the complexity of coupling the regional SEIR models through their travelers is not too high and because we want to be able to run the hybrid automaton on a population equaling the current world population of more than 8 billion people.

Consider for each region a SEIR model which possesses compartments for both unvaccinated and vaccinated individuals. This distinction is important since, for example, the vaccine may have the effect of different parameters for transmission, probability to develop severe symptoms, and probability of recovery. In addition, some travel restrictions might make a distinction between vaccinated and unvaccinated individuals. For an arbitrary region $i$, define the vector

where the subscript $i$ refers to the region number, $S$ means "susceptible" individuals, $E$ means "exposed", i.e., asymptomatic infected individuals, $I$ means "infected", i.e., symptomatic infected individuals, $R$ means "removed" and includes those who have either recovered or died from the disease, and the superscripts $u$ and $v$ refer to "unvaccinated" and "vaccinated" status, respectively. We then put all of the regions together to obtain

Altogether, these form our hybrid automaton's set of continuous states, $X\subseteq\mathbb{R}^d$ with dimension $d = 8n$, as in the Definition 2.1.

We assume for all values $t\ge 0$ that we have $S_i^u(t)+E_i^u(t)+I_i^u(t)+R_i^u(t)+S_i^v(t)+E_i^v(t)+I_i^v(t)+R_i^v(t) = N_i$, where $N_i$ is the population of region $i$, currently fixed in this description until we include travelers subsequently. In that regard, note that we neglect to include natural birth and natural death rates, as they are assumed to either balance or to be negligible relative to the spread of disease. We have eight equations (for simplicity of notation here, we have removed the regional dependence denoted by the subscript $i$):

where the following nine parameters dictate transfer between compartments:

● $\theta$ denotes the vaccination rate,

● $\beta^u$ denotes the rate by which the disease is transmitted from the unvaccinated asymptomatic to the susceptible group upon contact,

● $\beta^v$ denotes the rate by which the disease is transmitted from the vaccinated asymptomatic to the susceptible group upon contact,

● $\zeta^u$ denotes the rate by which the disease is transmitted from the unvaccinated symptomatic to the susceptible group upon contact,

● $\zeta^v$ denotes the rate by which the disease is transmitted from the vaccinated symptomatic to the susceptible group upon contact,

● $\epsilon^u$ is the rate of unvaccinated individuals becoming symptomatic infected after some latent incubation period,

● $\epsilon^v$ is the rate of vaccinated individuals becoming symptomatic infected after some latent incubation period,

● $\gamma^u$ is the recovery rate of unvaccinated infected individuals,

● $\gamma^v$ is the recovery rate of vaccinated infected individuals.

Importantly, the parameter values we considered above may now be region dependent, and so we designate them with a subscript, $\alpha_i: = [\theta_i, \beta_i^u, \beta_i^v, \zeta_i^u, \zeta_i^v, \epsilon_i^u, \epsilon_i^v, \gamma_i^u, \gamma_i^v]^T\in(\mathbb{R}_{\ge 0})^{9}$.

2.2.2.   Travel and SEIR Model's Discrete States

Consider an arbitrary region $i$. We include in our continuous dynamics the possibility that its mitigation policies, specifically, travel restrictions, allow travel to, and allow travelers from, any number of the other regions $j\ne i$. This potentiality couples the $i$th local (i.e., region specific) SEIR compartmental model with any number up to $n-1$ of the other $j\ne i$ region specific models. This is done via a directed graph whose nodes are the regions under consideration, and whose weighted edges correspond to proportions of travelers going to and coming from the various regions.

Definition 2.6. Let $i$ denote a node in the directed graph of regions. Using the typical vocabulary, the neighborhood of region $i$, denoted $Ngh(i)$, is the set of $j\ne i$ that can be reached by a single directed edge from $i$ to $j$. We make the simplifying assumption $j\in Ngh(i) \iff i\in Ngh(j)$. In other words, the decision to allow travel (or not) is always two-way, i.e., region $i$ and region $j$ both agree together to make the same type of mutual travel restrictions among themselves.

Note that, in general, $\text{N}gh(i)$ will vary with time, as the travel restrictions between two regions can change depending on the state of infection, a fact which is accounted for in the discrete dynamics.

Definition 2.7. In actuality, a neighborhood should exist for each compartment of the SEIR model, as individuals at various stages of infection will be allowed to travel differently. Hence, we have a set of eight neighborhoods $\{Ngh^{S^{u}}(i), Ngh^{S^{v}}(i), \dots, Ngh^{R^{v}}(i)\}$ of node $i$ for each of the eight compartments. These are then subject to the additional assumptions:

1. $Ngh^{I^u}(i) = Ngh^{I^v}(i) = \varnothing$

2. $Ngh^{R^u}(i) = Ngh^{S^u}(i)$ and $Ngh^{R^v}(i) = Ngh^{S^v}(i)$,

3. For each state variable $z_{i, s}\in\{S^{u}_i, S^{v}_i, \dots, R^{v}_i\}$,

where $\propto$ denotes "proportional to." Item 1) says that symptomatic individuals will never be allowed to travel. Item 2) says that recovered individuals are treated the same as susceptible individuals with regards to travel restriction policies. Item 3) indicates that we do not impose the symmetry $w_{ij} = w_{ji}$; rather, we balance the incoming and outgoing travelers in a different way, as described below in Eq (2.12).

Across all regions, we have $\binom{n}{2} = \frac{n(n-1)}{2}$ pairs $\{i, j\}$. Since our neighborhoods are symmetrically defined, and since we've effectively reduced (to four) the number of categories of neighborhood to be monitored, $\{Ngh^{S^{u}}(i), Ngh^{S^{v}}(i), Ngh^{E^{u}}(i), Ngh^{E^{v}}(i)\}$, we require $2n(n-1)$ neighborhoods to be defined in order to specify all of the travel policies. In total, these $2n(n-1)$ neighborhoods determine the discrete state $q\in Q$, as in Definition 2.1. These decisions will not be made in a vacuum or randomly; rather, they are based on strategies that model policies we saw implemented in the real world. In Section 2.2.2 below, we formulate $Q$ explicitly and calculate its cardinality.

In addition to the $2n(n-1)$ neighborhoods to be defined, there are $2\cdot\frac{n(n-1)}{2} = n(n-1)$ ordered pairs $(i, j)$ corresponding to the number of parameter weights $w_{ij}$ to be determined. These will be chosen so as to initially balance the number of incoming and outgoing travelers. Since only the infected individuals are not permitted to travel, we require:

More conditions need to be specified in order to uniquely determine all $n(n-1)$ weights. The above condition in Eq (2.12) reduces this number by half. For example, if $w_{ij}$ are chosen for $j > i$, then automatically all $w_{ji}$ would be uniquely determined. Indeed, this is the approach taken in subsequent discussions and in the simulations. What remains is specifying the $\frac{n(n-1)}{2}$ number of $w_{ij}$ for $j > i$. Two selection methods which seem natural are implemented. First is another type of balancing act, where outgoing travel numbers reflect the initial populations of the destinations, i.e., we assume that more people will travel to a destination which naturally supports a larger initial population ^[16]. In that regard, for a fixed region $i$, one may impose for all $j\ne k$ (and $j, k > i$) that

or, equivalently,

Second, we assume some known percentage of the regional population $N_i(0)$ is typically traveling at a time prior to any travel restrictions. For example, if $0.1\%$ of the population is typically traveling for work, vacation, or otherwise, we may impose $0.001 = \sum_{j\ne i}w_{ij}$. Together with Eqs (2.14) and (2.12), this final condition can be used to uniquely determine all values of the weights $w_{ij}$. One detail of note is, depending on the order in which these constraints are applied, it may happen that the last region with initial population $N_n(0)$ cannot enforce the final constraint of a known percentage of the population traveling at a given time, since if all $i\ne n$ have already done this and also the balancing Eqs (2.14) and 2.12 have been applied, then each $w_{ni}$ will already be determined and, therefore, no freedom of choice will remain to dictate the size of $\sum_{i = 1}^{n-1}w_{ni}$. In light of this, a final implementation strategy in the simulations will be to always index the regions from smallest initial population to largest, i.e., choose the region numbers such that $N_1(0)\le N_2(0)\le \dots\le N_{n-1}(0)\le N_n(0)$.

As witnessed during the COVID-19 pandemic, travel restrictions played an important role in containing the spread of the virus. They were shaped by policies that changed over time depending on the perceived severity of the pandemic. In the real world, travel restrictions are delicate to employ and cannot be modified too frequently for reasons of the practicality of their implementation, as well as the messaging, which must be delivered to the public. Indeed, they may change, but vary in a discrete way and not continuously. However, they do regulate some parts of the continuous dynamics, specifically with respect to the travelers, i.e., in particular, they determine the four neighborhoods $\text{N}gh^{S^u}(i)$, $\text{N}gh^{S^v}(i)$, $\text{N}gh^{E^u}(i)$, and $\text{N}gh^{E^v}(i)$. We describe the different travel restrictions that are possible in our scheme. They are based on what was actually implemented during the COVID-19 pandemic, but it is possible that others could also be considered.

Fix an arbitrary pair of regions $\{i, j\}$ and consider how they might arrange travel amongst themselves i.e. how might they choose their four neighborhoods. We notate the following $5$ strategies for selection:

1. open (O):

$j\in\text{N}gh^{S^u}(i)$, $j\in\text{N}gh^{E^u}(i)$, $j\in\text{N}gh^{S^v}(i)$, $j\in\text{N}gh^{E^v}(i)$.

An open region makes no restrictions on travel.

2. vax pass (P):

$j\notin\text{N}gh^{S^u}(i)$, $j\notin\text{N}gh^{E^u}(i)$, $j\in\text{N}gh^{S^v}(i)$, $j\in\text{N}gh^{E^v}(i)$.

A region implementing vaccine passports a.k.a. vax pass disallows unvaccinated travelers by requiring proof of vaccination upon arrival.

3. testing (T):

$j\in\text{N}gh^{S^u}(i)$, $j\notin\text{N}gh^{E^u}(i)$, $j\in\text{N}gh^{S^v}(i)$, $j\notin\text{N}gh^{E^v}(i)$.

A region implementing testing attempts to ensure no infected travelers arrive at all by testing whether they have been exposed and carry the virus without symptoms.

4. strict (PT):

$j\notin\text{N}gh^{S^u}(i)$, $j\notin\text{N}gh^{E^u}(i)$, $j\in\text{N}gh^{S^v}(i)$, $j\notin\text{N}gh^{E^v}(i)$.

A strict region implements both vaccine passports and testing.

5. closed (C):

$j\notin\text{N}gh^{S^u}(i)$, $j\notin\text{N}gh^{E^u}(i)$, $j\notin\text{N}gh^{S^v}(i)$, $j\notin\text{N}gh^{E^v}(i)$.

A closed region does not allow travelers at all.

So for the $\binom{n}{2}$ number of pairs $\{i, j\}$, their four neighborhood determinations about each other can come only from the five options $\{O, P, T, PT, C\}$ rather than the $16$ you might expect if we allowed the four decisions to be completely uncorrelated. Hence there are $5^{\binom{n}{2}} = 5^{\frac{n(n-1)}{2}} = (\sqrt{5})^{n^2-n}$ possible ways to form the four neighborhood designations to define the travel policies between all $n$ regions. Altogether, these neighborhood designations make up our hybrid automaton's set of discrete sets, $Q$.

2.2.3.   Hybrid dynamics, vector field $f$

The addition of travelers is the only change to the structure provided by Eqs (2.3)–(2.10):

Again, we assume $S_i^u(t)+E_i^u(t)+I_i^u(t)+R_i^u(t)+S_i^v(t)+E_i^v(t)+I_i^v(t)+R_i^v(t) = N_i(t)$ and furthermore we now have total population $N = \sum_{i = 1}^{n}N_i(t)$ for all $t\ge 0$.

The local dynamics and coupling in Equations (2.15)-(2.22) are of the general form seen in Equation 1.1, but the structure of $\mathcal{L}(x(t))$ is somewhat cumbersome to explicit due to the ordering of variables. Instead of $x_i = [S_i^u, E_i^u, I_i^u, R_i^u, S_i^v, E_i^v, I_i^v, R_i^v]^T\in\mathbb{R}_{\ge 0}^{8}$ and put altogether as $x = [x_1, \cdots, x_n]^T\in(\mathbb{R}_{\ge 0})^{8n}$, we use the following alternative: $S^u: = [S_1^u, \dots, S_n^u]^T, E^u: = [E_1^u, \dots, E_n^u]^T, I^u: = [I_1^u, \dots, I_n^u]^T,$ and $R^u: = [R_1^u, \dots, R_n^u]^T$ and similarly for the vaccinated individuals $S^v, E^v, I^v,$ and $R^v$. Although we already used previously in Eqs (2.3)–(2.10), the notation $S^u, E^u, \dots, I^v, R^v$, here they represent vector quantities instead. In context it should be clear that for one region they are scalar quantities and if there are $n > 1$ regions, they are vectors. We put all regions together as:

Similarly, we define vectors of parameter values for all regions $\theta: = [\theta_1, \dots, \theta_n]^T$, $\beta^u: = [\beta_1^u, \dots, \beta_n^u]^T$, $\zeta^u: = [\zeta_1^u, \dots, \zeta_n^u]^T$, $\epsilon^u: = [\epsilon_1^u, \dots, \epsilon_n^u]^T$, $\gamma^u: = [\gamma_1^u, \dots, \gamma_n^u]^T$, and likewise for $\beta^v, \zeta^v, \epsilon^v,$ and $\gamma^v$. Using $\circ$ notation for the Hadamard (also known as the Schur or entrywise) product, and suppressing time dependence on the righthand sides, we have succinctly

where the four matrices $L_{S^u}, L_{E^u}, L_{S^v}, L_{E^v}\in\mathbb{R}^{n\times n}$ are outgoing convention Laplacian matrices which, taken together, encode all of the travel policies afforded by the current discrete state $q$. The system can now be written as

with a matrix that is now conveniently block diagonal and $h$ is given by the local SEIR dynamics. The matrices $L_{S^u}, L_{E^u}, L_{S^v},$ and $L_{E^v}$ are not necessarily symmetric, as we have not imposed $w_{ij} = w_{ji}$.

To give an example of our four matrices and how they depend on $q$, suppose we have $n = 5$ regions, and each region allows travel to and from all four of the other regions, regardless of exposure or vaccination status i.e., each is in status open (O) with respect to all others. Furthermore for simplicity, assume we have weights $w_{ij} = w_{ji} = 0.05$ for all $i$ and $j$, then we would have

Now suppose that region $1$ reaches a severe level of infection, and so it may be the case that regions $2$ through $5$ want to enforce a vaccine passport (P) travel policy toward it. In that case, only vaccinated travelers would be allowed to travel, leaving neighborhoods $\text{N}gh^{S^v}(1)$ and $\text{N}gh^{E^v}(1)$ the same as before, but unvaccinated travelers would now be forbidden i.e., $\text{N}gh^{S^u}(1) = \text{N}gh^{E^u}(1) = \varnothing$. In that case $L_{S^v}$ and $L_{E^v}$ would stay the same as they are above, and instead we would have

Remember that the symmetry shown is not necessary in general, and only comes from our arbitrary example here that $w_{ij} = w_{ji} = 0.05$ for all $i$ and $j$. Altogether, the right sides of Eqs (2.24)–(2.31) are what defines our SEIR vector field $f(\cdot, \cdot):Q\times\mathbb{R}_{\ge 0}^{8n}\to\mathbb{R}_{\ge 0}^{8n}$ as in Definition 2.1. The initial states $\text{I}nit\subseteq Q\times\mathbb{R}_{\ge 0}^{8n}$ are simply the initial values at time zero, $y(0)$, along with the initial determination about travelers, $q_0\in Q$, presumably as in the first example above, with all regions open to all other regions (if we imagine the severe impact of the pandemic that has yet to occur at $t = 0$, then it is natural to suppose no travel restrictions have been made initially).

2.2.4.   SEIR model's domain and edges

The domain of our hybrid automaton $\text{D}om:Q\to 2^{\mathbb{R}_{\ge 0}^{8n}}$ is a map to all possible continuous states that are compatible with a discrete state $q$. In our case, some of the discrete states may be undesirable given a particular continuous state, but none are considered impossible. Therefore we allow all $N$ people to be in any of the $8$ compartments from any of the $n$ regions i.e. we have for all $q\in Q$,

To define a set of edges $E\subseteq Q\times Q$, we must consider when it makes sense to go from one discrete state $q$ to another.

The graph in Figure 3 encodes clearly the possible gradual changes we would expect for a fixed pair of regions $\{i, j\}$. Intuition tells us that changes to the policies on travel restrictions between regions $i$ and $j$ should be gradual, and should include at most two instances of turning "on" or "off" a type of travel eligibility. For a fixed pair $\{i, j\}$ consider the ordered quadruplet $W_{ij} = (W_{ij}^{S^u}, W_{ij}^{E^u}, W_{ij}^{S^v}, W_{ij}^{E^v})$ where we define

and similarly for the remaining travel groups $E^u, S^v$ and $E^v$. We only allow change from a policy at time $\tau_k'$, $W_{ij}(\tau_k')$, to a new policy at time $\tau_{k+1}$, $W_{ij}(\tau_{k+1})$ if the distance between the policies is at most two, $\|W_{ij}(\tau_k')-W_{ij}(\tau_{k+1})\|\leq 2$, where we take the norm $\|\cdot\|$ to be the taxicab metric i.e.

or given explicitly,

Along these same lines, define the matrix $W\in\{0, 1\}^{\binom{n}{2}\times 4}$ via

i.e. the $\frac{n(n-1)}{2}$ rows correspond to a choice from the set $\{O, P, T, PT, C\}$, which yields the four columns of zeros and ones. The matrix $W$ corresponds to a unique global travel strategy $q\in Q$ and encodes the zero and nonzero structure of the matrices $L_{S^u}, L_{E^u}, L_{S^v},$ and $L_{E^v}$. For example, consider again the scenario above when initially at times before some $\tau_0'$ the $n = 5$ regions were all open (O), but subsequently due to severe infection in region $1$, regions $2$ through $5$ switched to enforcing a vaccine passport (P) policy toward it at $\tau_{1}$. The corresponding matrices $W(\tau_0'), W(\tau_1)\in\{0, 1\}^{10\times 4}$ are given by

where one should note that the first four rows $W_{12}, W_{13}, W_{14},$ and $W_{15}$ changed from the open (O) policy $(1, 1, 1, 1)$ to the vaccine passport (P) policy $(0, 0, 1, 1)$.

An edge exists between two discrete states in $Q$ if their strategies are sufficiently similar i.e. if we consider going from one to the other to be a gradual change. To go from a discrete state at time $\tau_k'$, $W(\tau_k')$, to a new discrete state at time $\tau_{k+1}$, $W(\tau_{k+1})$, we require $\|W_{ij}(\tau_k')-W_{ij}(\tau_{k+1})\|\leq 2$ for all $\binom{n}{2}$ possible $\{i, j\}$ pairs. This is because it is not sufficient to impose the global condition $\|W(\tau_k')-W(\tau_{k+1})\|\leq 2\cdot\binom{n}{k} = n(n-1)$ since, for example, we could have $n = 3$ and think that the strategies $W(\tau_k') = [(1, 1, 1, 1);(1, 1, 1, 1);(1, 1, 1, 1)]$ and $W(\tau_{k+1}) = [(0, 0, 1, 0);(0, 0, 1, 0);(1, 1, 1, 1)]$ are sufficiently similar after calculating $\|W(\tau_k')-W(\tau_{k+1})\| = 6\leq 3(3-1)$. However, the physical meaning here is that region $1$ switched from an open (O) policy toward regions $2$ and $3$ to a strict (S) policy toward regions $2$ and $3$, without going through either vax pass (P) or testing (T) as we shall require, indicated clearly in Figure 3. Indeed, this counterexample is made evident by the fact that $\|W_{23}(\tau_k')-W_{23}(\tau_{k+1})\| = 0\leq 2$, but on the other hand, $\|W_{12}(\tau_k')-W_{12}(\tau_{k+1})\| = \|W_{13}(\tau_k')-W_{13}(\tau_{k+1})\| = 3 > 2$, which we will not allow.

2.2.5.   SEIR model's guard condition and reset map

We require the guard condition $\text{G}rd(\cdot):E\to 2^{\mathbb{R}_{\ge 0}^{8n}}$ to encode when we can potentially switch between discrete states. Suppose the discrete state is at $q_0\in Q$. If at some point in time, $t$, the continuous state $y(t)$ reaches the guard $\text{G}rd(q_0, q)\subseteq\mathbb{R}_{\ge 0}^{8n}$ for some edge $(q_0, q)\in E$, then the discrete state $q_0$ may change to the discrete state $q$. In the general theory of hybrid systems, we could simultaneously allow for the continuous state $y(t)$ to be reset to some value in $\text{R}es(q_0, q, y)\subseteq\mathbb{R}_{\ge 0}^{8n}$; however, in this particular application we do not need a reset, so we take the reset map to be the identity map, i.e., $\text{R}es(\cdot, y) = y$ regardless of the edge input $E$. We also wish to remove any ambiguity in the guard condition. Therefore, we construct $\text{G}rd(\cdot):E\to 2^{\mathbb{R}_{\ge 0}^{8n}}$ below such that one and only one transition is made eligible, and, furthermore, we assume that the transition is always taken.

The guard condition will take both the relative number of infected in region $i$ and the relative number of infected in another region $j\ne i$ into account. Moreover, it is reasonable to assume these decisions must be independent of the sign of the difference between these relative values. The simple option we implement here is to make travel policy decisions based on

Intuitively, we look at the numbers of symptomatic infected $I$ and not the asymptomatic infected $E$, since they would be even less likely to show up in any type of publicly reported data. Furthermore, we look only at the relative level of infected. Since a larger region could have a much larger population of infected without cause for concern, this value should be normalized. Consider thresholds that we might potentially reach for these values, $0 < \phi_1^u < \phi_2^u < 1$ and $0 < \phi_1^v < \phi_2^v < 1$, respectively. Suppose we say the values above are "good" if they fall within $[0, \phi_1^{(*)}]$, "bad" if they fall within $[\phi_2^{(*)}, 1]$, and "okay" if they are somewhere in the middle, $(\phi_1^{(*)}, \phi_2^{(*)})$. Recall an edge exists in $E$ if $\|W_{ij}(\tau_k')-W_{ij}(\tau_{k+1})\|\leq 2$ for all $\binom{n}{2}$ possible $\{i, j\}$ pairs. Assume that this holds, then for an $\{i, j\}$ pair, if the values $D_{ij}^u$ and $D_{ij}^v$ are both good, taking our new strategy to be open, $W_{ij}(\tau_{k+1}) = O$ makes sense. On the other hand, if the values are both bad, then taking a closed strategy, $W_{ij}(\tau_{k+1}) = C$ makes sense. If both values are okay, then we take the strict policy $W_{ij}(\tau_{k+1}) = PT$. If the value $D_{ij}^u$ is okay, but the value of $D_{ij}^v$ is good, then we use vaccine passports, $W_{ij}(\tau_{k+1}) = P$. Lastly, if the value of $D_{ij}^u$ is good, but the value of $D_{ij}^v$ is okay, then we implement testing, $W_{ij}(\tau_{k+1}) = T$.

We illustrate these inequalities in Figure 4 below, assuming that we have thresholds $\phi_1^u = 0.1$, $\phi_2^u = 0.2$ and double for vaccinated individuals, $\phi_1^v = 0.2$ and $\phi_2^v = 0.4$. Topologically, the green box is closed and corresponds to the open (O) policy. The red region is also topologically closed and corresponds to the closed (C) policy. The orange box is topologically open and corresponds to the strict (PT) policy. The yellow box above the green is neither topologically open nor closed, as it contains its sides, but not its top or bottom, and corresponds to the testing (T) policy. Similarly the yellow box to the right is neither open nor closed, contains its top and bottom but not its sides, and corresponds to the vaccine passport (P) policy. In the simulations these values are taken to be much smaller, but here are chosen to clearly view the figure.

3.   Simulations of disease spread utilizing the hybrid automaton

In this section we include a set of simulations to illustrate the functionality of our hybrid automaton model. The goal is not to recreate an especially realistic real world scenario, nor to fit any existing data, rather it is a proof of concept to demonstrate the feasibility of our model for computing at a scale of the global population, more than 8 billion people as of November 15, 2022 ^[17]. Furthermore, we highlight the automated travel restriction strategies as determined by the guard conditions. The model has been coded in MATLAB R2022b. When timing the code, we repeat the procedure some number of times (specified case by case below) and take the average. This is done in an attempt to reach a more accurate value, since variance or random noise can cause fluctuations in the time required to finish a single run. This is particularly helpful when running simulations on limited and/or outdated hardware, e.g., the first author's 2015 MacBook Pro.

3.1.   Computation time depending on population and number of regions

Here we analyze the computational time required to numerically integrate Eqs (2.24)–(2.31) as we vary the total population $N$ and the number of regions $n$. The difficulty in implementing naively the formalism of a hybrid automaton comes from the fact that the guard conditions potentially introduce discontinuities in the continuous dynamics. Although there is no reset map involved in our particular application here, the continuous state variables themselves do not ever change discontinuously, and their derivatives do experience discontinuities in the form of the discrete changes happening to the matrices $L_{S^u}, L_{E^u}, L_{S^v},$ and $L_{E^v}$, as they are impacted by the travel restrictions, triggered by the guard conditions that depend on the values of $D_{ij}^u(t)$ and $D_{ij}^v(t)$ defined in Eq (2.36). Hence the numerical integration must be stopped and then restarted each time a guard condition is reached, which naturally slows down the speed of the computation. The values for the corresponding thresholds are taken to be $(\phi_1^u, \phi_2^u) = (0.01, 0.02)$ for unvaccinated individuals and $(\phi_1^v, \phi_2^v) = (0.025, 0.05)$ for vaccinated. This is due to the fact that even a lesser level of infection is considered more alarming for unvaccinated individuals than for vaccinated individuals, since we assume the vaccine has the effect of minimizing the severe impact of symptomatic infection on an individual's health.

For all simulations we use a step size of $0.01$, which simplifies the output for plotting purposes. Note that by fixing the step size, a guard detection may be counted twice with precision defined between two subsequent steps. The solution we took was to throw out two guard conditions of the same type that were met in less time than the step size squared, $0.01^2 = 0.0001$. Adaptive step size numerical solvers should, in principle, largely avoid this issue.

3.1.1.   Varying the global population $N$

To begin, we test the computation time required to numerically integrate the single region SEIR model given by Eqs (2.3)–(2.10) over a three year timeline i.e., from $t = 0$ to $t = 1095$. This duration was chosen because without mitigation measures or a more sophisticated model that includes the possibility of reinfection, the most critical dynamics between interacting compartments will play out relatively quickly given our choice of parameters. Already after three years, solutions tend toward a steady state consisting of only disease free individuals, as seen in Figure 6. In this simulation, the discussion above in Section 3.1 about guard conditions does not yet apply, as these are not yet incorporated into the model since we have no travelers between different regions. The chosen parameter values are displayed in Table 2.

Here we averaged the computation time required over $15$ trials. The total population size is varied from $N = 1,000,000$ to $N = 8,041,000,000$, which is approximately the current world population. We increase the population in each iteration by $40,000,000$, for a total of $201$ data points, displayed in Figure 5. The initial population distribution is chosen in each case to start without any vaccinated individuals, but we add vaccinated individuals throughout via the parameter $\theta$. We start out with $99\%$ of the population susceptible, $\frac{1}{2}\%$ of the population exposed, and $\frac{1}{2}\%$ of the population infected, i.e. we have initial values shown in Table 1. The outcome of the final simulation's eight compartments for $N = 8,041,000,000$ is also shown in Figure 6. We conclude that the computational cost due to varying the total population $N$ is negligible in the single region case, as the best linear fit shown in Figure 5 only indicates, on average, an increase of roughly $0.05$ seconds to vary from a population of $N = 1,000,000$ to a population of $N = 8,041,000,000$.

Next we we test the computation time required to numerically integrate the hybrid model given by Eqs (2.24)–(2.31) for a world with $n = 6$ regions. Once again, we vary the world population from $N = 1,000,000$ to $N = 8,041,000,000$, increasing the population by $40,000,000$ each time, and we integrate over a three year timeline. Since the hybrid model incorporates the automated travel restrictions, now the guard conditions do come into play, and we must repeatedly stop and restart the numerical integration process each time a guard condition is satisfied. Hence we anticipate a longer computation time as compared to the single region $n = 1$ scenario above. We again average the computation time over $15$ trials. The initial population distributions are chosen in the same way as shown in Table 1. Notice now however, the subscript $i$ referring to the region number, and the explicit time dependence in the initial regional population $N_i(0)$, since $N_i(t)$ will in general vary in this case. We split up the total population $N$ into the six initial regional populations:

We have $N = \sum_{i = 1}^{6}N_i(0)$ as desired. The six values in Eq (3.1) were generated in the following way. First, using the MATLAB function normrnd, six random numbers were drawn from the normal distribution with mean zero and standard deviation one. Second, we take the absolute value of the random numbers and sum them. Lastly, we divide each by the sum so that the new sum is one as desired. We use this strategy again subsequently when we vary the number of regions $n$ while leaving the total population $N$ fixed. The parameter values are again those chosen in Table 2, except that the multipliers on the $\beta^{(*)}$ and $\zeta^{(*)}$ parameters are normalized by region i.e. we scale them each by $\frac{B}{N_i(0)}$ instead of $\frac{B}{N}$. The results of timing the computation are shown in Figure 7. Again, we conclude that the computational cost due to varying the total population $N$ is negligible in the hybrid model, as the best linear fit shown in Figure 7 only indicates, on average, an increase of less than $0.05$ seconds to vary $N = 1,000,000$ to a population of $N = 8,041,000,000$. However, it is notable that the average speed as compared to the single region model has increased by about $0.3$ seconds, a seemingly small change, but one that will become dramatically more pronounced as the number of regions $n$ is increased.

3.1.2.   Varying the number of regions $n$

Now we test the computation time required to numerically integrate the hybrid model given by Eqs (2.24)–(2.31) for a world with $N = 8,000,000,000$ people (just short of the current global population) as we vary the number of regions from $n = 3$ to $n = 303$ (the United Nations recognizes 195 countries worldwide). As indicated in the previous section, the dependence of the equations on $n$ contributes greatly to the computation time required. To decrease the computation time we integrated over two years instead of three, and similarly, we average the computation time over $5$ trials. Yet again, we take initial population distributions in the same way, with $99\%$ of the regional populations susceptible, $\frac{1}{2}\%$ exposed, and $\frac{1}{2}\%$ infected, and we use the parameter values in Table 2, except again the multipliers on the $\beta^{(*)}$ and $\zeta^{(*)}$ parameters are normalized by region i.e., we scale them each by $\frac{B}{N_i(0)}$ instead of $\frac{B}{N}$. The results of timing the computation are shown in Figure 8. We conclude that a global model is feasible with the simplified set of equations we are using. Using 195 regions corresponding to countries recognized by the United Nations will take on average 7.4 seconds per simulation. This is fast compared to, for instance, an agent-based model that will take hours to simulate the state of California alone. There are thousands of international airports, however, and so we hope to run the computation further in future work.

3.2.   Simulations to the scale of the world and role of guard conditions

As of April 27, 2023, the approximate populations of Oceania, North America, South America (includes here Latin American countries and the Caribbean), Europe, Africa, and Asia are respectively given in Table 3. We label these from $1$ to $6$ with a total of $N = \sum_{i = 1}^{6}N_i(0) = 8,031,247,037$ people in our model at the scale of the world. We would like to see the impact of the automated travel restrictions in simulations of hypothetical situations taking place over seven years. What if the largest region, Asia, is initially five times as infected (relatively) as the other five regions? How does this scenario compare to the same scenario where no travel restrictions are imposed i.e., all regions remain open (O) throughout the seven years? We call these scenarios $1$ and $2$, see Table 4.

We take initially the initial conditions found in Table 5.

In Asia $0.95\%$ exposed and $0.95\%$ infected are just below the critical $\phi_1^u = 0.01$ threshold. The parameters are again those in Table 2. Importantly, one final parameter is the assumed amount of travel at a given moment. We take from some online estimates six million people flying at a given time, out of a total population of roughly eight billion people yields $\frac{6\times 10^6}{8\times 10^9} = 0.00075$ i.e., we shall assume $0.00075 = \sum_{j\ne i}w_{ij}$ for $i = 1$ through $i = 5$. However as noted in Section 2.2.1, in order to balance incoming and outgoing travelers via Eqs (2.12)–(2.14), we have that $0.0003598 = \sum_{j = 1}^{5}w_{6j}$ in Asia. The $w_{ij}$ are given below in a matrix whose row number and column number correspond to the indices $i$ and $j$, respectively. For display purposes only the values here have been rounded to show four significant digits, however MATLAB itself utilizes $16$ digits of precision.

3.2.1.   Scenario 1–travel restrictions

There are several results worth mentioning about the scenario $1$ with travel restrictions. The entire world travel is closed (C) by day $326$ (a bit less than one year). By day $1721$ (4 years and 261 days), all regions are open (O) to each other except everyone is still closed to Asia. By day $1853$ (5 years and 28 days), full worldwide travel is restored to all regions. The full details of all of the changes in policy are given in Table 6.

Interestingly, the pair of smaller regions $1-$Oceania and $2-$North America had to change the policy between themselves at least as much as any other pair of regions in Scenario $1$ above. The result is shown below in Figure 9. One can see the color changes from forest green, to periwinkle, to brown relatively quickly, reflecting changes in policy from strict (PT), to vax pass (P), to open (O). These two changes occurred on days $1662$ and $1719$, nearly a two month span which is relatively realistic in the sense that there is no lower bound on how fast the changes can happen, a shortcoming of the model.

3.2.2.   Comparison of disease evolution between Scenario 1 and Scenario 2

Below we refer to subsequent figures and summarize our other findings. Any plot quantity which has been 'normalized' has been divided at each data point by its associated region's instantaneous total population.

1. Susceptible vs. Exposed: In Figure 10, we see the plots look nearly identical, except that the difference between Asia and the other five regions is much more exaggerated with travel restrictions. This is because with restrictions the populations are not allowed to mix as much, and so the virus is somewhat kept in Asia in that scenario.

2. Susceptible vs. Infected: Likewise, comparing plots below in Figure 11 we find Asia actually had more infected in Scenario 1 with travel restrictions.

3. Exposed vs. Infected: Comparing plots in Figure 12, staying open (Scenario 2) stratifies the curves of the $5$ smaller regions and shows noticeably more infected e.g., Africa's orange curve falls short of $0.1$ relative infected with travel restrictions, but exceeds this in the open scenario. Note that $1\%$ of Africa is more than $14$ million people, and therefore represents a significant number of individuals. Figures 13 and 14 further make clear these findings.

3.2.3.   Comparing Number of Exposed and Infected Individuals between Scenario 1 and Scenario 2

In Figure 13 we see the peak number of total normalized infected in the world is noticeably less in the travel restriction curve shown in blue as compared to the orange curve without travel restrictions. The small gap corresponds to about $16,037,175$ additional people infected at the peak if no travel restrictions are imposed. Another benefit is that the peak is delayed by nearly nine days. However, in Figure 14 we actually see the opposite. Asia experienced a greater peak by implementing travel restrictions. The gap corresponds to about $7,526,067$ additional people infected by not just staying open. The reason is that since we assume an initial balance of incoming and outgoing travelers in the model, a mixing with populations that are less infected would be beneficial. There is not a major payoff in terms of time difference, as the greater peak with travel restrictions occurs only a week earlier. Figure 14 shows a decidedly greater peak in relative infected in Africa by staying open and not implementing restrictions, corresponding to $6,150,493$ people. Indeed if some percentage of infected get severely ill and require hospitalization, just $1\%$ of this difference is more than sixty thousand people. Quite importantly, there is also nearly a $37$ day delay in when the peak occurs. The extra month could be incredibly important in preparing additional medical facilities or waiting for supplies to arrive.

All of the other four smaller regions not shown have the same exact relative curves as Africa in Figure 14, since they had the same homogeneous populations and disease spread parameters. Therefore the globally optimal strategy is to close off travel to and from Asia just short of $11$ months, even though this is worse for them in the short term. However, notice that both sets of curves flip eventually. In Asia, by the end of seven years, the open strategy ends up being on top with more infected. The opposite is true in Africa with more infected on the first day of the seventh year by having implemented the travel restrictions. However we know that flattening the curve was the general goal of most regions that implemented travel restrictions during the COVID-19 pandemic. In conclusion, the travel restrictions were a successful mitigation strategy in this scenario. They flattened the curve in the five smaller regions and also worldwide. They delayed the peak infection in these regions by more than a month.

This conclusion is further supported by the values displayed in Table 7. Consider that any of the $R_i(t)$ compartments are always strictly increasing in size, with contributions from the individuals in the $I_i(t)$ compartments. Therefore the total infected in each region after seven years is the final value in the recovered category, $R_i(2555)$ for $7\times 365 = 2555$ days. The outcome of the preferred scenario in each row is in green, due to less total infections, and the outcome with more infections is in red. Only Asia sees a benefit in staying open without travel restrictions.

3.2.4.   Comparing estimate impact on heath care system between Scenario 1 and Scenario 2

A metric that was used throughout the pandemic to assess the gravity of the situation at any moment was the count of active COVID-19 hospitalizations. A primary concern was to make sure the health system would not be overwhelmed, especially in the ICU with ventilators. A natural assumption to make is that the number of active hospitalizations in each region is proportional to the number of active infections. Assume for now only the special case where this percentage hospitalized is the same for each of the six regions. Simply comparing normalized values of total active infections, $\frac{I_i(t)}{N_i(t)}$ and $\frac{I_j(t)}{N_j(t)}$ for $i\ne j$, is a useful way of understanding the impact on different regions. Moreover, we wish to simultaneously consider the impact between the two scenarios. Specifically, we consider the curves

where the superscripts $(1)$ and $(2)$ refer to either Scenario 1 or Scenario 2, respectively. In this way, we can see that if $\phi_i(t) > 0$ then no travel restrictions would have been preferred since this means Scenario 1 is giving more infections, but if $\phi_i(t) < 0$, then travel restrictions will provide more protection. The $\phi_i(t)$ curves for $i = 6$ (Asia) and $i = 5$ (Africa) are shown in Figure 15. One can observe that without travel restrictions (Scenario 2), there is less strain on hospital systems early on as one would expect; however, post infection-peak, the strain on hospitals is higher. For Africa we see the opposite phenomenon. The most important feature to observe is difference in the peaks of these two curves. Asia has a positive value of about 2.5 and Africa a negative value of about -10. This implies that travel restrictions were most effective at mitigating the stress put on hospitals in Africa by 30% compared to Asia. It is a non-negligible amount, especially considering Africa's limited hospital resources.

4.   Conclusion

We demonstrated that a simplified model of disease spread can still be a useful tool when supplemented with the hybrid system, as we were able to run simulations with a global population. Furthermore, we were able to automate the travel restrictions in a way that allowed us to draw conclusions about how regions working together could help each other and lessen the impact of disease. Moreover, we demonstrated that a global model at the scale of the world's population was computationally feasible, even when incorporating the notion of automated travel restrictions. Making such a model that incorporates available data from the real world would be very challenging since during the COVID-19 pandemic, there was not a global protocol in place for coordinating travel restrictions.

5.   Discussion

Applying the guard conditions to control travel between regions and therefore the evolution of the disease demonstrated that it might be more efficient globally if some regions make decisions that are detrimental locally while keeping these regions' local spread within a reasonable threshold. This observation suggests that a global optimal control problem with constraints on the local dynamics could potentially provide a much better outcome. It is, however, much harder to formulate and then analyze. Regarding the guard conditions, note that the $(D^u_{ij}, D^v_{ij})$ curves are interesting to study in general. They are always continuous, but do not need to be differentiable, i.e., they can have kinks and corners, or even be self intersecting, for example, Figure 16 generated for a simulation with different parameters and initial populations. The difficulty in numerically implementing the guard condition methodology described in Section 2.2.5 cannot be overstated. All of the events that the MATLAB ODE Event Location option take in must be described in terms of a simple equality, for example $D^u_{12}-\phi^u_1 = 0$. Fortunately the detection through zero can be made in either direction i.e., crossing from positive to negative or vice versa. This helps simplify the implementation somewhat, but it is still not enough to completely specify the guard conditions. The curve passing through different colored regions in Figure 9 and Figure 16 corresponds to actual changes in travel policy that should be made, whereas the color changes in the curve itself are just numerical instances of an equality being met and the integration stopping. Therefore these must be further probed for whether or not they warrant an actual change in travel policy. Indeed, one can see in Figure 16 that every time the background color changes, so does the color of the curve, but not the other way around. Precisely, a color change to the curve occurs if any of $D^u_{12}-\phi^u_1, D^u_{12}-\phi^u_2, D^v_{12}-\phi^v_1,$ or $D^v_{12}-\phi^v_2$ equals zero. This can happen with no change of policy needed, for example, green to black, black to pink, and pink to yellow in Figure 9.

There are multiples ways to improve upon our epidemiological model. First, the local compartmental models could be made more complex, for instance by distinguishing types of individuals at different levels of risk like health care workers, students, etc. However, it is not clear if these alterations would have a significant impact on the global travel restriction strategies. Second, we could imagine travel restrictions not to be symmetrical between two regions i.e., using a graph Laplacian whose symmetry is broken not just by the weights, but by the existence of nonzero and zero entries. This was actually the case at some points during the pandemic, for instance, between the United States and Japan. Third, for the guard conditions, vaccination rates themselves could be used in addition to relative levels of infection. Indeed, a lower vaccination rate should be associated with a lower threshold of infections for travel restrictions, since the disease will presumably spread faster in that region. The model could also impose that guard conditions cannot occur unless a minimal time duration has elapsed, because changing strategies too often confuses the population and could be detrimental to compliance. In practice this implementation is straightforward in simulations, but it is not naturally described with the most basic hybrid automaton structure.

Use of AI tools declaration

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

Acknowledgments

This material is based upon work supported in part by the National Science Foundation under Grant CCF-2030789.

Conflict of interest

All authors declare no conflicts of interest in this paper.