A stochastic agent-based model to evaluate COVID-19 transmission influenced by human mobility

To investigate the impact of human mobility on the transmission of COVID-19, we develop a stochastic agent-based Mob-Cov model. The model includes three key components: (1) a synthetic population generated with the hierarchical geographical characteristics (Figure 1a); (2) realistic human travel rules, taking into account the power law behaviors of travelling within a container and transporting between containers (Figure 1b); (3) the infection and recovery with probabilities influenced by individual health conditions (Figure 1c). The first two components about human mobility are discussed in subsections II B and II C. The last component about disease infection and recovery is discussed in the subsection II D. The flow chart (framework) of the Mob-Cov model is shown in Figure 2.

This model framework includes a flexible individual-based approach that can capture not only the local interactions between individuals and microscopic disease transmission, but also the global spread of the disease influenced by the really long-distance cross-district travels. The framework allows policymakers to define measures at various spatial-temporal scales and evaluate the effective length scale of the interventions.

In the following sections, the three key components of the model and their mathematical foundations are first introduced. The selection of model parameters and the evaluation methods are then presented. Lastly, the pseudo code of the model and the computational procedures are given.

The hidden geographical structure of the population is modeled using nested containers proposed by Alessandretti et al.24 . Specifically, the physical space is assumed to contain hierarchy of $L$ levels. Level 1 is the smallest and level $L$ is the largest. At any level $l$ ($l=1,2,...,L$), the space is partitioned into $n_{l}$ topologically compact containers. For $l<L$, a container is fully included within a single parent container. Thus, each geographical location $k$ can be uniquely identified as a sequence of containers ($k=(k_{1},k_{2},...,k_{L})$), where the child container at $l$ level $k_{l}$ is fully included in the mother container at $l+1$ level $k_{l+1}$. For example, Figure 1a shows the schematics of a three-level hierarchical structure. At the largest level (level 3), there are 3 containers, as shown by the light orange polygons. Level 2 has 9 containers, as shown by the orange polygons. And level 3 contains 21 individual positions, represented by the black dots. The location of agent $j$ is defined as $j=(j_{1},j_{2},j_{3})$, where $j_{3}$ is the ID of the left large polygons in level 3, $j_{2}$ represents the middle polygon in $j_{3}$ and $j_{1}$ is the ID of one black dot in $j_{2}$.

Travelling in the nested containers follows a two-stage process: first deciding which level to travel at, second selecting a container at that level, and then a container within the chosen container and so on down to the lowest level. Specifically, suppose an agent is at position $j=(j_{1},j_{2},...,j_{L})$, the probability for the individual to travel to position $k=(k_{1},k_{2},...,k_{L})$ can be expressed as (Figure 1b)

where $P_{lv}(l)$ is the probability of the agent choosing level $l$. The agent stays in the same container if the container level is larger than $l$. In $1,2,...,l$ levels, the agent switches to different containers. Thus, $k_{l+1}=j_{l+1}$, $k_{l+2}=j_{l+2}$, …, $k_{L}=j_{L}$, but $k_{0}\neq j_{0}$, $k_{1}\neq j_{1}$, …, $k_{l}\neq j_{l}$. $P_{lv}(l)$ follows an exponential distribution with a parameter $d$, which is expressed as

where $l$ is an integer representing the level. Eq. (2) indicates that switching between large-level containers (e.g. between countries) is less frequent than switching between small-level containers (e.g. between houses or streets).

$P_{i}(k_{i})$ is the probability of choose container $k_{i}$ at level $i$. $P_{i}(k_{i})$ follows a long-tail distribution with parameters $c_{0}$ and $c_{1}$, which is expressed as

where $i$ is the level. $k$ is the container position. Eq. (3) suggests that the travel probability decreases over the distance, especially at a large level $i$. This indicates that, for example, traveling to a distant country is more scarcer than going to a remote road. But compared with the widely-used exponential distribution, the long-tail distribution allows the happening of long-distance movements at all levels.

Eqs. (1)-(3) calculate the probability of switching from location $j$ to location $k$ if the agent takes a movement. We consider the fact that people do not travel at every time moment but need to take a rest or stay at a place for a while. Thus, the probability for an agent to take a movement is modeled with a probability $P_{m}$, which is uniformly distributed between 0.2 and 0.4. A portion of people with low movement tendency (i.e. low mobility population) is further considered. The ratio of the low mobility people in the population $r_{lm}$ is a model parameter. The probability for the low mobility people to take a movement is only $P_{m}=0.03$.

If two agents $i$ and $j$ ($i$ is infected and $j$ is healthy) are in the same level-2 container (level-1 containers are the individual positions), infection happens with a probability $P_{ij}$ (Figure 1c). If there are more than one infected agents, the infection probability of a healthy agent $j$ is calculated as $P_{j}=1-\prod_{i}(1-P_{ij})$, for all infected agents $i$ in the same level-2 container with $j$. For the high risk population, $P_{ij}$ follows a uniform distribution between 0.1 and 0.3. For the low risk population, $P_{ij}=0.02$. The low risk people ratio in the population $r_{li}$ is a model parameter. Note that, $P_{ij}$ is chosen according to previous studies 15 ; 16 . Kerr et al. 15 set $P_{ij}=0.05$ for household-level transmission. Cuevas 16 set $P_{ij}=0.2$.

After being infected for a time period $T_{inf}$, the individual $i$ is able to recover to a healthy state with a probability $R_{i}$ (Figure 1c).

To evaluate the spread speed of the disease in the population, the outbreak endpoint is calculated. Using the outbreak endpoint instead of the time when all the agents are infected is because there are always few individuals having really low infectious probability $P_{ij}$. The time it takes to get the whole population infected is heavily biased by the last few individuals, which cannot directly and accurately reflect the system property. Therefore, an outbreak endpoint which considers the effects of both the infected population ratio and the infection rate is used in this study13 ; 30 .

The detection of the outbreak endpoint is similar to find the knee point in system engineering. The knee point usually represents the "right decision point", after which the system properties are no longer significantly influenced by the variation of the variable30 . Thus, in this study, we follow the previously proposed methods and calculate the maximum point of the cost function. The maximum point is defined as the outbreak endpoint of the system. The cost function is expressed as13

where $N_{inf}$ is the number of infected agents at time step (iteration step) $t$, $N_{p}$ is the total population number, $T_{max}$ is the total time (maximum iteration number), which is defined as the time when 90% of the population got infected plus 100 extra time steps.

It has been proved that the cost function $J$ can only have one global maximum point, which corresponds to the outbreak endpoint13 . Therefore, the outbreak time can be calculated as

Figure 3 shows a simulation example of the infection rate over time and the determination of the outbreak endpoint. Figure 3a shows that the number of newly infected people is the largest during the first 70 time steps. After 70 time steps, the infection rates gradually decrease to values close to zero. As expected, the proportion of infected people keeps increasing over time, but the increasing speed slows down (Figure 3b). By finding the maximum point of the cost function (Figure 3c), the outbreak point can be identified when $t$ is around 70, which corresponds to an infection ratio around 80%. The results suggest that the first 70 time steps are the most critical period when the number of infected people increases rapidly. After around 80% of the population gets infected, the infection dynamics slows down. And it can take a really long time (can be longer than 2000 time steps) for the whole population to get infected.

It should also be noted that the outbreak happens when only infection of individuals is included in the model. After considering the recovery of sick people, the system is likely to reach a dynamical equilibrium with both infected and healthy population, instead of reaching the outbreak endpoint.

The Mob-Cov model is implemented as an iterative scheme. The flowchart is shown in Figure 2 and the pseudo-code is shown below. In the pseudo-code, $L$ is the number of levels. $N_{cont}=(n_{1},n_{2},...,n_{L})$ is the number of containers in each level. $c_{0}$, $c_{1}$ and $d$ are parameters that determine the mobility, as discussed in section II2. $N_{p}$ is the number of agents. $N_{inf0}$ is the number of initially infected agents. $T_{inf}$ is the infection time period before the agent having a probability to recovery. $r_{lm}$ is the ratio of low mobility population. $r_{li}$ is the ratio of low infection risk population. NestedCont is the nested container structure, which is a nested random dictionary implemented in Python. $P_{mi}$ is the probability of taking a movement for agent $i$. $P_{ij}$ is the probability of agent $j$ getting infected by agent $i$. $R_{i}$ is the recovery rate of agent $i$. If agents $i$ and $j$ are in the same level-2 container, $j$ is defined as the neighboring agent of $i$.

In the simulation, the maximum iteration number is set as 2000, which is long enough for the system to reach the equilibrium.The number of agents is 1000. The number of initially infected people is 5. For each experiment, 50 simulations are performed, and the mean and variance of the results (including the outbreak time, infected population ratio when outbreak, infected ratio, zero-COVID ratio) are calculated.

Figure 4a shows a simulation result of the population distribution and the spread of COVID-19 disease. In the simulation, the number of levels is $L=5$. Parameters of mobility are $d=1.5$, $c_{0}=0.8$, $c_{1}=0.35$. The number of containers in each level is $N_{cont}=(8,4,9,9,250)$. Containers in level 5 are represented by colorful and transparent polygons. Individuals in different level-5 containers are represented by dots labeled with different colors. Figure 4b is the probability of selecting a container level, $P_{lv}(l)$. The grey bars are simulation results, which fit well with the analytical expression in Eq. (2). Figure 4c is the probability of selecting a container in level 2, $p_{2}(k)$. The simulation results represented by the black dots fit well to the analytical expression in Eq. (3).

At the beginning, only five individuals (four in the blue container and one in the red container) are infected. Due to the high population density in the blue container, the disease quickly spreads out in one cluster at time $t_{0}$. Many infected people then travels to the nearby cluster in the same blue container, accelerating the local transmission. Meanwhile, the spread of the disease is relatively slow in the red container due to the low population density. Afterwards, one or two infected people travel between the level-5 containers and reach the grey and brown areas, indicating that the disease is transmitted globally. After a further local spread at time $t_{2}$, a disease outbreak happens at $t_{3}$.

Overall, the spread of the disease is a multi-scale dynamic process influenced by the probability of switching between levels (Figure 4b) and the movement distance (Figure 4c). The results indicate that the total population density and the initial infected people number does not need to be really large. When the disease starts from a local population cluster with a high people density (e.g. a crowded market), a small probability of traveling between large-level containers (e.g. traveling cross cities or countries) can lead to a rapid global spread and outbreak. The local-to-global spread of the disease heavily relies on the long movements between containers. Thus, accurately modeling both the short walks and long sojourns using the hierarchical container structure in the Mob-Cov model ensures the precise prediction of the disease dynamics.

The spread and outbreak of COVID-19 is influenced by various complex factors, such as biological characteristics of the virus, health and immune conditions, population density, environmental conditions and policies. Using the Mob-Cov model, we mainly decipher the impact of microscopic and universal human mobility and individual health conditions on the disease spread and outbreak (Figure 5).

For the mobility, the results show that the moving distance and the probability of traveling between large-level containers have significant influence on the disease spread and the onset of the outbreak. When long movements in the same level happen more frequently, which corresponds to a smaller value of $c_{1}$, the local crowdedness is reduced. Thus, it takes a longer time to reach outbreak (Figure 5b). The infected people are likely to be distributed in every containers and the ratio of the infected population can be relatively small when the outbreak happens (Figure 5a). When individuals are likely to travel between large-level containers (e.g. cities and countries), which corresponds to a small value of $d$, the global transmission of the disease is enhanced. The time it takes for the onset of the outbreak becomes shorter. At the outbreak endpoint, more than 80% of the population have been infected. Restricting travels between large-level containers helps slow down the spread of the disease and reduce the size of the infected population (Figure 5c and Figure 5d). The ratio of the low mobility population over the total population does not have a significant influence on the infected people number when the outbreak happens (Figure 5e). But if most people prefer to stay at the same place than move, the time it takes for reaching the outbreak endpoint increases from 100 time steps to more than 600 time steps. In addition, the variation of the outbreak time is large when the mobility of most people is low, because the outbreak time now heavily relies on the random initial distribution of the population (Figure 5f).

For the individual health conditions, if a large percentage of the population have low infection risk (e.g. most people get a COVID-19 vaccination), the spread speed of the disease is slowed down (Figure 5h). The low risk people ratio does not have significant influence on the infected population size at the outbreak endpoint. Around 70% of the people are infected at the outbreak (Figure 5g). The total population size also does not have a significant impact on the infected population ratio (Figure 5i). But when the population number is large, the local crowdedness increases the infection probability of individuals and accelerates the spread of the disease (Figure 5i).

Lastly, we find that the hierarchical structure of the containers has a large impact on the infected population ratio and the outbreak time. If the hidden network of human travels contains many levels, the spread speed of the disease is slower and the outbreak time is larger (Figure 5l). In addition, less people get infected at the outbreak endpoint (Figure 5k). The result suggests that multi-layer geographical barriers and borders that restrict human mobility can help reduce the spread of the disease. For example, in the areas of high population density, the spread of COVID-19 tends to be fast. By zoning and adding multi-level physical borders, the transmission of the disease may be effectively reduced.

When the recovery from infection is considered, the system is no longer going to reach the outbreak endpoint. In contrast, depending on the infection rate, recovery speed and rate, and many other factors, two possible states may happen: (1) "zero-COVID" state, which means all the infected people have recovered; (2) "living with COVID" state, during which the system reaching an equilibrium state and the number of newly-infected people is equal to the number of newly-recovered people. In this section, the influential factors on the transition between "zero-COVID" state and "living with COVID" state are investigated.

First, the influence of moving distance determined by parameter $c_{1}$ is studied. When individuals tend to take long movements in same-level containers (i.e. mobility parameter $c_{1}$ is smaller than 0.25), the system is able to achieve "zero-COVID". Increasing $c_{1}$ leads to the local accumulation of people and facilitates the spread of the disease. Bifurcation happens when $c_{1}$ becomes larger than 0.25. The transition appears to be discontinuous, and a hysteresis loop appears in the bi-stable region where both "zero-COVID" state (black dots) and "living with COVID" state (orange dots) coexist. When $c_{1}$ is larger than 0.4, the percentage of infected people in the population increases to more than 30% (Figure 6a). Figure 6b shows that the probability of zero-COVID after a long enough time rapidly decreases to zero when $c_{1}$ is larger than 0.25. We further show the variation of the infected people ratio over time (Figure 6c). If the system eventually reaches a "zero-COVID" state, the fluctuation curves are labeled in pink colors. Otherwise, the curves are labeled in grey colors. Figure 6c shows that there is a peak of the infected ratio at the beginning, because healthy people keeps being infected but infected people have not recovered yet. Afterwards, the infected ratio decreases and then fluctuates between 3% to 5%. In some cases, the infected ratio decreases to zero, indicating that all the sick people become recovered.

Next, effects of the percentage of low mobility people in the population is evaluated. Increasing the percentage of low mobility people reduces the spread of the disease and the infected ratio. Thus, the system transitions from the "living with COVID" state to a bistable region, and then to the "zero-COVID" state (Figure 6d). The probability of achieving zero-COVID increases rapidly when the low mobility people ratio is larger than 60% (Figure 6e). In addition, the variation of the infected ratio is large and the time for the system to reach the "zero-COVID" state is longer when the low mobility people ratio is 70% (Figure 6f). This phenomenon may be caused by the stochastic locations and movements of infected people. When the infected people stay in the crowds for a while, it facilitates the disease transmission and increases the infected ratio. But if the infected people are located in a remote area, the low mobility characteristics helps reducing the infection.

Lastly, we also found that varying the population number leads to the bifurcation, and the bistable region with both "zero-COVID" state and "living with COVID" state appears at an intermediate population number (i.e. around 400 to 500 people) (Figure 6g). The probability of zero COVID decreases in the bistable region (Figure 6h). When the population size is larger than 600, the system is likely to reach an equilibrium and has more than 20% of the population got infected (Figure 6g). When the population number is less than 400, all the infected people can recover from the disease at the first 1000 time steps (Figure 6i) in most simulation cases.

In addition, effects of infection and recovery probabilities on the disease spread and the infected ratio have also been studied. As expected, increasing the infection probability leads to the increase of infected ratio, while increasing the recovery probability reduces infected ratio. These results are trivial and are not shown here.

When people take frequent and short local movements, or the population size is large, the spread of the disease is easier and faster. The infected people ratio is larger than 10%. Reducing local crowdedness helps decrease the infected ratio to around 5%. Then, fluctuations in the system originated from the random human travel, infection and recovery behaviors may lead to the disappearing of COVID-19 infection (Figures 6c, f and i). The probability to zero COVID-19 and the time it takes are influenced by many factors. In this section, we analyzed the phase diagrams of the infected ratio over time influenced by mobility parameter $c_{1}$, low mobility people ratio and population number.

Figure 7a shows that the system reaches equilibrium state quickly (i.e. either the "zero-COVID" state or the "living with COVID" state) when $c_{1}$ is larger than 0.4 or smaller than 0.25. When $c_{1}\in[0.25,0.4]$, the system is bistable and the probability of zero-COVID drops significantly with a small increase of $c_{1}$. The time for a thorough recovery from COVID-19 becomes longer. From Figure 7b, it can be found that having a large portion of low mobility people in the population (i.e. more than 30%) increases the time needed for zero-COVID. Though the recovery time is longer, the probability of reaching "zero-COVID" state is higher. When the low mobility people ratio is small, the system is likely to be in the "living with COVID" equilibrium state. Moreover, the population number also influences the transient COVID-19 transmission dynamics (Figure 7c). When the number of people in the population is small (i.e. less than 400), all the people can recover from the disease within the first 500 time steps. When the people number increases, the probability of reaching the zero-COVID state is lower and the population recovery takes a longer time. When the population contains more than 800 people, the newly infected people and recovered people achieves a dynamical equilibrium. The system is in the "live with COVID" state, and more than 20% of people in the population are infected.